1 | // $Id: SVD.h 43 2004-02-26 15:10:19Z jari $ |
---|
2 | |
---|
3 | #ifndef _theplu_cpptools_svd_ |
---|
4 | #define _theplu_cpptools_svd_ |
---|
5 | |
---|
6 | #include "matrix.h" |
---|
7 | #include "vector.h" |
---|
8 | |
---|
9 | namespace theplu { |
---|
10 | // Jari, this should probably go somewhere else for doxygen to process. |
---|
11 | /// |
---|
12 | /// All C++ tools functionality is to be defined within cpptools |
---|
13 | /// namespace. |
---|
14 | /// |
---|
15 | namespace cpptools { |
---|
16 | |
---|
17 | /** |
---|
18 | Class encapsulating methods for Singular Value Decomposition. |
---|
19 | \n\n |
---|
20 | A = Matrix to be decomposed, size MxN\n |
---|
21 | U = Orthogonal matrix, size MxN\n |
---|
22 | S = Diagonal matrix of singular values, size NxN\n |
---|
23 | V = Orthogonal matrix, size NxN\n |
---|
24 | A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n |
---|
25 | |
---|
26 | Shouldn't there be a reference to the algortithm here? |
---|
27 | |
---|
28 | @version 1.0 $Revision: 43 $ $Date: 2004-02-26 15:10:19 +0000 (Thu, 26 Feb 2004) $ |
---|
29 | */ |
---|
30 | |
---|
31 | class SVD |
---|
32 | { |
---|
33 | static const double MAXTOL = 1E-6; |
---|
34 | |
---|
35 | enum SVDtypes { |
---|
36 | Unmodified, |
---|
37 | Modified, |
---|
38 | Jacobi, |
---|
39 | Solver |
---|
40 | }; |
---|
41 | |
---|
42 | public: |
---|
43 | /** |
---|
44 | The default constructor should only be used for testing!!! |
---|
45 | */ |
---|
46 | SVD(void); |
---|
47 | |
---|
48 | |
---|
49 | /** |
---|
50 | Constructs an SVD object using the matrix A as only |
---|
51 | input. The input matrix is copied for further use in the |
---|
52 | object. |
---|
53 | */ |
---|
54 | SVD( const gslapi::matrix& ); |
---|
55 | |
---|
56 | /** |
---|
57 | Constructor initializing the SVD object for Solver. |
---|
58 | Solver requires an extra-vector paramter \a b: \f$Ax = b\f$ |
---|
59 | The \f$x\f$ vector will be reached using the get_x() |
---|
60 | function. |
---|
61 | @see get_x() function. |
---|
62 | */ |
---|
63 | |
---|
64 | SVD( const gslapi::matrix&, const gslapi::vector& b); |
---|
65 | ~SVD(); |
---|
66 | |
---|
67 | /** |
---|
68 | This function will run the method specified by the SVDtypes. Precently |
---|
69 | there are 3 SVD methods implemented and one solver. These are: \n\n |
---|
70 | Unmodified: (see GSL documentation about gsl_linalg_SV_decomp function)\n |
---|
71 | Modified: This method is faster when M>>N. (see GSL documentation |
---|
72 | about gsl_linalg_SV_decomp_mod function)\n |
---|
73 | Jacobi: Computes singular values to higer accuracy than default method |
---|
74 | (see GSL documentation about gsl_linalg_SV_decomp_jacobi.\n |
---|
75 | Solver: Using Unmodified to perform SVD and then solves the system \f$Ax=b\f$ |
---|
76 | */ |
---|
77 | bool process( SVDtypes stype = Unmodified ); |
---|
78 | |
---|
79 | /** |
---|
80 | Function will return the matrix U. Require that process() has been |
---|
81 | executed. |
---|
82 | */ |
---|
83 | gslapi::matrix get_u() const { return A_; } |
---|
84 | |
---|
85 | /** |
---|
86 | Function will return the matrix V. Require that process() has been |
---|
87 | executed. |
---|
88 | */ |
---|
89 | gslapi::matrix get_v() const { return V_; } |
---|
90 | |
---|
91 | |
---|
92 | /** |
---|
93 | Function will return the vector x containing the solution |
---|
94 | \f$ x=A^{-1}b \f$. Require that process( Solver ) has been |
---|
95 | executed. |
---|
96 | @see Read about Solver in the bool process() function. |
---|
97 | */ |
---|
98 | |
---|
99 | gslapi::vector get_x() const { return x_; } |
---|
100 | |
---|
101 | /** |
---|
102 | Function returning diagonal matrix S. Require that process() has been |
---|
103 | executed. |
---|
104 | */ |
---|
105 | gslapi::matrix get_s() const; |
---|
106 | |
---|
107 | |
---|
108 | /** |
---|
109 | Function returning diagonal matrix S in vector form. |
---|
110 | */ |
---|
111 | gslapi::vector get_s_vec() const { return s_; } |
---|
112 | |
---|
113 | |
---|
114 | /** |
---|
115 | This method will execute Default, Unmodified and Jacobi methods one at a time |
---|
116 | and calculate the error, \f$ e = \left\Vert A - USV^{T} \right\Vert_{2} \f$. |
---|
117 | The method will return "true" if \f$e < 10^{-6}\f$ else false. A test program |
---|
118 | is available that executes this method, c++_tools/test/svd_test. No test is |
---|
119 | performed for solver but the aim is to add one in the future. |
---|
120 | */ |
---|
121 | bool test(); |
---|
122 | |
---|
123 | private: |
---|
124 | |
---|
125 | gslapi::matrix A_, V_; |
---|
126 | gslapi::vector s_, b_, x_; |
---|
127 | bool loaded_, solver_; |
---|
128 | |
---|
129 | //SVD methods |
---|
130 | bool process_default(); |
---|
131 | bool process_modified(); |
---|
132 | bool process_jacobi(); |
---|
133 | bool process_solver(); |
---|
134 | |
---|
135 | double norm(gslapi::matrix& A) const; // move into class gslapi::matrix? |
---|
136 | }; |
---|
137 | |
---|
138 | }} // of namespace cpptools and namespace theplu |
---|
139 | |
---|
140 | #endif |
---|