source: trunk/c++_tools/utility/SVD.h @ 616

Last change on this file since 616 was 616, checked in by Jari Häkkinen, 15 years ago

Removed gslapi namespace and put the code into utility namespace.
Move #ifndef _header_ idiom to top of touched header files.
Removed unneccesary #includes, and added needed #includes.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.1 KB
Line 
1#ifndef _theplu_utility_svd_
2#define _theplu_utility_svd_
3
4// $Id: SVD.h 616 2006-08-31 08:52:02Z jari $
5
6#include <c++_tools/utility/matrix.h>
7#include <c++_tools/utility/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    inline SVD(const utility::matrix& Ain)
55      : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns()) {}
56
57    inline ~SVD(void) {}
58
59    ///
60    /// This function will perform SVD with the method specified by \a
61    /// algo.
62    ///
63    /// @return Whatever GSL returns.
64    ///
65    int decompose(SVDalgorithm algo=GolubReinsch);
66
67    ///
68    /// Access to the s vector.
69    ///
70    /// @return A copy of the s vector.
71    ///
72    /// @note If decompose() has not been run the outcome of the call
73    /// is undefined.
74    ///
75    inline const utility::vector& s(void) const { return s_; }
76
77    ///
78    /// Solve the system \f$ Ax=b \f$ using the decomposition of A.
79    ///
80    /// @note If decompose() has not been run the outcome of the call
81    /// is undefined.
82    ///
83    /// @return Whatever GSL returns.
84    ///
85    inline int solve(utility::vector b, utility::vector x)
86      { return gsl_linalg_SV_solve(U_.gsl_matrix_p(), V_.gsl_matrix_p(), 
87                                   s_.gsl_vector_p(), b.gsl_vector_p(),
88                                   x.gsl_vector_p()); }
89
90    ///
91    /// Access to the U matrix.
92    ///
93    /// @return A copy of the U matrix.
94    ///
95    /// @note If decompose() has not been run the outcome of the call
96    /// is undefined.
97    ///
98    inline const utility::matrix& U(void) const { return U_; }
99
100    ///
101    /// Access to the V matrix.
102    ///
103    /// @return A copy of the V matrix.
104    ///
105    /// @note If decompose() has not been run the outcome of the call
106    /// is undefined.
107    ///
108    inline const utility::matrix& V(void) const { return V_; }
109
110  private:
111    inline int jacobi(void)
112      { return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_p(), V_.gsl_matrix_p(), 
113                                           s_.gsl_vector_p()); }
114    int golub_reinsch(void);
115    int modified_golub_reinsch(void);
116
117    utility::matrix U_, V_;
118    utility::vector s_;
119  }; 
120
121}} // of namespace utility and namespace theplu
122
123#endif
Note: See TracBrowser for help on using the repository browser.