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

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

docs

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