source: trunk/test/svd_test.cc @ 792

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

Addresses #193. matrix now works as outlined in the ticket
discussion. Added support for const views. Added a clone function that
facilitates resizing of matrices. clone is needed since assignement
operator functionality is changed.

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