source: trunk/yat/utility/SVD.cc

Last change on this file was 3999, checked in by Peter, 10 months ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.0 KB
Line 
1// $Id: SVD.cc 3999 2020-10-08 23:22:32Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi
5  Copyright (C) 2004 Jari Häkkinen
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2006 Jari Häkkinen
8  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2009 Jari Häkkinen
10  Copyright (C) 2011, 2012, 2018, 2020 Peter Johansson
11
12  This file is part of the yat library, http://dev.thep.lu.se/yat
13
14  The yat library is free software; you can redistribute it and/or
15  modify it under the terms of the GNU General Public License as
16  published by the Free Software Foundation; either version 3 of the
17  License, or (at your option) any later version.
18
19  The yat library is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22  General Public License for more details.
23
24  You should have received a copy of the GNU General Public License
25  along with yat. If not, see <http://www.gnu.org/licenses/>.
26*/
27
28#include <config.h>
29
30#include "SVD.h"
31
32#include "Exception.h"
33#include "Vector.h"
34#include "VectorBase.h"
35
36#include <cassert>
37#include <sstream>
38#include <utility>
39
40namespace theplu {
41namespace yat {
42namespace utility {
43
44
45  SVD::SVD(const Matrix& Ain)
46    : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
47  {
48  }
49
50
51  SVD::SVD(Matrix&& Ain)
52    : U_(std::move(Ain)), V_(U_.columns(), U_.columns()), s_(U_.columns())
53  {
54  }
55
56
57  SVD::~SVD(void)
58  {
59  }
60
61
62  void SVD::decompose(SVDalgorithm algo)
63  {
64    assert(U_.rows()>=U_.columns());
65    int status=0;
66    switch (algo) {
67    case GolubReinsch:
68      status=golub_reinsch();
69      break;
70    case ModifiedGolubReinsch:
71      status=modified_golub_reinsch();
72      break;
73    case Jacobi:
74      status=jacobi();
75      break;
76    default:
77      std::stringstream ss("SVD::decompose ");
78      ss << algo << " is not a valid SVDalgorithm";
79      throw GSL_error(ss.str());
80    }
81    if (status)
82      throw GSL_error("SVD::decompose",status);
83  }
84
85
86  int SVD::golub_reinsch(void)
87  {
88    Vector w(U_.columns());
89    return gsl_linalg_SV_decomp(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
90                                s_.gsl_vector_p(), w.gsl_vector_p());
91  }
92
93
94  int SVD::jacobi(void)
95  {
96    return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
97                                       s_.gsl_vector_p());
98  }
99
100
101  int SVD::modified_golub_reinsch(void)
102  {
103    Vector w(U_.columns());
104    Matrix X(U_.columns(),U_.columns());
105    return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_p(), X.gsl_matrix_p(), 
106                                    V_.gsl_matrix_p(), s_.gsl_vector_p(),
107                                    w.gsl_vector_p());
108  }
109
110
111  const Vector& SVD::s(void) const
112  {
113    return s_;
114  }
115
116
117  void SVD::solve(const VectorBase& b, Vector& x)
118  {
119    int status=gsl_linalg_SV_solve(U_.gsl_matrix_p(), V_.gsl_matrix_p(), 
120                                   s_.gsl_vector_p(), b.gsl_vector_p(),
121                                   x.gsl_vector_p());
122    if (status)
123      throw GSL_error("SVD::solve",status);
124  }
125
126
127  const Matrix& SVD::U(void) const
128  {
129    return U_;
130  }
131
132
133  const Matrix& SVD::V(void) const
134  {
135    return V_;
136  }
137
138}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.