source: trunk/test/test_svd.cc @ 227

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

Started reimplementation of the vector class.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.2 KB
Line 
1// $Id: test_svd.cc 227 2005-02-01 00:52:27Z jari $
2//
3// Testing for a succesful reconstruction of a decomposed random
4// matrix, aswell that U and V are orthogonal.
5
6
7#include "SVD.h"
8#include "matrix.h"
9#include "random_singleton.h"
10#include "vector.h"
11
12using namespace std;
13using namespace theplu::cpptools;
14using namespace theplu::gslapi;
15
16
17
18double this_norm(const matrix& A)
19{
20  double sum=0.0;
21  for (size_t i=0; i<A.rows(); ++i)
22    for (size_t j=0; j<A.columns(); ++j)
23      sum += A(i,j)*A(i,j);
24  return sum;
25}
26
27
28
29bool test(size_t m, size_t n, SVD::SVDalgorithm algo) {
30
31  // accepted error, Jari: should be picked up from GSL
32  double MAXTOL=1e-13;
33  random_singleton* rng=random_singleton::get_instance(-1);
34
35  // initialise a random test-matrix
36  matrix A(m,n);
37  for (size_t i=0; i<m; ++i)
38    for(size_t j=0; j<n; ++j)
39      A(i,j)=1000*rng->get_uniform_double();
40
41  SVD svd(A);
42  svd.decompose(algo);
43  theplu::gslapi::vector s(svd.s());
44  matrix S(s.size(),s.size());
45  for (size_t i=0; i<s.size(); ++i)
46    S(i,i)=s[i];
47  matrix V=svd.V();
48  matrix U=svd.U();
49  double error = this_norm(A-U*S*V.transpose());
50  bool testerror=false;
51  if (error>MAXTOL) {
52    cerr << "test_svd: FAILED, algorithm " << algo << " recontruction error ("
53         << error << ") > tolerance (" << MAXTOL << "), matrix dimension ("
54         << m << ',' << n << ')' << endl;
55    testerror=true;
56  }
57
58  error=this_norm(V.transpose()*V)-n;
59  if (error>MAXTOL) {
60    cerr << "test_svd: FAILED, algorithm " << algo << " V orthogonality error ("
61         << error << ") > tolerance (" << MAXTOL << ')' << endl;
62    testerror=true;
63  }
64
65  error=this_norm(U.transpose()*U)-n;
66  if (error>MAXTOL) {
67    cerr << "test_svd: FAILED, algorithm " << algo << " U orthogonality error ("
68         << error << ") > tolerance (" << MAXTOL << ')' << endl;
69    testerror=true;
70  }
71  return testerror;
72}
73
74
75
76int main(const int argc,const char* argv[])
77{
78  bool testfail=false;
79
80  // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
81  // implementations supports rows>=columns matrix dimensions only
82  testfail|=test(12,12,SVD::GolubReinsch);
83  testfail|=test(12,4,SVD::GolubReinsch);
84  testfail|=test(12,12,SVD::ModifiedGolubReinsch);
85  testfail|=test(12,4,SVD::ModifiedGolubReinsch);
86  testfail|=test(12,12,SVD::Jacobi);
87  testfail|=test(12,4,SVD::Jacobi);
88
89  return (testfail ? -1 : 0);
90}
Note: See TracBrowser for help on using the repository browser.