From 453f4e2e4c2af7c4cfebd4461a85616aa54ffed0 Mon Sep 17 00:00:00 2001 From: Benjamin Brands Date: Tue, 20 Feb 2018 19:41:24 +0100 Subject: [PATCH] LAPACKFullMatrix::add() uses BLAS routine --- doc/news/changes/minor/20180220BenjaminBrands | 4 ++++ include/deal.II/lac/lapack_templates.h | 24 +++++++++---------- source/lac/lapack_full_matrix.cc | 20 +++++++++++----- 3 files changed, 30 insertions(+), 18 deletions(-) create mode 100644 doc/news/changes/minor/20180220BenjaminBrands diff --git a/doc/news/changes/minor/20180220BenjaminBrands b/doc/news/changes/minor/20180220BenjaminBrands new file mode 100644 index 0000000000..e5d09899cb --- /dev/null +++ b/doc/news/changes/minor/20180220BenjaminBrands @@ -0,0 +1,4 @@ +New: LAPACKFullMatrix: make operators *= and /= use Lapack function, +LAPACKFullMatrix::add() uses BLAS 1 routine instead of handwritten loops +
+(Benjamin Brands, 2018/02/20) diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index 9e62ea97a2..d0491cef2f 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -317,25 +317,25 @@ extern "C" 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 dealii::types::blas_int *kl, + const dealii::types::blas_int *ku, const double *cfrom, const double *cto, - const int *m, - const int *n, + const dealii::types::blas_int *m, + const dealii::types::blas_int *n, double *A, - const int *lda, - int *info); + const dealii::types::blas_int *lda, + dealii::types::blas_int *info); void slascl_(const char *type, - const int *kl, - const int *ku, + const dealii::types::blas_int *kl, + const dealii::types::blas_int *ku, const float *cfrom, const float *cto, - const int *m, - const int *n, + const dealii::types::blas_int *m, + const dealii::types::blas_int *n, float *A, - const int *lda, - int *info); + const dealii::types::blas_int *lda, + dealii::types::blas_int *info); } DEAL_II_NAMESPACE_OPEN diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 8e71fd8d33..cd23a761cd 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -195,12 +195,13 @@ LAPACKFullMatrix::operator*= (const number factor) 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) + // 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); + // Negative return value implies a wrong argument. This should be internal. Assert(info >= 0, ExcInternalError()); return *this; @@ -224,12 +225,13 @@ LAPACKFullMatrix::operator/= (const number factor) 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) + // 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); + // Negative return value implies a wrong argument. This should be internal. Assert(info >= 0, ExcInternalError()); return *this; @@ -239,7 +241,7 @@ LAPACKFullMatrix::operator/= (const number factor) template void -LAPACKFullMatrix::add (const number a, +LAPACKFullMatrix::add (const number a, const LAPACKFullMatrix &A) { Assert(state == LAPACKSupport::matrix || @@ -251,9 +253,15 @@ LAPACKFullMatrix::add (const number a, AssertIsFinite(a); - for (size_type i=0; i use BLAS 1 for adding vectors + const types::blas_int n = this->m() * this->n(); + const types::blas_int inc=1; + number *values = &this->values[0]; + const number *values_A = &A.values[0]; + + axpy(&n,&a,values_A,&inc,values,&inc); } -- 2.39.5