test/svd.cc

Code
Comments
Other
Rev Date Author Line
4 25 Feb 03 daniel 1 // $Id$
4 25 Feb 03 daniel 2
675 10 Oct 06 jari 3 /*
831 27 Mar 07 peter 4   Copyright (C) 2003 Daniel Dalevi
2119 12 Dec 09 peter 5   Copyright (C) 2004 Jari Häkkinen
2119 12 Dec 09 peter 6   Copyright (C) 2005 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 7   Copyright (C) 2006 Jari Häkkinen
2119 12 Dec 09 peter 8   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
4207 26 Aug 22 peter 9   Copyright (C) 2012, 2022 Peter Johansson
4 25 Feb 03 daniel 10
1437 25 Aug 08 peter 11   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 12
675 10 Oct 06 jari 13   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 14   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 15   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 16   License, or (at your option) any later version.
675 10 Oct 06 jari 17
675 10 Oct 06 jari 18   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 19   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 21   General Public License for more details.
675 10 Oct 06 jari 22
675 10 Oct 06 jari 23   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 24   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 25 */
675 10 Oct 06 jari 26
2881 18 Nov 12 peter 27 #include <config.h>
2881 18 Nov 12 peter 28
1248 19 Mar 08 peter 29 #include "Suite.h"
1248 19 Mar 08 peter 30
675 10 Oct 06 jari 31 #include "yat/random/random.h"
1121 22 Feb 08 peter 32 #include "yat/utility/Matrix.h"
675 10 Oct 06 jari 33 #include "yat/utility/SVD.h"
1120 21 Feb 08 peter 34 #include "yat/utility/Vector.h"
675 10 Oct 06 jari 35
680 11 Oct 06 jari 36 using namespace theplu::yat;
4 25 Feb 03 daniel 37
1248 19 Mar 08 peter 38 void do_test(size_t m, size_t n, utility::SVD::SVDalgorithm algo,
1248 19 Mar 08 peter 39              test::Suite& suite)
4 25 Feb 03 daniel 40 {
1248 19 Mar 08 peter 41   // initialise a random test-matrix
1248 19 Mar 08 peter 42   theplu::yat::random::ContinuousUniform rnd;
227 01 Feb 05 jari 43
1248 19 Mar 08 peter 44   // set seed to make test deterministic
1248 19 Mar 08 peter 45   theplu::yat::random::RNG::instance()->seed(20);
227 01 Feb 05 jari 46
1121 22 Feb 08 peter 47   utility::Matrix A(m,n);
227 01 Feb 05 jari 48   for (size_t i=0; i<m; ++i)
227 01 Feb 05 jari 49     for(size_t j=0; j<n; ++j)
377 07 Aug 05 jari 50       A(i,j)=1000*rnd();
227 01 Feb 05 jari 51
301 30 Apr 05 peter 52   utility::SVD svd(A);
227 01 Feb 05 jari 53    svd.decompose(algo);
1120 21 Feb 08 peter 54   theplu::yat::utility::Vector s(svd.s());
1121 22 Feb 08 peter 55   utility::Matrix S(s.size(),s.size());
1476 04 Sep 08 peter 56   for (size_t i=1; i<s.size(); ++i)
1476 04 Sep 08 peter 57     if (!suite.add(s(i-1)>=s(i)))
1476 04 Sep 08 peter 58       suite.err() << "single values not in non-increasing order\n";
227 01 Feb 05 jari 59   for (size_t i=0; i<s.size(); ++i)
1151 25 Feb 08 peter 60     S(i,i)=s(i);
1121 22 Feb 08 peter 61   utility::Matrix Vtranspose=svd.V();
420 02 Dec 05 jari 62   Vtranspose.transpose();
420 02 Dec 05 jari 63   // Reconstructing A = U*S*Vtranspose
1121 22 Feb 08 peter 64   utility::Matrix Areconstruct(svd.U());
420 02 Dec 05 jari 65   Areconstruct*=S;
420 02 Dec 05 jari 66   Areconstruct*=Vtranspose;
1251 03 Apr 08 peter 67   if (!suite.equal_range(A.begin(), A.end(), Areconstruct.begin(), 10000)) {
1248 19 Mar 08 peter 68      suite.err() << "test_svd: FAILED, algorithm " << algo << std::endl;
1248 19 Mar 08 peter 69      suite.err() << "reconstruction error" << std::endl;
1248 19 Mar 08 peter 70     suite.add(false);
227 01 Feb 05 jari 71   }
227 01 Feb 05 jari 72
420 02 Dec 05 jari 73   Vtranspose*=svd.V();  // Expect unity matrix
1248 19 Mar 08 peter 74   Vtranspose += 1.0;
1248 19 Mar 08 peter 75   utility::Matrix eye(Vtranspose.rows(), Vtranspose.columns(), 1.0);
1248 19 Mar 08 peter 76   for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
1248 19 Mar 08 peter 77     eye(i,i)=2.0;
1248 19 Mar 08 peter 78   if (!suite.equal_range(eye.begin(), eye.end(), Vtranspose.begin(), 100) ) {
4200 19 Aug 22 peter 79     suite.err() << "test_svd: FAILED, algorithm " << algo
1248 19 Mar 08 peter 80                 << " V orthogonality error " << std::endl;
1248 19 Mar 08 peter 81     suite.add(false);
227 01 Feb 05 jari 82   }
227 01 Feb 05 jari 83
1121 22 Feb 08 peter 84   utility::Matrix Utranspose(svd.U());
420 02 Dec 05 jari 85   Utranspose.transpose();
420 02 Dec 05 jari 86   Utranspose*=svd.U();  // Expect unity matrix
1248 19 Mar 08 peter 87   Utranspose += 1.0;
1248 19 Mar 08 peter 88   eye = utility::Matrix(Utranspose.rows(), Utranspose.columns(), 1.0);
1248 19 Mar 08 peter 89   for (size_t i=0; i<eye.rows() && i<eye.columns(); ++i)
1248 19 Mar 08 peter 90     eye(i,i)=2.0;
1248 19 Mar 08 peter 91   if (!suite.equal_range(eye.begin(), eye.end(), Utranspose.begin(), 100) ) {
4200 19 Aug 22 peter 92     suite.err() << "test_svd: FAILED, algorithm " << algo
1248 19 Mar 08 peter 93                 << " U orthogonality error " << std::endl;
1248 19 Mar 08 peter 94     suite.add(false);
227 01 Feb 05 jari 95   }
227 01 Feb 05 jari 96 }
227 01 Feb 05 jari 97
227 01 Feb 05 jari 98
227 01 Feb 05 jari 99
1248 19 Mar 08 peter 100 int main(int argc, char* argv[])
227 01 Feb 05 jari 101 {
1248 19 Mar 08 peter 102   theplu::yat::test::Suite suite(argc, argv);
227 01 Feb 05 jari 103
227 01 Feb 05 jari 104   // The GSL Jacobi, Golub-Reinsch, and modified Golub-Reinsch
227 01 Feb 05 jari 105   // implementations supports rows>=columns matrix dimensions only
1248 19 Mar 08 peter 106   do_test(12,12,utility::SVD::GolubReinsch, suite);
1248 19 Mar 08 peter 107   do_test(12,4,utility::SVD::GolubReinsch, suite);
1248 19 Mar 08 peter 108   do_test(12,12,utility::SVD::ModifiedGolubReinsch, suite);
1248 19 Mar 08 peter 109   do_test(12,4,utility::SVD::ModifiedGolubReinsch, suite);
1248 19 Mar 08 peter 110   do_test(12,12,utility::SVD::Jacobi, suite);
1248 19 Mar 08 peter 111   do_test(12,4,utility::SVD::Jacobi, suite);
227 01 Feb 05 jari 112
1248 19 Mar 08 peter 113   return suite.return_value();
227 01 Feb 05 jari 114 }

Generated on Tue Apr 25 01:20:39 2023 (UTC) by svndigest 0.10.4pre (dev build 1607)