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

Last change on this file since 1566 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
RevLine 
[4]1// $Id: SVD.cc 1487 2008-09-10 08:41:36Z jari $
2
[675]3/*
[831]4  Copyright (C) 2003 Daniel Dalevi
5  Copyright (C) 2004 Jari Häkkinen
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
[1275]7  Copyright (C) 2006 Jari Häkkinen
8  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2008 Peter Johansson
[4]10
[1437]11  This file is part of the yat library, http://dev.thep.lu.se/yat
[42]12
[675]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
[1486]15  published by the Free Software Foundation; either version 3 of the
[675]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
[1487]24  along with yat. If not, see <http://www.gnu.org/licenses/>.
[675]25*/
26
[680]27#include "SVD.h"
[1120]28#include "Vector.h"
[1016]29#include "VectorBase.h"
30
[1478]31#include <cassert>
[751]32#include <sstream>
[675]33
[42]34namespace theplu {
[680]35namespace yat {
[301]36namespace utility {
[42]37
[4]38
[1121]39  SVD::SVD(const utility::Matrix& Ain)
[703]40    : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
41  {
42  }
43
44
45  SVD::~SVD(void)
46  {
47  }
48
49
[751]50  void SVD::decompose(SVDalgorithm algo)
[227]51  {
[1478]52    assert(U_.rows()>=U_.columns());
[751]53    int status=0;
[227]54    switch (algo) {
[751]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());
[616]68    }
[751]69    if (status)
70      throw utility::GSL_error(std::string("SVD::decompose",status));
[227]71  }
[4]72
73
[227]74  int SVD::golub_reinsch(void)
75  {
[1120]76    utility::Vector w(U_.columns());
[420]77    return gsl_linalg_SV_decomp(U_.gsl_matrix_p(), V_.gsl_matrix_p(),
78                                s_.gsl_vector_p(), w.gsl_vector_p());
[227]79  }
[4]80
81
[715]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  }
[4]87
[715]88
[227]89  int SVD::modified_golub_reinsch(void)
90  {
[1120]91    utility::Vector w(U_.columns());
[1121]92    utility::Matrix X(U_.columns(),U_.columns());
[420]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());
[227]96  }
[4]97
98
[1120]99  const utility::Vector& SVD::s(void) const
[715]100  {
101    return s_;
102  }
103
104
[1120]105  void SVD::solve(const utility::VectorBase& b, utility::Vector& x)
[715]106  {
[751]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));
[715]112  }
113
114
[1121]115  const utility::Matrix& SVD::U(void) const
[715]116  {
117    return U_;
118  }
119
120
[1121]121  const utility::Matrix& SVD::V(void) const
[715]122  {
123    return V_;
124  }
125
[687]126}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.