source: trunk/test/svd_test.cc @ 301

Last change on this file since 301 was 301, checked in by Peter, 17 years ago

modified includes in tests

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