From b56c5049bed27d446221af7bd15d6be38ddec620 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 19 Nov 2008 21:29:51 +0000 Subject: [PATCH] Merge the new inversion code into the general template. git-svn-id: https://svn.dealii.org/trunk@17650 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/full_matrix.templates.h | 62 ++++++++- deal.II/lac/source/full_matrix.cc | 124 ------------------ 2 files changed, 59 insertions(+), 127 deletions(-) diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 08f678ce9a..643bc966bc 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -1241,9 +1242,64 @@ FullMatrix::gauss_jordan () { Assert (!this->empty(), ExcEmptyMatrix()); Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic()); - - // Gauss-Jordan-Algorithmus - // cf. Stoer I (4th Edition) p. 153 + + // see if we can use Lapack + // algorithms for this and if the + // type for 'number' works for us: +#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && \ + defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_) + if (types_are_equal::value + || + types_are_equal::value) + { + // In case we have the LAPACK functions + // getrf and getri detected at configure, + // we use these algorithms for inversion + // since they provide better performance + // than the deal.II native functions. + // + // Note that BLAS/LAPACK stores matrix + // elements column-wise (i.e., all values in + // one column, then all in the next, etc.), + // whereas the FullMatrix stores them + // row-wise. + // We ignore that difference, and give our + // row-wise data to LAPACK, + // let LAPACK build the inverse of the + // transpose matrix, and read the result as + // if it were row-wise again. In other words, + // we just got ((A^T)^{-1})^T, which is + // A^{-1}. + + const int nn = this->n(); + ipiv.resize(nn); + int info; + + // Use the LAPACK function getrf for + // calculating the LU factorization. + getrf(&nn, &nn, this->data(), &nn, &ipiv[0], &info); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + inv_work.resize (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); + + Assert(info >= 0, ExcInternalError()); + Assert(info == 0, LACExceptions::ExcSingular()); + + return; + } + +#endif + + // otherwise do it by hand. use the + // Gauss-Jordan-Algorithmus from + // Stoer & Bulirsch I (4th Edition) + // p. 153 const unsigned int N = n(); // first get an estimate of the diff --git a/deal.II/lac/source/full_matrix.cc b/deal.II/lac/source/full_matrix.cc index bf7d74b8b4..0fa3339cdb 100644 --- a/deal.II/lac/source/full_matrix.cc +++ b/deal.II/lac/source/full_matrix.cc @@ -13,133 +13,10 @@ #include -#include #include DEAL_II_NAMESPACE_OPEN - - // Need to explicitly state the Lapack - // inversion since it only works with - // floats and doubles in case LAPACK was - // detected by configure. -#if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_) - -template <> -void -FullMatrix::gauss_jordan () -{ - Assert (!this->empty(), ExcEmptyMatrix()); - Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic()); - - // In case we have the LAPACK functions - // getrf and getri detected at configure, - // we use these algorithms for inversion - // since they provide better performance - // than the deal.II native functions. - // - // Note that BLAS/LAPACK stores matrix - // elements column-wise (i.e., all values in - // one column, then all in the next, etc.), - // whereas the FullMatrix stores them - // row-wise. - // We ignore that difference, and give our - // row-wise data to LAPACK, - // let LAPACK build the inverse of the - // transpose matrix, and read the result as - // if it were row-wise again. In other words, - // we just got ((A^T)^{-1})^T, which is - // A^{-1}. - - const int nn = this->n(); - float* values = const_cast (this->data()); - ipiv.resize(nn); - int info; - - // Use the LAPACK function getrf for - // calculating the LU factorization. - getrf(&nn, &nn, values, &nn, &ipiv[0], &info); - - Assert(info >= 0, ExcInternalError()); - Assert(info == 0, LACExceptions::ExcSingular()); - - inv_work.resize (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); - - Assert(info >= 0, ExcInternalError()); - Assert(info == 0, LACExceptions::ExcSingular()); -} - -template <> -void -FullMatrix::gauss_jordan () -{ - Assert (!this->empty(), ExcEmptyMatrix()); - Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic()); - - // In case we have the LAPACK functions - // getrf and getri detected at configure, - // we use these algorithms for inversion - // since they provide better performance - // than the deal.II native functions. - // - // Note that BLAS/LAPACK stores matrix - // elements column-wise (i.e., all values in - // one column, then all in the next, etc.), - // whereas the FullMatrix stores them - // row-wise. - // We ignore that difference, and give our - // row-wise data to LAPACK, - // let LAPACK build the inverse of the - // transpose matrix, and read the result as - // if it were row-wise again. In other words, - // we just got ((A^T)^{-1})^T, which is - // A^{-1}. - - const int nn = this->n(); - double* values = const_cast (this->data()); - ipiv.resize(nn); - int info; - - // Use the LAPACK function getrf for - // calculating the LU factorization. - getrf(&nn, &nn, values, &nn, &ipiv[0], &info); - - Assert(info >= 0, ExcInternalError()); - Assert(info == 0, LACExceptions::ExcSingular()); - - inv_work.resize (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); - - Assert(info >= 0, ExcInternalError()); - Assert(info == 0, LACExceptions::ExcSingular()); -} - - // ... and now the usual instantiations - // of gauss_jordan() and all the rest. -template void FullMatrix::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); - -#else - -template void FullMatrix::gauss_jordan (); -template void FullMatrix::gauss_jordan (); -template void FullMatrix::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); -template void FullMatrix >::gauss_jordan (); - -#endif - - #include "full_matrix.inst" @@ -173,5 +50,4 @@ TEMPL_OP_EQ(std::complex,std::complex); #undef TEMPL_OP_EQ - DEAL_II_NAMESPACE_CLOSE -- 2.39.5