# source:branches/better_matrix_class/lib/utility/SVD.h@415

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

Removed some gslapi copy intensive operators.
Made gslapi assignment operators ignore views, and change them
into normal vectors if assigned.
Cleaning up of code aswell.

• 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 415 2005-12-01 15:52:36Z jari \$
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    inline SVD(const gslapi::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    ///
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 gslapi::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(gslapi::vector b, gslapi::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    ///
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 gslapi::matrix& U(void) const { return U_; }
99
100    ///
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 gslapi::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    gslapi::matrix U_, V_;
118    gslapi::vector s_;
119  };
120
121}} // of namespace utility and namespace theplu
122
123#endif
Note: See TracBrowser for help on using the repository browser.