source: trunk/test/svd_test.cc @ 1619

Last change on this file since 1619 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.6 KB
Line 
1// $Id: svd_test.cc 1487 2008-09-10 08:41:36Z jari $
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 Jari Häkkinen
8  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2008 Peter Johansson
10
11  This file is part of the yat library, http://dev.thep.lu.se/yat
12
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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
25*/
26
27#include "Suite.h"
28
29#include "yat/random/random.h"
30#include "yat/utility/Matrix.h"
31#include "yat/utility/SVD.h"
32#include "yat/utility/Vector.h"
33
34using namespace theplu::yat;
35
36void do_test(size_t m, size_t n, utility::SVD::SVDalgorithm algo,
37             test::Suite& suite)
38{
39  // initialise a random test-matrix
40  theplu::yat::random::ContinuousUniform rnd;
41
42  // set seed to make test deterministic
43  theplu::yat::random::RNG::instance()->seed(20);
44
45  utility::Matrix A(m,n);
46  for (size_t i=0; i<m; ++i)
47    for(size_t j=0; j<n; ++j)
48      A(i,j)=1000*rnd();
49
50  utility::SVD svd(A);
51  svd.decompose(algo);
52  theplu::yat::utility::Vector s(svd.s());
53  utility::Matrix S(s.size(),s.size());
54  for (size_t i=1; i<s.size(); ++i)
55    if (!suite.add(s(i-1)>=s(i)))
56      suite.err() << "single values not in non-increasing order\n";
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  if (!suite.equal_range(A.begin(), A.end(), Areconstruct.begin(), 10000)) {
66    suite.err() << "test_svd: FAILED, algorithm " << algo << std::endl;
67    suite.err() << "reconstruction error" << std::endl;
68    suite.add(false);
69  }
70
71  Vtranspose*=svd.V();  // Expect unity matrix
72  Vtranspose += 1.0;
73  utility::Matrix eye(Vtranspose.rows(), Vtranspose.columns(), 1.0);
74  for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
75    eye(i,i)=2.0;
76  if (!suite.equal_range(eye.begin(), eye.end(), Vtranspose.begin(), 100) ) {
77    suite.err() << "test_svd: FAILED, algorithm " << algo
78                << " V orthogonality error " << std::endl;
79    suite.add(false);
80  }
81
82  utility::Matrix Utranspose(svd.U());
83  Utranspose.transpose();
84  Utranspose*=svd.U();  // Expect unity matrix
85  Utranspose += 1.0;
86  eye = utility::Matrix(Utranspose.rows(), Utranspose.columns(), 1.0);
87  for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
88    eye(i,i)=2.0;
89  if (!suite.equal_range(eye.begin(), eye.end(), Utranspose.begin(), 100) ) {
90    suite.err() << "test_svd: FAILED, algorithm " << algo
91                << " U orthogonality error " << std::endl;
92    suite.add(false);
93  }
94}
95
96
97
98int main(int argc, char* argv[])
99{
100  theplu::yat::test::Suite suite(argc, argv);
101
102  // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
103  // implementations supports rows>=columns matrix dimensions only
104  do_test(12,12,utility::SVD::GolubReinsch, suite);
105  do_test(12,4,utility::SVD::GolubReinsch, suite);
106  do_test(12,12,utility::SVD::ModifiedGolubReinsch, suite);
107  do_test(12,4,utility::SVD::ModifiedGolubReinsch, suite);
108  do_test(12,12,utility::SVD::Jacobi, suite);
109  do_test(12,4,utility::SVD::Jacobi, suite);
110
111  return suite.return_value();
112}
Note: See TracBrowser for help on using the repository browser.