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

Last change on this file since 865 was 865, checked in by Peter, 14 years ago

changing URL to http://trac.thep.lu.se/trac/yat

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