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

Last change on this file since 1437 was 1437, checked in by Peter, 13 years ago

merge patch release 0.4.2 to trunk. Delta 0.4.2-0.4.1

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.0 KB
Line 
1#ifndef _theplu_yat_utility_svd_
2#define _theplu_yat_utility_svd_
3
4// $Id: SVD.h 1437 2008-08-25 17:55:00Z 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 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 2 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 this program; if not, write to the Free Software
28  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
29  02111-1307, USA.
30*/
31
32#include "Matrix.h"
33#include "Vector.h"
34
35#include <gsl/gsl_linalg.h>
36
37namespace theplu {
38namespace yat {
39namespace utility {
40
41  class VectorBase;
42
43  /**
44     @brief Singular Value Decomposition
45
46     Class encapsulating GSL methods for singular value
47     decomposition, SVD.
48
49     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
50
51     A = Matrix to be decomposed, size MxN\n
52     U = Orthogonal matrix, size MxN\n
53     S = Diagonal matrix of singular values, size NxN\n
54     V = Orthogonal matrix, size NxN\n
55  */
56  class SVD
57  {
58  public:
59
60    /**
61       A number of SVD algorithms are implemented in GSL. They have
62       their strengths and weaknesses, check the GSL documentation.
63   
64       There are restrictions on the matrix dimensions depending on
65       which SVD algorithm is used. From the GSL's SVD source code one
66       finds that the Golub-Reinsch algorithm implementation will not
67       work on matrices with fewer rows than columns, the same is also
68       true for the modified Golub-Reinsch algorithm.
69   
70       \see GSL's SVD documentation.
71    */
72    enum SVDalgorithm {
73      GolubReinsch,
74      ModifiedGolubReinsch,
75      Jacobi
76    };
77
78    /**
79       \brief Constructs an SVD object using the matrix Ain as only
80       input. The input matrix is copied for further use in the
81       object.
82    */
83    SVD(const utility::Matrix& Ain);
84
85    /**
86       \brief The destructor
87    */
88    ~SVD(void);
89
90    /**
91       \brief This function will perform SVD with the method specified
92       by \a algo.
93
94       \throw GSL_error if the underlying GSL function fails.
95    */
96    void decompose(SVDalgorithm algo=GolubReinsch);
97
98    /**
99       \brief Access to the s vector.
100
101       \return A copy of the s vector.
102
103       \note If decompose() has not been run the outcome of the call
104       is undefined.
105    */
106    const utility::Vector& s(void) const;
107
108    /**
109       \brief Solve the system \f$ Ax=b \f$ using the decomposition of
110       A.
111
112       \note If decompose() has not been run the outcome of the call
113       is undefined.
114
115       \throw GSL_error if the underlying GSL function fails.
116    */
117    void solve(const utility::VectorBase& b, utility::Vector& x);
118
119    /**
120       \brief Access to the U matrix.
121
122       \return A copy of the U matrix.
123
124       \note If decompose() has not been run the outcome of the call
125       is undefined.
126    */
127    const utility::Matrix& U(void) const;
128
129    /**
130       \brief Access to the V matrix.
131
132       \return A copy of the V matrix.
133
134       \note If decompose() has not been run the outcome of the call
135       is undefined.
136    */
137    const utility::Matrix& V(void) const;
138
139  private:
140    /**
141       \brief Call GSL's Jacobi algorithm for SVD.
142
143       \return gsl_error status.
144    */
145    int jacobi(void);
146
147    /**
148       \brief Call GSL's Golub-Reinsch algorithm for SVD.
149
150       \return gsl_error status.
151    */
152    int golub_reinsch(void);
153
154    /**
155       \brief Call GSL's modified Golub-Reinch algorithm for SVD.
156
157       \return gsl_error status.
158    */
159    int modified_golub_reinsch(void);
160
161    utility::Matrix U_, V_;
162    utility::Vector s_;
163  }; 
164
165}}} // of namespace utility, yat, and theplu
166
167#endif
Note: See TracBrowser for help on using the repository browser.