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

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

Addresses #2. Continued adding GSL_error exceptions.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1#ifndef _theplu_yat_utility_svd_
2#define _theplu_yat_utility_svd_
3
4// $Id: SVD.h 751 2007-02-17 12:25:36Z jari $
5
6/*
7  Copyright (C) The authors contributing to this file.
8
9  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "matrix.h"
28#include "vector.h"
29
30#include <gsl/gsl_linalg.h>
31
32namespace theplu {
33namespace yat {
34namespace utility {
35
36  /**
37     Class encapsulating GSL methods for singular value decomposition,
38     SVD.
39
40     A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n   
41
42     A = Matrix to be decomposed, size MxN\n
43     U = Orthogonal matrix, size MxN\n
44     S = Diagonal matrix of singular values, size NxN\n
45     V = Orthogonal matrix, size NxN\n
46  */
47  class SVD
48  {
49  public:
50
51    /**
52       A number of SVD algorithms are implemented in GSL. They have
53       their strengths and weaknesses, check the GSL documentation.
54   
55       There are restrictions on the matrix dimensions depending on
56       which SVD algorithm is used. From the GSL's SVD source code one
57       finds that the Golub-Reinsch algorithm implementation will not
58       work on matrices with fewer rows than columns, the same is also
59       true for the modified Golub-Reinsch algorithm.
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    SVD(const utility::matrix& Ain);
75
76    /**
77       \brief The destructor
78    */
79    ~SVD(void);
80
81    /**
82       \brief This function will perform SVD with the method specified
83       by \a algo.
84
85       \throw GSL_error if the underlying GSL function fails.
86    */
87    void decompose(SVDalgorithm algo=GolubReinsch);
88
89    /**
90       \brief Access to the s vector.
91
92       \return A copy of the s vector.
93
94       \note If decompose() has not been run the outcome of the call
95       is undefined.
96    */
97    const utility::vector& s(void) const;
98
99    /**
100       \brief Solve the system \f$ Ax=b \f$ using the decomposition of
101       A.
102
103       \note If decompose() has not been run the outcome of the call
104       is undefined.
105
106       \throw GSL_error if the underlying GSL function fails.
107    */
108    void solve(const utility::vector& b, utility::vector& x);
109
110    /**
111       \brief Access to the U matrix.
112
113       \return A copy of the U matrix.
114
115       \note If decompose() has not been run the outcome of the call
116       is undefined.
117    */
118    const utility::matrix& U(void) const;
119
120    /**
121       \brief Access to the V matrix.
122
123       \return A copy of the V matrix.
124
125       \note If decompose() has not been run the outcome of the call
126       is undefined.
127    */
128    const utility::matrix& V(void) const;
129
130  private:
131    /**
132       \brief Call GSL's Jacobi algorithm for SVD.
133
134       \return gsl_error status.
135    */
136    int jacobi(void);
137
138    /**
139       \brief Call GSL's Golub-Reinsch algorithm for SVD.
140
141       \return gsl_error status.
142    */
143    int golub_reinsch(void);
144
145    /**
146       \brief Call GSL's modified Golub-Reinch algorithm for SVD.
147
148       \return gsl_error status.
149    */
150    int modified_golub_reinsch(void);
151
152    utility::matrix U_, V_;
153    utility::vector s_;
154  }; 
155
156}}} // of namespace utility, yat, and theplu
157
158#endif
Note: See TracBrowser for help on using the repository browser.