[680] | 1 | #ifndef _theplu_yat_utility_svd_ |
---|
| 2 | #define _theplu_yat_utility_svd_ |
---|
[4] | 3 | |
---|
[616] | 4 | // $Id: SVD.h 1487 2008-09-10 08:41:36Z jari $ |
---|
[4] | 5 | |
---|
[675] | 6 | /* |
---|
[831] | 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 |
---|
[1275] | 12 | Copyright (C) 2008 Peter Johansson |
---|
[616] | 13 | |
---|
[1437] | 14 | This file is part of the yat library, http://dev.thep.lu.se/yat |
---|
[675] | 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 |
---|
[1486] | 18 | published by the Free Software Foundation; either version 3 of the |
---|
[675] | 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 |
---|
[1487] | 27 | along with yat. If not, see <http://www.gnu.org/licenses/>. |
---|
[675] | 28 | */ |
---|
| 29 | |
---|
[1121] | 30 | #include "Matrix.h" |
---|
[1120] | 31 | #include "Vector.h" |
---|
[675] | 32 | |
---|
[227] | 33 | #include <gsl/gsl_linalg.h> |
---|
| 34 | |
---|
[42] | 35 | namespace theplu { |
---|
[680] | 36 | namespace yat { |
---|
[301] | 37 | namespace utility { |
---|
[4] | 38 | |
---|
[1016] | 39 | class VectorBase; |
---|
| 40 | |
---|
[227] | 41 | /** |
---|
[767] | 42 | @brief Singular Value Decomposition |
---|
[227] | 43 | |
---|
[767] | 44 | Class encapsulating GSL methods for singular value |
---|
| 45 | decomposition, SVD. |
---|
| 46 | |
---|
[1478] | 47 | A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n |
---|
[227] | 48 | |
---|
[7] | 49 | A = Matrix to be decomposed, size MxN\n |
---|
[1478] | 50 | U = Orthogonal matrix, size MxN\n |
---|
| 51 | S = Diagonal matrix of singular values, size NxN\n |
---|
[7] | 52 | V = Orthogonal matrix, size NxN\n |
---|
| 53 | */ |
---|
[4] | 54 | class SVD |
---|
[42] | 55 | { |
---|
[227] | 56 | public: |
---|
[4] | 57 | |
---|
[751] | 58 | /** |
---|
| 59 | A number of SVD algorithms are implemented in GSL. They have |
---|
[1478] | 60 | their strengths and weaknesses. |
---|
[751] | 61 | |
---|
[1478] | 62 | \see GSL's SVD documentation. |
---|
[751] | 63 | */ |
---|
[227] | 64 | enum SVDalgorithm { |
---|
| 65 | GolubReinsch, |
---|
| 66 | ModifiedGolubReinsch, |
---|
| 67 | Jacobi |
---|
[42] | 68 | }; |
---|
[4] | 69 | |
---|
[751] | 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. |
---|
[1478] | 74 | |
---|
| 75 | \note Number of rows must be equal or larger than number of columns. |
---|
[751] | 76 | */ |
---|
[1121] | 77 | SVD(const utility::Matrix& Ain); |
---|
[11] | 78 | |
---|
[751] | 79 | /** |
---|
| 80 | \brief The destructor |
---|
| 81 | */ |
---|
[703] | 82 | ~SVD(void); |
---|
[7] | 83 | |
---|
[751] | 84 | /** |
---|
| 85 | \brief This function will perform SVD with the method specified |
---|
| 86 | by \a algo. |
---|
[7] | 87 | |
---|
[751] | 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 | */ |
---|
[1120] | 100 | const utility::Vector& s(void) const; |
---|
[7] | 101 | |
---|
[751] | 102 | /** |
---|
| 103 | \brief Solve the system \f$ Ax=b \f$ using the decomposition of |
---|
| 104 | A. |
---|
[7] | 105 | |
---|
[751] | 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 | */ |
---|
[1120] | 111 | void solve(const utility::VectorBase& b, utility::Vector& x); |
---|
[751] | 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 | */ |
---|
[1121] | 121 | const utility::Matrix& U(void) const; |
---|
[7] | 122 | |
---|
[751] | 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 | */ |
---|
[1121] | 131 | const utility::Matrix& V(void) const; |
---|
[7] | 132 | |
---|
[4] | 133 | private: |
---|
[751] | 134 | /** |
---|
| 135 | \brief Call GSL's Jacobi algorithm for SVD. |
---|
| 136 | |
---|
| 137 | \return gsl_error status. |
---|
| 138 | */ |
---|
[715] | 139 | int jacobi(void); |
---|
[751] | 140 | |
---|
| 141 | /** |
---|
| 142 | \brief Call GSL's Golub-Reinsch algorithm for SVD. |
---|
| 143 | |
---|
| 144 | \return gsl_error status. |
---|
| 145 | */ |
---|
[227] | 146 | int golub_reinsch(void); |
---|
[751] | 147 | |
---|
| 148 | /** |
---|
| 149 | \brief Call GSL's modified Golub-Reinch algorithm for SVD. |
---|
| 150 | |
---|
| 151 | \return gsl_error status. |
---|
| 152 | */ |
---|
[227] | 153 | int modified_golub_reinsch(void); |
---|
[7] | 154 | |
---|
[1121] | 155 | utility::Matrix U_, V_; |
---|
[1120] | 156 | utility::Vector s_; |
---|
[4] | 157 | }; |
---|
| 158 | |
---|
[687] | 159 | }}} // of namespace utility, yat, and theplu |
---|
[42] | 160 | |
---|
[4] | 161 | #endif |
---|