source: trunk/yat/utility/SVD.h

Last change on this file was 3938, checked in by Peter, 11 months ago

let configure fail if compiler is not a C++11 compiler. YAT_HAVE_RVALUE and friends are no longer defined in 'config.h' and guards araound declarations and implementations are removed. refs #949

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.8 KB
Line 
1#ifndef _theplu_yat_utility_svd_
2#define _theplu_yat_utility_svd_
3
4// $Id: SVD.h 3938 2020-07-16 13:16:56Z 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  Copyright (C) 2018 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 "config_public.h"
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.
63   
64       \see GSL's SVD documentation.
65    */
66    enum SVDalgorithm {
67      GolubReinsch,
68      ModifiedGolubReinsch,
69      Jacobi
70    };
71
72    /**
73       \brief Constructs an SVD object using the matrix Ain as only
74       input. The input matrix is copied for further use in the
75       object.
76
77       \note Number of rows must be equal or larger than number of columns.
78    */
79    explicit SVD(const Matrix& Ain);
80
81    /**
82       Same behaviour as const& constructor, but \a Ain is moved
83       rather than copied.
84
85       \since new in yat 0.16
86     */
87    explicit SVD(Matrix&& Ain);
88
89    /**
90       \brief The destructor
91    */
92    ~SVD(void);
93
94    /**
95       \brief This function will perform SVD with the method specified
96       by \a algo.
97
98       \throw GSL_error if the underlying GSL function fails.
99    */
100    void decompose(SVDalgorithm algo=GolubReinsch);
101
102    /**
103       \brief Access to the s vector.
104
105       \return A copy of the s vector.
106
107       \note If decompose() has not been run the outcome of the call
108       is undefined.
109    */
110    const Vector& s(void) const;
111
112    /**
113       \brief Solve the system \f$ Ax=b \f$ using the decomposition of
114       A.
115
116       \note If decompose() has not been run the outcome of the call
117       is undefined.
118
119       \throw GSL_error if the underlying GSL function fails.
120    */
121    void solve(const VectorBase& b, Vector& x);
122
123    /**
124       \brief Access to the U matrix.
125
126       \return A copy of the U matrix.
127
128       \note If decompose() has not been run the outcome of the call
129       is undefined.
130    */
131    const Matrix& U(void) const;
132
133    /**
134       \brief Access to the V matrix.
135
136       \return A copy of the V matrix.
137
138       \note If decompose() has not been run the outcome of the call
139       is undefined.
140    */
141    const Matrix& V(void) const;
142
143  private:
144    /**
145       \brief Call GSL's Jacobi algorithm for SVD.
146
147       \return gsl_error status.
148    */
149    int jacobi(void);
150
151    /**
152       \brief Call GSL's Golub-Reinsch algorithm for SVD.
153
154       \return gsl_error status.
155    */
156    int golub_reinsch(void);
157
158    /**
159       \brief Call GSL's modified Golub-Reinch algorithm for SVD.
160
161       \return gsl_error status.
162    */
163    int modified_golub_reinsch(void);
164
165    Matrix U_, V_;
166    Vector s_;
167  };
168
169}}} // of namespace utility, yat, and theplu
170
171#endif
Note: See TracBrowser for help on using the repository browser.