source: trunk/lib/utility/SVD.h @ 301

Last change on this file since 301 was 301, checked in by Peter, 17 years ago

modified includes in tests

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.2 KB
Line 
1// $Id: SVD.h 301 2005-04-30 13:39:27Z peter $
2
3#ifndef _theplu_cpptools_svd_
4#define _theplu_cpptools_svd_
5
6#include <c++_tools/gslapi/matrix.h>
7#include <c++_tools/gslapi/vector.h>
8
9#include <gsl/gsl_linalg.h>
10
11namespace theplu {
12// Jari, this should probably go somewhere else for doxygen to process.
13///
14/// All C++ tools functionality is to be defined within cpptools
15/// namespace.
16///
17namespace utility {
18
19  // Jari check that interface is complete
20
21  /**
22     Class encapsulating GSL methods for singular value decomposition,
23     SVD.
24
25     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
26
27     A = Matrix to be decomposed, size MxN\n
28     U = Orthogonal matrix, size MxN\n
29     S = Diagonal matrix of singular values, size NxN\n
30     V = Orthogonal matrix, size NxN\n
31  */
32
33  class SVD
34  {
35  public:
36
37    ///
38    /// A number of SVD algorithms are implemented in GSL. They have
39    /// their strengths and weaknesses, check the GSL documentation.
40    ///
41    /// There are restrictions on the matrix dimensions depending on
42    /// which SVD algorithm is used. From the GSL's SVD source code
43    /// one finds that the Golub-Reinsch algorithm implementation will
44    /// not work on matrices with fewer rows than columns, the same is
45    /// also true for the modified Golub-Reinsch algorithm.
46    ///
47    /// @see GSL's SVD documentation.
48    ///
49    enum SVDalgorithm {
50      GolubReinsch,
51      ModifiedGolubReinsch,
52      Jacobi
53    };
54
55    ///
56    /// Constructs an SVD object using the matrix A as only input. The
57    /// input matrix is copied for further use in the object.
58    ///
59    SVD(const gslapi::matrix& );
60
61    ~SVD(void);
62
63    ///
64    /// This function will perform SVD with the method specified by \a
65    /// algo.
66    ///
67    /// @return Whatever GSL returns.
68    ///
69    int decompose(SVDalgorithm algo=GolubReinsch);
70
71    ///
72    /// Access to the s vector.
73    ///
74    /// @return A copy of the s vector.
75    ///
76    /// @note If decompose() has not been run the outcome of the call
77    /// is undefined.
78    ///
79    gslapi::vector s(void) const { return s_; }
80
81    ///
82    /// Solve the system \f$Ax=b\f$ using the decomposition of A.
83    ///
84    /// @note If decompose() has not been run the outcome of the call
85    /// is undefined.
86    ///
87    /// @return Whatever GSL returns.
88    ///
89    inline int solve(gslapi::vector b, gslapi::vector x)
90      { return gsl_linalg_SV_solve(U_.gsl_matrix_pointer(),
91                                   V_.gsl_matrix_pointer(), 
92                                   s_.gsl_vector_pointer(),
93                                   b.gsl_vector_pointer(),
94                                   x.TEMP_gsl_vector_pointer());
95      }
96
97    ///
98    /// Access to the U matrix.
99    ///
100    /// @return A copy of the U matrix.
101    ///
102    /// @note If decompose() has not been run the outcome of the call
103    /// is undefined.
104    ///
105    gslapi::matrix U(void) const { return U_; }
106
107    ///
108    /// Access to the V matrix.
109    ///
110    /// @return A copy of the V matrix.
111    ///
112    /// @note If decompose() has not been run the outcome of the call
113    /// is undefined.
114    ///
115    gslapi::matrix V(void) const { return V_; }
116
117  private:
118    inline int jacobi(void)
119      { return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_pointer(),
120                                           V_.gsl_matrix_pointer(), 
121                                           s_.TEMP_gsl_vector_pointer());
122      }
123    int golub_reinsch(void);
124    int modified_golub_reinsch(void);
125
126    gslapi::matrix U_, V_;
127    gslapi::vector s_;
128  }; 
129
130}} // of namespace utility and namespace theplu
131
132#endif
Note: See TracBrowser for help on using the repository browser.