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

Last change on this file since 1437 was 1437, checked in by Peter, 13 years ago

merge patch release 0.4.2 to trunk. Delta 0.4.2-0.4.1

  • 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 1437 2008-08-25 17:55:00Z peter $
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
15  published by the Free Software Foundation; either version 2 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 this program; if not, write to the Free Software
25  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
26  02111-1307, USA.
27*/
28
[680]29#include "SVD.h"
[1120]30#include "Vector.h"
[1016]31#include "VectorBase.h"
32
[751]33#include <sstream>
[675]34
[42]35namespace theplu {
[680]36namespace yat {
[301]37namespace utility {
[42]38
[4]39
[1121]40  SVD::SVD(const utility::Matrix& Ain)
[703]41    : U_(Ain), V_(Ain.columns(),Ain.columns()), s_(Ain.columns())
42  {
43  }
44
45
46  SVD::~SVD(void)
47  {
48  }
49
50
[751]51  void SVD::decompose(SVDalgorithm algo)
[227]52  {
[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.