From 6d055a291895bf267a902c8b4d04a719b57523b5 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 19 Nov 2008 22:04:18 +0000 Subject: [PATCH] Correctly implement gauss_jordan using LAPACK. git-svn-id: https://svn.dealii.org/trunk@17655 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/full_matrix.h | 20 -------- .../lac/include/lac/full_matrix.templates.h | 49 +++++++++++++++++-- 2 files changed, 45 insertions(+), 24 deletions(-) diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index 7a27ce59dc..bf4aaa34df 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -1153,26 +1153,6 @@ class FullMatrix : public Table<2,number> //@} friend class Accessor; - - private: - -#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_) - /** - * The vector storing the - * permutations applied for - * pivoting in the - * LU-factorization. - */ - std::vector ipiv; - - /** - * Workspace for calculating the - * inverse matrix from an LU - * factorization. - */ - std::vector inv_work; -#endif - }; /**@}*/ diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 643bc966bc..a6282b63f4 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -1235,6 +1235,43 @@ FullMatrix::print_formatted ( } +namespace internal +{ + // a namespace into which we import + // the global definitions of the + // LAPACK functions getrf for + // various data types and add + // general templates for all other + // types. this allows to call the + // getrf function for all data + // types but generated an exception + // if something is called for a + // data type not supported by + // LAPACK. + namespace getrf_switch + { + using ::getrf; + + template + void + getrf (const int*, const int*, T*, const int*, int*, int*) + { + Assert (false, LAPACKSupport::ExcMissing("dgetr for this data type")); + } + + + using ::getri; + + template + void + getri (const int*, T*, const int*, int*, T*, const int*, int*) + { + Assert (false, LAPACKSupport::ExcMissing("dgetri for this data type")); + } + } +} + + template void @@ -1272,21 +1309,25 @@ FullMatrix::gauss_jordan () // A^{-1}. const int nn = this->n(); - ipiv.resize(nn); + + // workspace for permutations + std::vector ipiv(nn); int info; // Use the LAPACK function getrf for // calculating the LU factorization. - getrf(&nn, &nn, this->data(), &nn, &ipiv[0], &info); + internal::getrf_switch::getrf(&nn, &nn, this->val, &nn, &ipiv[0], &info); Assert(info >= 0, ExcInternalError()); Assert(info == 0, LACExceptions::ExcSingular()); - inv_work.resize (nn); + // scratch array + std::vector inv_work (nn); + // Use the LAPACK function getri for // calculating the actual inverse using // the LU factorization. - getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info); + internal::getrf_switch::getri(&nn, this->val, &nn, &ipiv[0], &inv_work[0], &nn, &info); Assert(info >= 0, ExcInternalError()); Assert(info == 0, LACExceptions::ExcSingular()); -- 2.39.5