yat/utility/SVD.h

Code
Comments
Other
Rev Date Author Line
680 11 Oct 06 jari 1 #ifndef _theplu_yat_utility_svd_
680 11 Oct 06 jari 2 #define _theplu_yat_utility_svd_
4 25 Feb 03 daniel 3
616 31 Aug 06 jari 4 // $Id$
4 25 Feb 03 daniel 5
675 10 Oct 06 jari 6 /*
2119 12 Dec 09 peter 7   Copyright (C) 2003 Daniel Dalevi, Jari Häkkinen
2119 12 Dec 09 peter 8   Copyright (C) 2004 Jari Häkkinen
2119 12 Dec 09 peter 9   Copyright (C) 2005 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 10   Copyright (C) 2006 Jari Häkkinen
2119 12 Dec 09 peter 11   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 12   Copyright (C) 2018, 2022, 2023 Peter Johansson
616 31 Aug 06 jari 13
1437 25 Aug 08 peter 14   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 15
675 10 Oct 06 jari 16   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 17   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 18   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 19   License, or (at your option) any later version.
675 10 Oct 06 jari 20
675 10 Oct 06 jari 21   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 22   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 23   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 24   General Public License for more details.
675 10 Oct 06 jari 25
675 10 Oct 06 jari 26   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 27   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 28 */
675 10 Oct 06 jari 29
3736 22 May 18 peter 30 #include "config_public.h"
3736 22 May 18 peter 31
1121 22 Feb 08 peter 32 #include "Matrix.h"
1120 21 Feb 08 peter 33 #include "Vector.h"
675 10 Oct 06 jari 34
227 01 Feb 05 jari 35 #include <gsl/gsl_linalg.h>
227 01 Feb 05 jari 36
42 26 Feb 04 jari 37 namespace theplu {
680 11 Oct 06 jari 38 namespace yat {
301 30 Apr 05 peter 39 namespace utility {
4 25 Feb 03 daniel 40
1016 01 Feb 08 peter 41   class VectorBase;
1016 01 Feb 08 peter 42
227 01 Feb 05 jari 43   /**
767 22 Feb 07 peter 44      @brief Singular Value Decomposition
227 01 Feb 05 jari 45
767 22 Feb 07 peter 46      Class encapsulating GSL methods for singular value
767 22 Feb 07 peter 47      decomposition, SVD.
767 22 Feb 07 peter 48
4200 19 Aug 22 peter 49      A = U S V' = (MxN)(NxN)(NxN) = (MxN)\n
227 01 Feb 05 jari 50
7 28 Feb 03 daniel 51      A = Matrix to be decomposed, size MxN\n
1478 06 Sep 08 peter 52      U = Orthogonal matrix, size MxN\n
1478 06 Sep 08 peter 53      S = Diagonal matrix of singular values, size NxN\n
7 28 Feb 03 daniel 54      V = Orthogonal matrix, size NxN\n
7 28 Feb 03 daniel 55   */
4 25 Feb 03 daniel 56   class SVD
42 26 Feb 04 jari 57   {
227 01 Feb 05 jari 58   public:
4 25 Feb 03 daniel 59
751 17 Feb 07 jari 60     /**
751 17 Feb 07 jari 61        A number of SVD algorithms are implemented in GSL. They have
1478 06 Sep 08 peter 62        their strengths and weaknesses.
4200 19 Aug 22 peter 63
4200 19 Aug 22 peter 64        \see GSL's SVD documentation.
751 17 Feb 07 jari 65     */
227 01 Feb 05 jari 66     enum SVDalgorithm {
227 01 Feb 05 jari 67        GolubReinsch,
227 01 Feb 05 jari 68       ModifiedGolubReinsch,
227 01 Feb 05 jari 69       Jacobi
42 26 Feb 04 jari 70     };
4 25 Feb 03 daniel 71
751 17 Feb 07 jari 72     /**
751 17 Feb 07 jari 73        \brief Constructs an SVD object using the matrix Ain as only
751 17 Feb 07 jari 74        input. The input matrix is copied for further use in the
751 17 Feb 07 jari 75        object.
1478 06 Sep 08 peter 76
1478 06 Sep 08 peter 77        \note Number of rows must be equal or larger than number of columns.
751 17 Feb 07 jari 78     */
4125 14 Jan 22 peter 79     explicit SVD(const MatrixBase& Ain);
11 03 Jun 03 daniel 80
751 17 Feb 07 jari 81     /**
3736 22 May 18 peter 82        Same behaviour as const& constructor, but \a Ain is moved
3736 22 May 18 peter 83        rather than copied.
3736 22 May 18 peter 84
3736 22 May 18 peter 85        \since new in yat 0.16
3736 22 May 18 peter 86      */
4125 14 Jan 22 peter 87     explicit SVD(MatrixMutable&& Ain);
3736 22 May 18 peter 88
3736 22 May 18 peter 89     /**
751 17 Feb 07 jari 90        \brief The destructor
751 17 Feb 07 jari 91     */
703 18 Dec 06 jari 92     ~SVD(void);
7 28 Feb 03 daniel 93
751 17 Feb 07 jari 94     /**
751 17 Feb 07 jari 95        \brief This function will perform SVD with the method specified
751 17 Feb 07 jari 96        by \a algo.
7 28 Feb 03 daniel 97
751 17 Feb 07 jari 98        \throw GSL_error if the underlying GSL function fails.
751 17 Feb 07 jari 99     */
751 17 Feb 07 jari 100     void decompose(SVDalgorithm algo=GolubReinsch);
751 17 Feb 07 jari 101
751 17 Feb 07 jari 102     /**
4286 30 Jan 23 peter 103        Calculate the inverse of the Matrix provided in the
4286 30 Jan 23 peter 104        constructor. Behaviour is undefined if the input matrix is not
4286 30 Jan 23 peter 105        square.
4286 30 Jan 23 peter 106
4286 30 Jan 23 peter 107        \since New in yat 0.21
4286 30 Jan 23 peter 108
4286 30 Jan 23 peter 109        \see inverse_svd(const MatrixBase&, Matrix&)
4286 30 Jan 23 peter 110      */
4286 30 Jan 23 peter 111     void inverse(Matrix& result) const;
4286 30 Jan 23 peter 112
4286 30 Jan 23 peter 113     /**
751 17 Feb 07 jari 114        \brief Access to the s vector.
751 17 Feb 07 jari 115
751 17 Feb 07 jari 116        \return A copy of the s vector.
751 17 Feb 07 jari 117
751 17 Feb 07 jari 118        \note If decompose() has not been run the outcome of the call
751 17 Feb 07 jari 119        is undefined.
751 17 Feb 07 jari 120     */
3736 22 May 18 peter 121     const Vector& s(void) const;
7 28 Feb 03 daniel 122
751 17 Feb 07 jari 123     /**
751 17 Feb 07 jari 124        \brief Solve the system \f$ Ax=b \f$ using the decomposition of
751 17 Feb 07 jari 125        A.
7 28 Feb 03 daniel 126
751 17 Feb 07 jari 127        \note If decompose() has not been run the outcome of the call
751 17 Feb 07 jari 128        is undefined.
751 17 Feb 07 jari 129
751 17 Feb 07 jari 130        \throw GSL_error if the underlying GSL function fails.
751 17 Feb 07 jari 131     */
3736 22 May 18 peter 132     void solve(const VectorBase& b, Vector& x);
751 17 Feb 07 jari 133
751 17 Feb 07 jari 134     /**
751 17 Feb 07 jari 135        \brief Access to the U matrix.
751 17 Feb 07 jari 136
4367 12 Sep 23 peter 137        The column vectors of U are orthonormal, i.e., \f$ U' * U = I
4367 12 Sep 23 peter 138        \f$. The row vectors are typically not orthonormal, i.e., \f$ U
4367 12 Sep 23 peter 139        * U' \f$ does typically not equal the identity matrix, I.
751 17 Feb 07 jari 140
4367 12 Sep 23 peter 141        \return A reference to the U matrix.
4367 12 Sep 23 peter 142
751 17 Feb 07 jari 143        \note If decompose() has not been run the outcome of the call
751 17 Feb 07 jari 144        is undefined.
751 17 Feb 07 jari 145     */
3736 22 May 18 peter 146     const Matrix& U(void) const;
7 28 Feb 03 daniel 147
751 17 Feb 07 jari 148     /**
751 17 Feb 07 jari 149        \brief Access to the V matrix.
751 17 Feb 07 jari 150
4367 12 Sep 23 peter 151        The V matrix is square and orthonormal, i.e., \f$ V' * V = V *
4367 12 Sep 23 peter 152        V' = I \f$
751 17 Feb 07 jari 153
4367 12 Sep 23 peter 154        \return A reference to the V matrix.
4367 12 Sep 23 peter 155
751 17 Feb 07 jari 156        \note If decompose() has not been run the outcome of the call
751 17 Feb 07 jari 157        is undefined.
751 17 Feb 07 jari 158     */
3736 22 May 18 peter 159     const Matrix& V(void) const;
7 28 Feb 03 daniel 160
4 25 Feb 03 daniel 161   private:
751 17 Feb 07 jari 162     /**
751 17 Feb 07 jari 163        \brief Call GSL's Jacobi algorithm for SVD.
751 17 Feb 07 jari 164
751 17 Feb 07 jari 165        \return gsl_error status.
751 17 Feb 07 jari 166     */
715 22 Dec 06 jari 167     int jacobi(void);
751 17 Feb 07 jari 168
751 17 Feb 07 jari 169     /**
751 17 Feb 07 jari 170        \brief Call GSL's Golub-Reinsch algorithm for SVD.
751 17 Feb 07 jari 171
751 17 Feb 07 jari 172        \return gsl_error status.
751 17 Feb 07 jari 173     */
227 01 Feb 05 jari 174     int golub_reinsch(void);
751 17 Feb 07 jari 175
751 17 Feb 07 jari 176     /**
751 17 Feb 07 jari 177        \brief Call GSL's modified Golub-Reinch algorithm for SVD.
751 17 Feb 07 jari 178
751 17 Feb 07 jari 179        \return gsl_error status.
751 17 Feb 07 jari 180     */
227 01 Feb 05 jari 181     int modified_golub_reinsch(void);
7 28 Feb 03 daniel 182
3736 22 May 18 peter 183     Matrix U_, V_;
3736 22 May 18 peter 184     Vector s_;
3736 22 May 18 peter 185   };
4 25 Feb 03 daniel 186
687 16 Oct 06 jari 187 }}} // of namespace utility, yat, and theplu
42 26 Feb 04 jari 188
4 25 Feb 03 daniel 189 #endif

Generated on Sat Oct 14 02:25:15 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)