yat/utility/SVD.cc

Code
Comments
Other
Rev Date Author Line
4 25 Feb 03 daniel 1 // $Id$
4 25 Feb 03 daniel 2
675 10 Oct 06 jari 3 /*
831 27 Mar 07 peter 4   Copyright (C) 2003 Daniel Dalevi
2119 12 Dec 09 peter 5   Copyright (C) 2004 Jari Häkkinen
2119 12 Dec 09 peter 6   Copyright (C) 2005 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 7   Copyright (C) 2006 Jari Häkkinen
2119 12 Dec 09 peter 8   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 9   Copyright (C) 2009 Jari Häkkinen
4359 23 Aug 23 peter 10   Copyright (C) 2012, 2018, 2020, 2022, 2023 Peter Johansson
4 25 Feb 03 daniel 11
1437 25 Aug 08 peter 12   This file is part of the yat library, http://dev.thep.lu.se/yat
42 26 Feb 04 jari 13
675 10 Oct 06 jari 14   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 15   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 16   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 17   License, or (at your option) any later version.
675 10 Oct 06 jari 18
675 10 Oct 06 jari 19   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 20   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 21   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 22   General Public License for more details.
675 10 Oct 06 jari 23
675 10 Oct 06 jari 24   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 25   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 26 */
675 10 Oct 06 jari 27
2881 18 Nov 12 peter 28 #include <config.h>
2881 18 Nov 12 peter 29
680 11 Oct 06 jari 30 #include "SVD.h"
1723 15 Jan 09 jari 31
4286 30 Jan 23 peter 32 #include "DiagonalMatrix.h"
2414 15 Jan 11 peter 33 #include "Exception.h"
1120 21 Feb 08 peter 34 #include "Vector.h"
1016 01 Feb 08 peter 35 #include "VectorBase.h"
1016 01 Feb 08 peter 36
1478 06 Sep 08 peter 37 #include <cassert>
751 17 Feb 07 jari 38 #include <sstream>
3736 22 May 18 peter 39 #include <utility>
675 10 Oct 06 jari 40
42 26 Feb 04 jari 41 namespace theplu {
680 11 Oct 06 jari 42 namespace yat {
301 30 Apr 05 peter 43 namespace utility {
42 26 Feb 04 jari 44
4 25 Feb 03 daniel 45
4125 14 Jan 22 peter 46   SVD::SVD(const MatrixBase& Ain)
703 18 Dec 06 jari 47     : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
703 18 Dec 06 jari 48   {
703 18 Dec 06 jari 49   }
703 18 Dec 06 jari 50
703 18 Dec 06 jari 51
4125 14 Jan 22 peter 52   SVD::SVD(MatrixMutable&& Ain)
3977 23 Aug 20 peter 53     : U_(std::move(Ain)), V_(U_.columns(), U_.columns()), s_(U_.columns())
3736 22 May 18 peter 54   {
3736 22 May 18 peter 55   }
3736 22 May 18 peter 56
3736 22 May 18 peter 57
703 18 Dec 06 jari 58   SVD::~SVD(void)
703 18 Dec 06 jari 59   {
703 18 Dec 06 jari 60   }
703 18 Dec 06 jari 61
703 18 Dec 06 jari 62
751 17 Feb 07 jari 63   void SVD::decompose(SVDalgorithm algo)
227 01 Feb 05 jari 64   {
1478 06 Sep 08 peter 65     assert(U_.rows()>=U_.columns());
751 17 Feb 07 jari 66     int status=0;
227 01 Feb 05 jari 67     switch (algo) {
751 17 Feb 07 jari 68     case GolubReinsch:
751 17 Feb 07 jari 69       status=golub_reinsch();
751 17 Feb 07 jari 70       break;
751 17 Feb 07 jari 71     case ModifiedGolubReinsch:
751 17 Feb 07 jari 72       status=modified_golub_reinsch();
751 17 Feb 07 jari 73       break;
751 17 Feb 07 jari 74     case Jacobi:
751 17 Feb 07 jari 75       status=jacobi();
751 17 Feb 07 jari 76       break;
751 17 Feb 07 jari 77     default:
751 17 Feb 07 jari 78       std::stringstream ss("SVD::decompose ");
751 17 Feb 07 jari 79       ss << algo << " is not a valid SVDalgorithm";
751 17 Feb 07 jari 80       throw GSL_error(ss.str());
616 31 Aug 06 jari 81     }
751 17 Feb 07 jari 82     if (status)
3743 12 Jul 18 peter 83       throw GSL_error("SVD::decompose",status);
227 01 Feb 05 jari 84   }
4 25 Feb 03 daniel 85
4 25 Feb 03 daniel 86
227 01 Feb 05 jari 87   int SVD::golub_reinsch(void)
227 01 Feb 05 jari 88   {
3736 22 May 18 peter 89     Vector w(U_.columns());
420 02 Dec 05 jari 90     return gsl_linalg_SV_decomp(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
420 02 Dec 05 jari 91                                 s_.gsl_vector_p(), w.gsl_vector_p());
227 01 Feb 05 jari 92   }
4 25 Feb 03 daniel 93
4 25 Feb 03 daniel 94
715 22 Dec 06 jari 95   int SVD::jacobi(void)
715 22 Dec 06 jari 96   {
715 22 Dec 06 jari 97     return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
715 22 Dec 06 jari 98                                        s_.gsl_vector_p());
715 22 Dec 06 jari 99   }
4 25 Feb 03 daniel 100
715 22 Dec 06 jari 101
227 01 Feb 05 jari 102   int SVD::modified_golub_reinsch(void)
227 01 Feb 05 jari 103   {
3736 22 May 18 peter 104     Vector w(U_.columns());
3736 22 May 18 peter 105     Matrix X(U_.columns(),U_.columns());
4200 19 Aug 22 peter 106     return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_p(), X.gsl_matrix_p(),
420 02 Dec 05 jari 107                                     V_.gsl_matrix_p(), s_.gsl_vector_p(),
420 02 Dec 05 jari 108                                     w.gsl_vector_p());
227 01 Feb 05 jari 109   }
4 25 Feb 03 daniel 110
4 25 Feb 03 daniel 111
4286 30 Jan 23 peter 112   void SVD::inverse(Matrix& result) const
4286 30 Jan 23 peter 113   {
4286 30 Jan 23 peter 114     assert(U_.rows()==U_.columns() && "input matrix must be square");
4286 30 Jan 23 peter 115     // A = U S V'
4286 30 Jan 23 peter 116
4286 30 Jan 23 peter 117     // result = A^-1 = V * S^-1 * U'
4286 30 Jan 23 peter 118
4286 30 Jan 23 peter 119     DiagonalMatrix D(s_.size(), s_.size());
4286 30 Jan 23 peter 120     // If eigenvalue is zero, A is not invertable, perhaps that should
4286 30 Jan 23 peter 121     // trigger an exception
4286 30 Jan 23 peter 122     for (size_t i=0; i<D.rows(); ++i)
4286 30 Jan 23 peter 123       D(i) = 1.0 / s_(i);
4286 30 Jan 23 peter 124
4286 30 Jan 23 peter 125     result = V_ * D * transpose(U_);
4286 30 Jan 23 peter 126   }
4286 30 Jan 23 peter 127
4286 30 Jan 23 peter 128
3736 22 May 18 peter 129   const Vector& SVD::s(void) const
715 22 Dec 06 jari 130   {
715 22 Dec 06 jari 131     return s_;
715 22 Dec 06 jari 132   }
715 22 Dec 06 jari 133
715 22 Dec 06 jari 134
3736 22 May 18 peter 135   void SVD::solve(const VectorBase& b, Vector& x)
715 22 Dec 06 jari 136   {
4200 19 Aug 22 peter 137     int status=gsl_linalg_SV_solve(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
751 17 Feb 07 jari 138                                    s_.gsl_vector_p(), b.gsl_vector_p(),
751 17 Feb 07 jari 139                                    x.gsl_vector_p());
751 17 Feb 07 jari 140     if (status)
3743 12 Jul 18 peter 141       throw GSL_error("SVD::solve",status);
715 22 Dec 06 jari 142   }
715 22 Dec 06 jari 143
715 22 Dec 06 jari 144
3736 22 May 18 peter 145   const Matrix& SVD::U(void) const
715 22 Dec 06 jari 146   {
715 22 Dec 06 jari 147     return U_;
715 22 Dec 06 jari 148   }
715 22 Dec 06 jari 149
715 22 Dec 06 jari 150
3736 22 May 18 peter 151   const Matrix& SVD::V(void) const
715 22 Dec 06 jari 152   {
715 22 Dec 06 jari 153     return V_;
715 22 Dec 06 jari 154   }
715 22 Dec 06 jari 155
687 16 Oct 06 jari 156 }}} // of namespace utility, yat, and theplu

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