From 453f4e2e4c2af7c4cfebd4461a85616aa54ffed0 Mon Sep 17 00:00:00 2001
From: Benjamin Brands <benjamin.brands@fau.de>
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
+<br>
+(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<number>::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<number>::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<number>::operator/= (const number factor)
 
 template <typename number>
 void
-LAPACKFullMatrix<number>::add (const number              a,
+LAPACKFullMatrix<number>::add (const number a,
                                const LAPACKFullMatrix<number> &A)
 {
   Assert(state == LAPACKSupport::matrix ||
@@ -251,9 +253,15 @@ LAPACKFullMatrix<number>::add (const number              a,
 
   AssertIsFinite(a);
 
-  for (size_type i=0; i<m(); ++i)
-    for (size_type j=0; j<n(); ++j)
-      (*this)(i,j) += a * A(i,j);
+  // BLAS does not offer functions to add matrices.
+  // LapackFullMatrix is stored in contiguous array
+  // ==> 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