From fe21ac951928fac80759c6bb07aa96835fd9e9ac Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Tue, 20 Feb 2018 18:21:56 +0100 Subject: [PATCH] LAPACKFullMatrix: make operators *= and /= use Lapack function --- include/deal.II/lac/lapack_templates.h | 88 +++++++++++++++++++++++++- source/lac/lapack_full_matrix.cc | 33 ++++++++-- 2 files changed, 114 insertions(+), 7 deletions(-) diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index 45932cdb4a..9e62ea97a2 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -315,7 +315,27 @@ extern "C" const float *alpha, const float *x, const dealii::types::blas_int *incx, float *A, const dealii::types::blas_int *lda); - +// Multiply rectangular mxn real matrix by real scalar CTO/CFROM + void dlascl_(const char *type, + const int *kl, + const int *ku, + const double *cfrom, + const double *cto, + const int *m, + const int *n, + double *A, + const int *lda, + int *info); + void slascl_(const char *type, + const int *kl, + const int *ku, + const float *cfrom, + const float *cto, + const int *m, + const int *n, + float *A, + const int *lda, + int *info); } DEAL_II_NAMESPACE_OPEN @@ -1534,6 +1554,72 @@ stev (const char *, const types::blas_int *, float *, float *, float *, const ty #endif +/// Template wrapper for LAPACK functions dlascl and slascl +template +inline void +lascl(const char *, + const types::blas_int *, + const types::blas_int *, + const number *, + const number *, + const types::blas_int *, + const types::blas_int *, + number *, + const types::blas_int *, + types::blas_int *) +{ + Assert (false, ExcNotImplemented()); +} + +#ifdef DEAL_II_WITH_LAPACK +inline void +lascl(const char *type, + const types::blas_int *kl, + const types::blas_int *ku, + const double *cfrom, + const double *cto, + const types::blas_int *m, + const types::blas_int *n, + double *A, + const types::blas_int *lda, + types::blas_int *info) +{ + dlascl_(type,kl,ku,cfrom,cto,m,n,A,lda,info); +} +#else +inline void +lascl(const char *, const types::blas_int *, const types::blas_int *, const double *, const double *, + const types::blas_int *, const types::blas_int *, double *, const types::blas_int *, types::blas_int *) +{ + Assert (false, LAPACKSupport::ExcMissing("dlascl")); +} +#endif + +#ifdef DEAL_II_WITH_LAPACK +inline void +lascl(const char *type, + const types::blas_int *kl, + const types::blas_int *ku, + const float *cfrom, + const float *cto, + const types::blas_int *m, + const types::blas_int *n, + float *A, + const types::blas_int *lda, + types::blas_int *info) +{ + slascl_(type,kl,ku,cfrom,cto,m,n,A,lda,info); +} +#else +inline void +lascl(const char *, const types::blas_int *, const types::blas_int *, const float *, const float *, + const types::blas_int *, const types::blas_int *, float *, const types::blas_int *, types::blas_int *) +{ + Assert (false, LAPACKSupport::ExcMissing("slascl")); +} +#endif + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index d7fafdb05c..8e71fd8d33 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -188,9 +188,20 @@ LAPACKFullMatrix::operator*= (const number factor) state == LAPACKSupport::inverse_matrix, ExcState(state)); - for (size_type column = 0; columnn(); ++column) - for (size_type row = 0; rowm(); ++row) - (*this)(row,column) *= factor; + AssertIsFinite(factor); + const char type = 'G'; + const number cfrom = 1.; + const types::blas_int m = this->m(); + const types::blas_int n = this->n(); + const types::blas_int lda = this->m(); + types::blas_int info; + // kl and ku will not be referenced for type = G (dense matrices) + const types::blas_int kl=0; + number *values = &this->values[0]; + + lascl(&type,&kl,&kl,&cfrom,&factor,&m,&n,values,&lda,&info); + + Assert(info >= 0, ExcInternalError()); return *this; } @@ -207,9 +218,19 @@ LAPACKFullMatrix::operator/= (const number factor) AssertIsFinite(factor); Assert (factor != number(0.), ExcZero() ); - for (size_type column = 0; columnn(); ++column) - for (size_type row = 0; rowm(); ++row) - (*this)(row,column) /= factor; + const char type = 'G'; + const number cto = 1.; + const types::blas_int m = this->m(); + const types::blas_int n = this->n(); + const types::blas_int lda = this->m(); + types::blas_int info; + // kl and ku will not be referenced for type = G (dense matrices) + const types::blas_int kl=0; + number *values = &this->values[0]; + + lascl(&type,&kl,&kl,&factor,&cto,&m,&n,values,&lda,&info); + + Assert(info >= 0, ExcInternalError()); return *this; } -- 2.39.5