test/cholesky_decomposer.cc

Code
Comments
Other
Rev Date Author Line
4280 27 Jan 23 peter 1 // $Id$
4280 27 Jan 23 peter 2
4280 27 Jan 23 peter 3 /*
4280 27 Jan 23 peter 4   Copyright (C) 2023 Peter Johansson
4280 27 Jan 23 peter 5
4280 27 Jan 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4280 27 Jan 23 peter 7
4280 27 Jan 23 peter 8   The yat library is free software; you can redistribute it and/or
4280 27 Jan 23 peter 9   modify it under the terms of the GNU General Public License as
4280 27 Jan 23 peter 10   published by the Free Software Foundation; either version 3 of the
4280 27 Jan 23 peter 11   License, or (at your option) any later version.
4280 27 Jan 23 peter 12
4280 27 Jan 23 peter 13   The yat library is distributed in the hope that it will be useful,
4280 27 Jan 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4280 27 Jan 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4280 27 Jan 23 peter 16   General Public License for more details.
4280 27 Jan 23 peter 17
4280 27 Jan 23 peter 18   You should have received a copy of the GNU General Public License
4280 27 Jan 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4280 27 Jan 23 peter 20 */
4280 27 Jan 23 peter 21
4280 27 Jan 23 peter 22 #include <config.h>
4280 27 Jan 23 peter 23
4280 27 Jan 23 peter 24 #include "Suite.h"
4280 27 Jan 23 peter 25
4280 27 Jan 23 peter 26 #include "yat/utility/CholeskyDecomposer.h"
4280 27 Jan 23 peter 27 #include "yat/utility/Matrix.h"
4280 27 Jan 23 peter 28 #include "yat/utility/Vector.h"
4280 27 Jan 23 peter 29
4280 27 Jan 23 peter 30 using namespace theplu::yat;
4280 27 Jan 23 peter 31
4280 27 Jan 23 peter 32 int main(int argc, char* argv[])
4280 27 Jan 23 peter 33 {
4280 27 Jan 23 peter 34   test::Suite suite(argc, argv);
4280 27 Jan 23 peter 35
4280 27 Jan 23 peter 36   utility::Matrix L = test::generate_Matrix(10, 10, 0.1);
4280 27 Jan 23 peter 37   for (size_t i=0; i<L.rows(); ++i)
4280 27 Jan 23 peter 38     for (size_t j=i+1; j<L.columns(); ++j)
4280 27 Jan 23 peter 39       L(i,j) = 0.0;
4280 27 Jan 23 peter 40   utility::Matrix A = L * transpose(L);
4280 27 Jan 23 peter 41   utility::CholeskyDecomposer decomposer(A);
4280 27 Jan 23 peter 42   for (size_t i=0; i<decomposer.L().rows(); ++i)
4280 27 Jan 23 peter 43     for (size_t j=i+1; j<decomposer.L().columns(); ++j)
4280 27 Jan 23 peter 44       if (decomposer.L()(i, j))
4280 27 Jan 23 peter 45         suite.err() << "L(" << i << ", " << j << ") non-zero: "
4280 27 Jan 23 peter 46                     << decomposer.L()(i, j) << "\n";
4280 27 Jan 23 peter 47
4280 27 Jan 23 peter 48   if (decomposer.L().rows() != 10 || decomposer.L().columns()!=10) {
4280 27 Jan 23 peter 49     suite.add(false);
4280 27 Jan 23 peter 50     suite.err() << "error: L incorrect dimensions: "
4280 27 Jan 23 peter 51                 << decomposer.L().rows() << " x " << decomposer.L().columns()
4280 27 Jan 23 peter 52                 << "\n";
4280 27 Jan 23 peter 53   }
4280 27 Jan 23 peter 54   else {
4280 27 Jan 23 peter 55     utility::Matrix Y = decomposer.L() * transpose(decomposer.L());
4341 16 Apr 23 peter 56     if (!suite.equal_matrix(Y, A, 2)) {
4280 27 Jan 23 peter 57       suite.add(false);
4280 27 Jan 23 peter 58       suite.err() << "error: X != L * transpose(L)\n";
4280 27 Jan 23 peter 59     }
4280 27 Jan 23 peter 60   }
4280 27 Jan 23 peter 61
4280 27 Jan 23 peter 62   utility::Matrix B;
4280 27 Jan 23 peter 63   decomposer.invert(B);
4280 27 Jan 23 peter 64   if (B.rows() != A.rows() || B.columns() != A.columns()) {
4280 27 Jan 23 peter 65     suite.add(false);
4280 27 Jan 23 peter 66     suite.err() << "inverse: wrong dimensions: " << B.rows()
4280 27 Jan 23 peter 67                 << " x " << B.columns() << "\n";
4280 27 Jan 23 peter 68   }
4280 27 Jan 23 peter 69   else {
4280 27 Jan 23 peter 70     utility::Matrix Eye = B * A;
4280 27 Jan 23 peter 71     for (size_t i=0; i<Eye.rows(); ++i) {
4280 27 Jan 23 peter 72       for (size_t j=0; j<Eye.columns(); ++j) {
4280 27 Jan 23 peter 73         if (!suite.add(suite.equal_fix(Eye(i,j), i==j ? 1.0 : 0.0, 1e-9)))
4280 27 Jan 23 peter 74           suite.err() << "error: row: " << i << " column: " << j << "\n";
4280 27 Jan 23 peter 75       }
4280 27 Jan 23 peter 76     }
4280 27 Jan 23 peter 77   }
4280 27 Jan 23 peter 78
4280 27 Jan 23 peter 79   utility::Vector b(10, 1.0);
4280 27 Jan 23 peter 80   b(0) = 0.5;
4280 27 Jan 23 peter 81   b(1) = 2.13;
4280 27 Jan 23 peter 82   b(2) = 3.14;
4280 27 Jan 23 peter 83   utility::Vector x(10);
4280 27 Jan 23 peter 84   decomposer.solve(b, x);
4280 27 Jan 23 peter 85   if (x.size() != 10) {
4280 27 Jan 23 peter 86     suite.add(false);
4280 27 Jan 23 peter 87     suite.err() << "::solve: incorrect x size: " << x.size() << "\n";
4280 27 Jan 23 peter 88   }
4280 27 Jan 23 peter 89   else {
4280 27 Jan 23 peter 90     utility::Vector Ax = A * x;
4378 11 Oct 23 peter 91     if (!suite.equal_range_fix(Ax.begin(), Ax.end(), b.begin(), 1.24e-10)) {
4280 27 Jan 23 peter 92       suite.add(false);
4280 27 Jan 23 peter 93       suite.err() << "b: " << b << "\n";
4280 27 Jan 23 peter 94       suite.err() << "A*x: " << Ax << "\n";
4280 27 Jan 23 peter 95     }
4280 27 Jan 23 peter 96   }
4280 27 Jan 23 peter 97   return suite.return_value();
4280 27 Jan 23 peter 98 }

Generated on Sat Oct 14 02:20:44 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)