From fe21ac951928fac80759c6bb07aa96835fd9e9ac Mon Sep 17 00:00:00 2001
From: Benjamin Brands <benjamin.brands@fau.de>
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 <typename number>
+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<number>::operator*= (const number factor)
          state == LAPACKSupport::inverse_matrix,
          ExcState(state));
 
-  for (size_type column = 0; column<this->n(); ++column)
-    for (size_type row = 0; row<this->m(); ++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<number>::operator/= (const number factor)
   AssertIsFinite(factor);
   Assert (factor != number(0.), ExcZero() );
 
-  for (size_type column = 0; column<this->n(); ++column)
-    for (size_type row = 0; row<this->m(); ++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