From c938044e1876ed415dbc9c63a9963272e3022a31 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 8 Nov 2010 00:19:59 +0000 Subject: [PATCH] introduce templates for LAPACK functions git-svn-id: https://svn.dealii.org/trunk@22631 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/common/scripts/lapack_templates.pl | 22 +++++ .../deal.II/lac/full_matrix.templates.h | 89 ++--------------- .../include/deal.II/lac/lapack_full_matrix.h | 98 ++++++++++++++++--- .../lac/precondition_block.templates.h | 8 ++ .../deal.II/lac/precondition_block_base.h | 72 +++++++++++++- .../deal.II/lac/relaxation_block.templates.h | 5 + deal.II/source/lac/lapack_full_matrix.cc | 8 ++ 7 files changed, 204 insertions(+), 98 deletions(-) diff --git a/deal.II/common/scripts/lapack_templates.pl b/deal.II/common/scripts/lapack_templates.pl index 2fd0419614..886e9f9ee5 100644 --- a/deal.II/common/scripts/lapack_templates.pl +++ b/deal.II/common/scripts/lapack_templates.pl @@ -105,7 +105,29 @@ while(<>) $args0 = $args; $args0 =~ s/\*[^,]*,/\*,/g; $args0 =~ s/\*[^,]*$/\*/g; + # First, do the general template None of these functions is + # implemented, but they allow us to link for instance with + # long double lapack support + my $numbers = 1; + my $argst = $args0; + my $typet = $type; + while ($argst =~ s/double/number$numbers/) + { + $numbers++; + } + $typet =~ s/double/number1/g; + $templates .= "\n\n//--------------------------------------------------------------------------------//\ntemplate::equ (const number a, -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 gemm_switch - { - using ::dealii::gemm; - - template - void - gemm (const char*, const char*, const int*, const int*, const int*, - const T*, const T2*, const int*, const T*, const int*, - const T*, T2*, const int*) - { - Assert (false, LAPACKSupport::ExcMissing("dgemm for this data type")); - } - - } -} - - - template template void FullMatrix::mmult (FullMatrix &dst, @@ -616,8 +585,8 @@ void FullMatrix::mmult (FullMatrix &dst, // Use the BLAS function gemm for // calculating the matrix-matrix // product. - internal::gemm_switch::gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m, - &this->values[0], &k, &beta, &dst(0,0), &m); + gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m, + &this->values[0], &k, &beta, &dst(0,0), &m); return; } @@ -704,8 +673,8 @@ void FullMatrix::Tmmult (FullMatrix &dst, // Use the BLAS function gemm for // calculating the matrix-matrix // product. - internal::gemm_switch::gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m, - &this->values[0], &n, &beta, &dst(0,0), &m); + gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m, + &this->values[0], &n, &beta, &dst(0,0), &m); return; } @@ -795,8 +764,8 @@ void FullMatrix::mTmult (FullMatrix &dst, // Use the BLAS function gemm for // calculating the matrix-matrix // product. - internal::gemm_switch::gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k, - &this->values[0], &k, &beta, &dst(0,0), &m); + gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k, + &this->values[0], &k, &beta, &dst(0,0), &m); return; } @@ -880,8 +849,8 @@ void FullMatrix::TmTmult (FullMatrix &dst, // Use the BLAS function gemm for // calculating the matrix-matrix // product. - internal::gemm_switch::gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k, - &this->values[0], &n, &beta, &dst(0,0), &m); + gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k, + &this->values[0], &n, &beta, &dst(0,0), &m); return; } @@ -1714,44 +1683,6 @@ 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 ::dealii::getrf; - - template - void - getrf (const int*, const int*, T*, const int*, int*, int*) - { - Assert (false, LAPACKSupport::ExcMissing("dgetrf for this data type")); - } - - - using ::dealii::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 FullMatrix::gauss_jordan () @@ -1798,7 +1729,7 @@ FullMatrix::gauss_jordan () // Use the LAPACK function getrf for // calculating the LU factorization. - internal::getrf_switch::getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info); + getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info); Assert(info >= 0, ExcInternalError()); Assert(info == 0, LACExceptions::ExcSingular()); @@ -1809,7 +1740,7 @@ FullMatrix::gauss_jordan () // Use the LAPACK function getri for // calculating the actual inverse using // the LU factorization. - internal::getrf_switch::getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info); + getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info); Assert(info >= 0, ExcInternalError()); Assert(info == 0, LACExceptions::ExcSingular()); diff --git a/deal.II/include/deal.II/lac/lapack_full_matrix.h b/deal.II/include/deal.II/lac/lapack_full_matrix.h index d80c651a2d..dfafccc881 100644 --- a/deal.II/include/deal.II/lac/lapack_full_matrix.h +++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h @@ -161,6 +161,7 @@ class LAPACKFullMatrix : public TransposeTable const number factor = 1., const bool transpose = false); + /** * Matrix-vector-multiplication. * @@ -177,19 +178,30 @@ class LAPACKFullMatrix : public TransposeTable * * Source and destination must * not be the same vector. + * + * @note This template only + * exists for compile-time + * compatibility with + * FullMatrix. Implementation is + * only available for VECTOR=Vector<number> */ - void vmult (Vector &w, - const Vector &v, - const bool adding=false) const; + template + void vmult(VECTOR& dst, const VECTOR& src, const bool adding = false) const; /** * Adding Matrix-vector-multiplication. * w += A*v * * Source and destination must * not be the same vector. + * + * @note This template only + * exists for compile-time + * compatibility with + * FullMatrix. Implementation is + * only available for VECTOR=Vector<number> */ - void vmult_add (Vector &w, - const Vector &v) const; + template + void vmult_add (VECTOR& w, const VECTOR& v) const; /** * Transpose @@ -209,9 +221,15 @@ class LAPACKFullMatrix : public TransposeTable * * Source and destination must * not be the same vector. - */ - void Tvmult (Vector &w, - const Vector &v, + * + * @note This template only + * exists for compile-time + * compatibility with + * FullMatrix. Implementation is + * only available for VECTOR=Vector<number> + */ + template + void Tvmult (VECTOR& w, const VECTOR& v, const bool adding=false) const; /** @@ -221,10 +239,26 @@ class LAPACKFullMatrix : public TransposeTable * * Source and destination must * not be the same vector. - */ + * + * @note This template only + * exists for compile-time + * compatibility with + * FullMatrix. Implementation is + * only available for VECTOR=Vector<number> + */ + template + void Tvmult_add (VECTOR& w, const VECTOR& v) const; + + void vmult (Vector &w, + const Vector &v, + const bool adding=false) const; + void vmult_add (Vector &w, + const Vector &v) const; + void Tvmult (Vector &w, + const Vector &v, + const bool adding=false) const; void Tvmult_add (Vector &w, const Vector &v) const; - /** * Compute the LU factorization * of the matrix using LAPACK @@ -576,7 +610,7 @@ class PreconditionLU template template -void +inline void LAPACKFullMatrix::copy_from (const MATRIX& M) { this->reinit (M.m(), M.n()); @@ -592,7 +626,7 @@ LAPACKFullMatrix::copy_from (const MATRIX& M) template template -void +inline void LAPACKFullMatrix::fill ( const MATRIX& M, const unsigned int dst_offset_i, @@ -620,7 +654,43 @@ LAPACKFullMatrix::fill ( template -std::complex +template +inline void +LAPACKFullMatrix::vmult(VECTOR&, const VECTOR&, bool) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +template +inline void +LAPACKFullMatrix::vmult_add(VECTOR&, const VECTOR&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +template +inline void +LAPACKFullMatrix::Tvmult(VECTOR&, const VECTOR&, bool) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +template +inline void +LAPACKFullMatrix::Tvmult_add(VECTOR&, const VECTOR&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +inline std::complex LAPACKFullMatrix::eigenvalue (const unsigned int i) const { Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState()); @@ -633,7 +703,7 @@ LAPACKFullMatrix::eigenvalue (const unsigned int i) const template -number +inline number LAPACKFullMatrix::singular_value (const unsigned int i) const { Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state)); diff --git a/deal.II/include/deal.II/lac/precondition_block.templates.h b/deal.II/include/deal.II/lac/precondition_block.templates.h index 76eb4cacd3..4e28fa9d23 100644 --- a/deal.II/include/deal.II/lac/precondition_block.templates.h +++ b/deal.II/include/deal.II/lac/precondition_block.templates.h @@ -138,6 +138,10 @@ void PreconditionBlock::invert_permuted_diagblocks( case PreconditionBlockBase::householder: this->inverse_householder(0).initialize(M_cell); break; + case PreconditionBlockBase::svd: + this->inverse_svd(0) = M_cell; + this->inverse_svd(0).compute_inverse_svd(0.); + break; default: Assert(false, ExcNotImplemented()); @@ -189,6 +193,10 @@ void PreconditionBlock::invert_permuted_diagblocks( case PreconditionBlockBase::householder: this->inverse_householder(cell).initialize(M_cell); break; + case PreconditionBlockBase::svd: + this->inverse_svd(cell) = M_cell; + this->inverse_svd(cell).compute_inverse_svd(0.); + break; default: Assert(false, ExcNotImplemented()); diff --git a/deal.II/include/deal.II/lac/precondition_block_base.h b/deal.II/include/deal.II/lac/precondition_block_base.h index d8bc501f1b..0b15a13d43 100644 --- a/deal.II/include/deal.II/lac/precondition_block_base.h +++ b/deal.II/include/deal.II/lac/precondition_block_base.h @@ -20,6 +20,7 @@ #include #include #include +#include #include @@ -68,7 +69,12 @@ class PreconditionBlockBase * Use QR decomposition of * the Householder class. */ - householder + householder, + /** + * Use the singular value + * decomposition of LAPACKFullMatrix. + */ + svd }; /** @@ -187,6 +193,12 @@ class PreconditionBlockBase */ Householder& inverse_householder (unsigned int i); + /** + * Access to the inverse diagonal + * blocks if Inversion is #householder. + */ + LAPACKFullMatrix& inverse_svd (unsigned int i); + /** * Access to the inverse diagonal * blocks. @@ -243,7 +255,8 @@ class PreconditionBlockBase * #gauss_jordan is used. Using * number=float saves * memory in comparison with - * number=double. + * number=double, but + * may introduce numerical instability. */ std::vector > var_inverse_full; @@ -251,15 +264,30 @@ class PreconditionBlockBase * Storage of the inverse * matrices of the diagonal * blocks matrices as - * Householder + * Householder * matrices if Inversion * #householder is used. Using * number=float saves * memory in comparison with - * number=double. + * number=double, but + * may introduce numerical instability. */ std::vector > var_inverse_householder; + /** + * Storage of the inverse + * matrices of the diagonal + * blocks matrices as + * LAPACKFullMatrix + * matrices if Inversion + * #svd is used. Using + * number=float saves + * memory in comparison with + * number=double, but + * may introduce numerical instability. + */ + std::vector > var_inverse_svd; + /** * Storage of the original diagonal blocks. * @@ -310,8 +338,10 @@ PreconditionBlockBase::clear() { if (var_inverse_full.size()!=0) var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end()); - if (var_inverse_full.size()!=0) + if (var_inverse_householder.size()!=0) var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end()); + if (var_inverse_svd.size()!=0) + var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end()); if (var_diagonal.size()!=0) var_diagonal.erase(var_diagonal.begin(), var_diagonal.end()); var_same_diagonal = false; @@ -341,6 +371,10 @@ Inversion method) case householder: var_inverse_householder.resize(1); break; + case svd: + var_inverse_svd.resize(1); + var_inverse_svd[0].reinit(b,b); + break; default: Assert(false, ExcNotImplemented()); } @@ -380,6 +414,13 @@ Inversion method) case householder: var_inverse_householder.resize(n); break; + case svd: + { + std::vector > + tmp(n, LAPACKFullMatrix(b)); + var_inverse_svd.swap (tmp); + break; + } default: Assert(false, ExcNotImplemented()); } @@ -414,6 +455,10 @@ PreconditionBlockBase::inverse_vmult( AssertIndexRange (ii, var_inverse_householder.size()); var_inverse_householder[ii].vmult(dst, src); break; + case svd: + AssertIndexRange (ii, var_inverse_svd.size()); + var_inverse_svd[ii].vmult(dst, src); + break; default: Assert(false, ExcNotImplemented()); } @@ -439,6 +484,10 @@ PreconditionBlockBase::inverse_Tvmult( AssertIndexRange (ii, var_inverse_householder.size()); var_inverse_householder[ii].Tvmult(dst, src); break; + case svd: + AssertIndexRange (ii, var_inverse_svd.size()); + var_inverse_svd[ii].Tvmult(dst, src); + break; default: Assert(false, ExcNotImplemented()); } @@ -499,6 +548,19 @@ PreconditionBlockBase::inverse_householder(unsigned int i) } +template +inline +LAPACKFullMatrix& +PreconditionBlockBase::inverse_svd(unsigned int i) +{ + if (same_diagonal()) + return var_inverse_svd[0]; + + AssertIndexRange (i, var_inverse_svd.size()); + return var_inverse_svd[i]; +} + + template inline FullMatrix& diff --git a/deal.II/include/deal.II/lac/relaxation_block.templates.h b/deal.II/include/deal.II/lac/relaxation_block.templates.h index 003ea60b6d..6298b4bfb3 100644 --- a/deal.II/include/deal.II/lac/relaxation_block.templates.h +++ b/deal.II/include/deal.II/lac/relaxation_block.templates.h @@ -132,6 +132,11 @@ RelaxationBlock::invert_diagblocks () case PreconditionBlockBase::householder: this->inverse_householder(block).initialize(M_cell); break; + case PreconditionBlockBase::svd: + this->inverse_svd(block).reinit(bs, bs); + this->inverse_svd(block) = M_cell; + this->inverse_svd(block).compute_inverse_svd(0.); + break; default: Assert(false, ExcNotImplemented()); } diff --git a/deal.II/source/lac/lapack_full_matrix.cc b/deal.II/source/lac/lapack_full_matrix.cc index 407a63e014..a8ed9c210d 100644 --- a/deal.II/source/lac/lapack_full_matrix.cc +++ b/deal.II/source/lac/lapack_full_matrix.cc @@ -683,6 +683,14 @@ LAPACKFullMatrix::operator = (const FullMatrix& M); template LAPACKFullMatrix & LAPACKFullMatrix::operator = (const FullMatrix& M); +template class LAPACKFullMatrix; +template LAPACKFullMatrix & +LAPACKFullMatrix::operator = (const FullMatrix& M); +template LAPACKFullMatrix & +LAPACKFullMatrix::operator = (const FullMatrix& M); +template LAPACKFullMatrix & +LAPACKFullMatrix::operator = (const FullMatrix& M); + template class PreconditionLU; template class PreconditionLU; -- 2.39.5