From 7f272457aed89adc2b56e7eac9c17e1cffcc48c6 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 6 Dec 2017 18:23:18 +0100 Subject: [PATCH] lapack: test solve() with Cholesky factorization, improve documentation --- include/deal.II/lac/lapack_full_matrix.h | 15 +++-- tests/lapack/full_matrix_25.cc | 79 ++++++++++++++++++++++++ tests/lapack/full_matrix_25.output | 2 + 3 files changed, 90 insertions(+), 6 deletions(-) create mode 100644 tests/lapack/full_matrix_25.cc create mode 100644 tests/lapack/full_matrix_25.output diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index daf294a1f1..a1be34d185 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -445,8 +445,10 @@ public: number reciprocal_condition_number(const number l1_norm) const; /** - * Same as above but for triangular matrices. The matrix has to have - * the LAPACKSupport::Property set, see set_property(). + * Estimate the reciprocal of the condition number $1/k(A)$ in $L_1$ norm + * for triangular matrices. The matrix has to have + * the LAPACKSupport::Property set to either LAPACKSupport::Property::upper_triangular + * or LAPACKSupport::Property::lower_triangular, see set_property(). */ number reciprocal_condition_number() const; @@ -485,21 +487,22 @@ public: void invert (); /** - * Solve the linear system with right hand side. The matrix should - * be either triangular or LU factorization should be previously computed. + * Solve the linear system with right hand side @p v and put the solution + * back to @p v. The matrix should be either triangular or LU/Cholesky + * factorization should be previously computed. * * The flag transposed indicates whether the solution of the transposed * system is to be performed. */ void solve(Vector &v, - const bool transposed) const; + const bool transposed = false) const; /** * Same as above but for multiple right hand sides (as many as there * are columns in the matrix @p B). */ void solve(LAPACKFullMatrix &B, - const bool transposed) const; + const bool transposed = false) const; /** * Solve the linear system with right hand side given by applying diff --git a/tests/lapack/full_matrix_25.cc b/tests/lapack/full_matrix_25.cc new file mode 100644 index 0000000000..0b5fc92a89 --- /dev/null +++ b/tests/lapack/full_matrix_25.cc @@ -0,0 +1,79 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::solve() for Cholesky factorization + +/* MWE for size=3 in Octave: +R = [10,2,3; 2, 20, 6; 3, 6, 90] +x = [2; -7; 1] + +> R\x +ans = + + 0.268693 + -0.385220 + 0.027836 +*/ + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test() +{ + const unsigned int size = 3; + LAPACKFullMatrix M(size); + M.set_property(LAPACKSupport::symmetric); + + M(0,0) = 10; + M(0,1) = 2; + M(1,0) = 2; + M(0,2) = 3; + M(2,0) = 3; + M(1,1) = 20; + M(1,2) = 6; + M(2,1) = 6; + M(2,2) = 90; + + Vector x(size), y(size); + x[0] = 2; + x[1] = -7; + x[2] = 1; + + y = x; + M.compute_cholesky_factorization(); + M.solve(y); + y.print(deallog.get_file_stream(), 6, false); +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + test(); + +} diff --git a/tests/lapack/full_matrix_25.output b/tests/lapack/full_matrix_25.output new file mode 100644 index 0000000000..4f8485e70d --- /dev/null +++ b/tests/lapack/full_matrix_25.output @@ -0,0 +1,2 @@ + +0.268693 -0.385220 0.027836 -- 2.39.5