From 6b9f04834e425728533ed688ce07ea322d8a7ef7 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 6 Dec 2017 16:00:36 +0100 Subject: [PATCH] lapack: add solve() for triangular matrices or factorized (LU/Cholesky) states --- include/deal.II/lac/lapack_full_matrix.h | 23 ++++++ include/deal.II/lac/lapack_templates.h | 60 ++++++++++++++ source/lac/lapack_full_matrix.cc | 101 +++++++++++++++++++---- tests/lapack/full_matrix_24.cc | 90 ++++++++++++++++++++ tests/lapack/full_matrix_24.output | 3 + 5 files changed, 263 insertions(+), 14 deletions(-) create mode 100644 tests/lapack/full_matrix_24.cc create mode 100644 tests/lapack/full_matrix_24.output diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 83660f69b9..daf294a1f1 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -484,6 +484,23 @@ public: */ void invert (); + /** + * Solve the linear system with right hand side. The matrix should + * be either triangular or LU 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; + + /** + * 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; + /** * Solve the linear system with right hand side given by applying * forward/backward substitution to the previously computed LU @@ -491,7 +508,10 @@ public: * * The flag transposed indicates whether the solution of the transposed * system is to be performed. + * + * @deprecated use solve() instead. */ + DEAL_II_DEPRECATED void apply_lu_factorization (Vector &v, const bool transposed) const; @@ -503,7 +523,10 @@ public: * * The flag transposed indicates whether the solution of the transposed * system is to be performed. + * + * @deprecated use solve() instead. */ + DEAL_II_DEPRECATED void apply_lu_factorization (LAPACKFullMatrix &B, const bool transposed) const; diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index a48ee3739a..8696a5d382 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -76,6 +76,15 @@ extern "C" const int *lda, int *info); void spotrf_ (const char *uplo, const int *n, float *A, const int *lda, int *info); +// Apply forward/backward substitution to Cholesky factorization + void dpotrs_ (const char *uplo, const int *n, const int *nrhs, + const double *A, const int *lda, + double *B, const int *ldb, + int *info); + void spotrs_ (const char *uplo, const int *n, const int *nrhs, + const float *A, const int *lda, + float *B, const int *ldb, + int *info); // Estimate the reciprocal of the condition number in 1-norm from Cholesky void dpocon_ (const char *uplo, const int *n, const double *A, const int *lda, const double *anorm, double *rcond, @@ -819,6 +828,57 @@ getrs (const char *, const int *, const int *, const float *, const int *, const #endif + +/// Template wrapper for LAPACK functions dpotrs and spotrs +template +inline void +potrs (const char *, const int *, const int *, + const number *, const int *, + number *, const int *, + int *) +{ + Assert (false, ExcNotImplemented()); +} + +#ifdef DEAL_II_WITH_LAPACK +inline void +potrs (const char *uplo, const int *n, const int *nrhs, + const double *A, const int *lda, + double *B, const int *ldb, + int *info) +{ + dpotrs_(uplo,n,nrhs,A,lda,B,ldb,info); +} +inline void +potrs (const char *uplo, const int *n, const int *nrhs, + const float *A, const int *lda, + float *B, const int *ldb, + int *info) +{ + spotrs_(uplo,n,nrhs,A,lda,B,ldb,info); +} +#else +inline void +potrs (const char *, const int *, const int *, + const double *, const int *, + double *, const int *, + int *) +{ + Assert (false, LAPACKSupport::ExcMissing("dpotrs")); +} +inline void +potrs (const char *, const int *, const int *, + const float *, const int *, + flaot *, const int *, + int *) +{ + Assert (false, LAPACKSupport::ExcMissing("spotrs")); +} +#endif + + + + /// Template wrapper for LAPACK functions dgetri and sgetri template inline void diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 3cfc34d741..a12e833979 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -888,50 +888,123 @@ LAPACKFullMatrix::invert() } + template void -LAPACKFullMatrix::apply_lu_factorization(Vector &v, - const bool transposed) const +LAPACKFullMatrix::solve(Vector &v, + const bool transposed) const { - Assert(state == lu, ExcState(state)); Assert(this->n_rows() == this->n_cols(), LACExceptions::ExcNotQuadratic()); AssertDimension(this->n_rows(), v.size()); - const char *trans = transposed ? &T : &N; const int nn = this->n_cols(); const number *values = &this->values[0]; + const int n_rhs = 1; int info = 0; - getrs(trans, &nn, &one, values, &nn, ipiv.data(), - v.begin(), &nn, &info); + if (state == lu) + { + getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(), + v.begin(), &nn, &info); + } + else if (state == cholesky) + { + potrs(&LAPACKSupport::L, &nn, &n_rhs, + values, &nn, v.begin(), &nn, &info); + } + else if (property == upper_triangular || + property == lower_triangular) + { + const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L); + + const int lda = nn; + const int ldb = nn; + trtrs (&uplo, trans, "N", + &nn, &n_rhs, + values, &lda, + v.begin(), &ldb, &info); + } + else + { + Assert(false, + ExcMessage("The matrix has to be either factorized or triangular.")); + } Assert(info == 0, ExcInternalError()); } + template void -LAPACKFullMatrix::apply_lu_factorization(LAPACKFullMatrix &B, - const bool transposed) const +LAPACKFullMatrix::solve(LAPACKFullMatrix &B, + const bool transposed) const { - Assert(state == lu, ExcState(state)); - Assert(B.state == matrix, ExcState(state)); - Assert(this->n_rows() == this->n_cols(), LACExceptions::ExcNotQuadratic()); - AssertDimension(this->n_rows(), B.n_rows()); + Assert(B.state == matrix, ExcState(B.state)); + Assert(this->n_rows() == this->n_cols(), + LACExceptions::ExcNotQuadratic()); + AssertDimension(this->n_rows(), B.n_rows()); const char *trans = transposed ? &T : &N; const int nn = this->n_cols(); - const int kk = B.n_cols(); const number *values = &this->values[0]; + const int n_rhs = B.n_cols(); int info = 0; - getrs(trans, &nn, &kk, values, &nn, ipiv.data(), &B.values[0], &nn, &info); + if (state == lu) + { + getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(), + &B.values[0], &nn, &info); + } + else if (state == cholesky) + { + potrs(&LAPACKSupport::L, &nn, &n_rhs, + values, &nn, &B.values[0], &nn, &info); + } + else if (property == upper_triangular || + property == lower_triangular) + { + const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L); + + const int lda = nn; + const int ldb = nn; + trtrs (&uplo, trans, "N", + &nn, &n_rhs, + values, &lda, + &B.values[0], &ldb, &info); + } + else + { + Assert(false, + ExcMessage("The matrix has to be either factorized or triangular.")); + } Assert(info == 0, ExcInternalError()); } + +template +void +LAPACKFullMatrix::apply_lu_factorization(Vector &v, + const bool transposed) const +{ + solve(v,transposed); +} + + + +template +void +LAPACKFullMatrix::apply_lu_factorization(LAPACKFullMatrix &B, + const bool transposed) const +{ + solve(B,transposed); +} + + + template number LAPACKFullMatrix::determinant() const diff --git a/tests/lapack/full_matrix_24.cc b/tests/lapack/full_matrix_24.cc new file mode 100644 index 0000000000..72a63c0e47 --- /dev/null +++ b/tests/lapack/full_matrix_24.cc @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// 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 triangular matrices + +/* MWE for size=3 in Octave: +R = [1,2,3; 0, 5, 6; 0, 0, 9] +x = [2; -7; 1] + +> R\x +ans = + + 4.73333 + -1.53333 + 0.11111 + +> R'\x +ans = + + 2.00000 + -2.20000 + 0.91111 +*/ + +#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::upper_triangular); + + M = 0.; + unsigned int counter = 1; + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = 0; j < size; ++j) + { + if (j >= i) + M(i,j) = counter; + + counter++; + } + + Vector x(size), y(size); + x[0] = 2; + x[1] = -7; + x[2] = 1; + + y = x; + M.solve(y, false); + y.print(deallog.get_file_stream(), 6, false); + + y = x; + M.solve(y, true); + 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_24.output b/tests/lapack/full_matrix_24.output new file mode 100644 index 0000000000..551292b038 --- /dev/null +++ b/tests/lapack/full_matrix_24.output @@ -0,0 +1,3 @@ + +4.733333 -1.533333 0.111111 +2.000000 -2.200000 0.911111 -- 2.39.5