source: trunk/test/svd_test.cc @ 1151

Last change on this file since 1151 was 1151, checked in by Peter, 15 years ago

remove operator [] in Vector.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.4 KB
Line 
1// $Id: svd_test.cc 1151 2008-02-25 22:32:04Z 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, 2007 Jari Häkkinen
8
9  This file is part of the yat library, http://trac.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 2 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include "yat/random/random.h"
28#include "yat/utility/Matrix.h"
29#include "yat/utility/SVD.h"
30#include "yat/utility/Vector.h"
31
32using namespace theplu::yat;
33
34double this_norm(const utility::Matrix& A)
35{
36  double sum=0.0;
37  for (size_t i=0; i<A.rows(); ++i)
38    for (size_t j=0; j<A.columns(); ++j)
39      sum += A(i,j)*A(i,j);
40  return sum;
41}
42
43
44
45bool test(size_t m, size_t n, utility::SVD::SVDalgorithm algo)
46{
47  double MAXTOL=1e-13; // accepted error
48
49  // initialise a random test-matrix
50  theplu::yat::random::ContinuousUniform rnd;
51  utility::Matrix A(m,n);
52  for (size_t i=0; i<m; ++i)
53    for(size_t j=0; j<n; ++j)
54      A(i,j)=1000*rnd();
55
56  utility::SVD svd(A);
57  svd.decompose(algo);
58  theplu::yat::utility::Vector s(svd.s());
59  utility::Matrix S(s.size(),s.size());
60  for (size_t i=0; i<s.size(); ++i)
61    S(i,i)=s(i);
62  utility::Matrix Vtranspose=svd.V();
63  Vtranspose.transpose();
64  // Reconstructing A = U*S*Vtranspose
65  utility::Matrix Areconstruct(svd.U());
66  Areconstruct*=S;
67  Areconstruct*=Vtranspose;
68  Areconstruct-=A;  // Expect null matrix
69  double error = this_norm(Areconstruct);
70  bool testerror=false;
71  if (error>MAXTOL) {
72    std::cerr << "test_svd: FAILED, algorithm " << algo
73              << " recontruction error ("
74              << error << ") > tolerance (" << MAXTOL << "), matrix dimension ("
75              << m << ',' << n << ')' << std::endl;
76    testerror=true;
77  }
78
79  Vtranspose*=svd.V();  // Expect unity matrix
80  error=this_norm(Vtranspose)-n;
81  if (error>MAXTOL) {
82    std::cerr << "test_svd: FAILED, algorithm " << algo
83              << " V orthogonality error ("
84              << error << ") > tolerance (" << MAXTOL << ')' << std::endl;
85    testerror=true;
86  }
87
88  utility::Matrix Utranspose(svd.U());
89  Utranspose.transpose();
90  Utranspose*=svd.U();  // Expect unity matrix
91  error=this_norm(Utranspose)-n;
92  if (error>MAXTOL) {
93    std::cerr << "test_svd: FAILED, algorithm " << algo
94              << " U orthogonality error ("
95              << error << ") > tolerance (" << MAXTOL << ')' << std::endl;
96    testerror=true;
97  }
98  return testerror;
99}
100
101
102
103int main(const int argc,const char* argv[])
104{
105  bool testfail=false;
106
107  // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
108  // implementations supports rows>=columns matrix dimensions only
109  testfail|=test(12,12,utility::SVD::GolubReinsch);
110  testfail|=test(12,4,utility::SVD::GolubReinsch);
111  testfail|=test(12,12,utility::SVD::ModifiedGolubReinsch);
112  testfail|=test(12,4,utility::SVD::ModifiedGolubReinsch);
113  testfail|=test(12,12,utility::SVD::Jacobi);
114  testfail|=test(12,4,utility::SVD::Jacobi);
115
116  return (testfail ? -1 : 0);
117}
Note: See TracBrowser for help on using the repository browser.