yat/utility/CholeskyDecomposer.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 "CholeskyDecomposer.h"
4280 27 Jan 23 peter 25 #include "Exception.h"
4280 27 Jan 23 peter 26 #include "Matrix.h"
4280 27 Jan 23 peter 27 #include "VectorMutable.h"
4280 27 Jan 23 peter 28
4280 27 Jan 23 peter 29 #include <gsl/gsl_linalg.h>
4280 27 Jan 23 peter 30
4280 27 Jan 23 peter 31 #include <cassert>
4280 27 Jan 23 peter 32
4280 27 Jan 23 peter 33 namespace theplu {
4280 27 Jan 23 peter 34 namespace yat {
4280 27 Jan 23 peter 35 namespace utility {
4280 27 Jan 23 peter 36
4280 27 Jan 23 peter 37   CholeskyDecomposer::CholeskyDecomposer(const MatrixBase& A)
4280 27 Jan 23 peter 38   {
4280 27 Jan 23 peter 39     assert(A.rows() == A.columns());
4280 27 Jan 23 peter 40     decompose(Matrix(A));
4280 27 Jan 23 peter 41   }
4280 27 Jan 23 peter 42
4280 27 Jan 23 peter 43
4280 27 Jan 23 peter 44   CholeskyDecomposer::CholeskyDecomposer(Matrix&& A)
4280 27 Jan 23 peter 45   {
4280 27 Jan 23 peter 46     assert(A.rows() == A.columns());
4280 27 Jan 23 peter 47     decompose(std::move(A));
4280 27 Jan 23 peter 48   }
4280 27 Jan 23 peter 49
4280 27 Jan 23 peter 50
4280 27 Jan 23 peter 51   void CholeskyDecomposer::decompose(Matrix&& A)
4280 27 Jan 23 peter 52   {
4280 27 Jan 23 peter 53
4280 27 Jan 23 peter 54     L_ = std::move(A);
4280 27 Jan 23 peter 55     if (int status = gsl_linalg_cholesky_decomp1(L_.gsl_matrix_p()))
4280 27 Jan 23 peter 56       throw GSL_error("CholeskyDecomposer failed", status);
4280 27 Jan 23 peter 57     // zero elements above diagonal
4280 27 Jan 23 peter 58     for (size_t row=0; row+1<L_.rows(); ++row)
4280 27 Jan 23 peter 59       std::fill(L_.begin_row(row) + (row+1), L_.end_row(row), 0.0);
4280 27 Jan 23 peter 60     assert(L_.rows() == L_.columns());
4280 27 Jan 23 peter 61   }
4280 27 Jan 23 peter 62
4280 27 Jan 23 peter 63
4280 27 Jan 23 peter 64   void CholeskyDecomposer::invert(Matrix& X) const
4280 27 Jan 23 peter 65   {
4280 27 Jan 23 peter 66     X = L_;
4280 27 Jan 23 peter 67     if (int status = gsl_linalg_cholesky_invert(X.gsl_matrix_p()))
4280 27 Jan 23 peter 68       throw GSL_error("CholeskyDecomposer::invert failed", status);
4280 27 Jan 23 peter 69   }
4280 27 Jan 23 peter 70
4280 27 Jan 23 peter 71
4280 27 Jan 23 peter 72   const Matrix& CholeskyDecomposer::L(void) const
4280 27 Jan 23 peter 73   {
4280 27 Jan 23 peter 74     return L_;
4280 27 Jan 23 peter 75   }
4280 27 Jan 23 peter 76
4280 27 Jan 23 peter 77
4280 27 Jan 23 peter 78   void CholeskyDecomposer::solve(const VectorBase& b, VectorMutable& x) const
4280 27 Jan 23 peter 79   {
4280 27 Jan 23 peter 80     if (int status=gsl_linalg_cholesky_solve(L_.gsl_matrix_p(),
4280 27 Jan 23 peter 81                                              b.gsl_vector_p(),
4280 27 Jan 23 peter 82                                              x.gsl_vector_p()))
4280 27 Jan 23 peter 83       throw GSL_error("CholeskyDecomposer::solve failed", status);
4280 27 Jan 23 peter 84   }
4280 27 Jan 23 peter 85
4280 27 Jan 23 peter 86 }}} // of namespace utility, yat, and theplu

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