1 | // $Id: svd_test.cc 759 2007-02-19 19:41:25Z 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 | |
---|
29 | using namespace theplu::yat; |
---|
30 | |
---|
31 | double 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 | |
---|
42 | bool 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 | |
---|
100 | int 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 | } |
---|