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

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

Added missing include

  • 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 1723 2009-01-15 13:36:32Z 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  Copyright (C) 2009 Jari Häkkinen
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 "SVD.h"
29
30#include "yat/utility/Exception.h"
31#include "Vector.h"
32#include "VectorBase.h"
33
34#include <cassert>
35#include <sstream>
36
37namespace theplu {
38namespace yat {
39namespace utility {
40
41
42  SVD::SVD(const utility::Matrix& Ain)
43    : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
44  {
45  }
46
47
48  SVD::~SVD(void)
49  {
50  }
51
52
53  void SVD::decompose(SVDalgorithm algo)
54  {
55    assert(U_.rows()>=U_.columns());
56    int status=0;
57    switch (algo) {
58    case GolubReinsch:
59      status=golub_reinsch();
60      break;
61    case ModifiedGolubReinsch:
62      status=modified_golub_reinsch();
63      break;
64    case Jacobi:
65      status=jacobi();
66      break;
67    default:
68      std::stringstream ss("SVD::decompose ");
69      ss << algo << " is not a valid SVDalgorithm";
70      throw GSL_error(ss.str());
71    }
72    if (status)
73      throw utility::GSL_error(std::string("SVD::decompose",status));
74  }
75
76
77  int SVD::golub_reinsch(void)
78  {
79    utility::Vector w(U_.columns());
80    return gsl_linalg_SV_decomp(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
81                                s_.gsl_vector_p(), w.gsl_vector_p());
82  }
83
84
85  int SVD::jacobi(void)
86  {
87    return gsl_linalg_SV_decomp_jacobi(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
88                                       s_.gsl_vector_p());
89  }
90
91
92  int SVD::modified_golub_reinsch(void)
93  {
94    utility::Vector w(U_.columns());
95    utility::Matrix X(U_.columns(),U_.columns());
96    return gsl_linalg_SV_decomp_mod(U_.gsl_matrix_p(), X.gsl_matrix_p(), 
97                                    V_.gsl_matrix_p(), s_.gsl_vector_p(),
98                                    w.gsl_vector_p());
99  }
100
101
102  const utility::Vector& SVD::s(void) const
103  {
104    return s_;
105  }
106
107
108  void SVD::solve(const utility::VectorBase& b, utility::Vector& x)
109  {
110    int status=gsl_linalg_SV_solve(U_.gsl_matrix_p(), V_.gsl_matrix_p(), 
111                                   s_.gsl_vector_p(), b.gsl_vector_p(),
112                                   x.gsl_vector_p());
113    if (status)
114      throw utility::GSL_error(std::string("SVD::solve",status));
115  }
116
117
118  const utility::Matrix& SVD::U(void) const
119  {
120    return U_;
121  }
122
123
124  const utility::Matrix& SVD::V(void) const
125  {
126    return V_;
127  }
128
129}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.