source: trunk/yat/utility/SVD.cc @ 1706

Last change on this file since 1706 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

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