source: trunk/yat/utility/SVD.h @ 1797

Last change on this file since 1797 was 1797, checked in by Peter, 12 years ago

updating copyright statements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.6 KB
Line 
1#ifndef _theplu_yat_utility_svd_
2#define _theplu_yat_utility_svd_
3
4// $Id: SVD.h 1797 2009-02-12 18:07:10Z peter $
5
6/*
7  Copyright (C) 2003 Daniel Dalevi, Jari Häkkinen
8  Copyright (C) 2004 Jari Häkkinen
9  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
10  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
11  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
12
13  This file is part of the yat library, http://dev.thep.lu.se/yat
14
15  The yat library is free software; you can redistribute it and/or
16  modify it under the terms of the GNU General Public License as
17  published by the Free Software Foundation; either version 3 of the
18  License, or (at your option) any later version.
19
20  The yat library is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
23  General Public License for more details.
24
25  You should have received a copy of the GNU General Public License
26  along with yat. If not, see <http://www.gnu.org/licenses/>.
27*/
28
29#include "Matrix.h"
30#include "Vector.h"
31
32#include <gsl/gsl_linalg.h>
33
34namespace theplu {
35namespace yat {
36namespace utility {
37
38  class VectorBase;
39
40  /**
41     @brief Singular Value Decomposition
42
43     Class encapsulating GSL methods for singular value
44     decomposition, SVD.
45
46     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
47
48     A = Matrix to be decomposed, size MxN\n
49     U = Orthogonal matrix, size MxN\n
50     S = Diagonal matrix of singular values, size NxN\n
51     V = Orthogonal matrix, size NxN\n
52  */
53  class SVD
54  {
55  public:
56
57    /**
58       A number of SVD algorithms are implemented in GSL. They have
59       their strengths and weaknesses.
60   
61       \see GSL's SVD documentation.
62    */
63    enum SVDalgorithm {
64      GolubReinsch,
65      ModifiedGolubReinsch,
66      Jacobi
67    };
68
69    /**
70       \brief Constructs an SVD object using the matrix Ain as only
71       input. The input matrix is copied for further use in the
72       object.
73
74       \note Number of rows must be equal or larger than number of columns.
75    */
76    SVD(const utility::Matrix& Ain);
77
78    /**
79       \brief The destructor
80    */
81    ~SVD(void);
82
83    /**
84       \brief This function will perform SVD with the method specified
85       by \a algo.
86
87       \throw GSL_error if the underlying GSL function fails.
88    */
89    void decompose(SVDalgorithm algo=GolubReinsch);
90
91    /**
92       \brief Access to the s vector.
93
94       \return A copy of the s vector.
95
96       \note If decompose() has not been run the outcome of the call
97       is undefined.
98    */
99    const utility::Vector& s(void) const;
100
101    /**
102       \brief Solve the system \f$ Ax=b \f$ using the decomposition of
103       A.
104
105       \note If decompose() has not been run the outcome of the call
106       is undefined.
107
108       \throw GSL_error if the underlying GSL function fails.
109    */
110    void solve(const utility::VectorBase& b, utility::Vector& x);
111
112    /**
113       \brief Access to the U matrix.
114
115       \return A copy of the U matrix.
116
117       \note If decompose() has not been run the outcome of the call
118       is undefined.
119    */
120    const utility::Matrix& U(void) const;
121
122    /**
123       \brief Access to the V matrix.
124
125       \return A copy of the V matrix.
126
127       \note If decompose() has not been run the outcome of the call
128       is undefined.
129    */
130    const utility::Matrix& V(void) const;
131
132  private:
133    /**
134       \brief Call GSL's Jacobi algorithm for SVD.
135
136       \return gsl_error status.
137    */
138    int jacobi(void);
139
140    /**
141       \brief Call GSL's Golub-Reinsch algorithm for SVD.
142
143       \return gsl_error status.
144    */
145    int golub_reinsch(void);
146
147    /**
148       \brief Call GSL's modified Golub-Reinch algorithm for SVD.
149
150       \return gsl_error status.
151    */
152    int modified_golub_reinsch(void);
153
154    utility::Matrix U_, V_;
155    utility::Vector s_;
156  }; 
157
158}}} // of namespace utility, yat, and theplu
159
160#endif
Note: See TracBrowser for help on using the repository browser.