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

Last change on this file since 1703 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

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