LAPACKFullMatrix::add() uses BLAS routine

diff --git a/doc/news/changes/minor/20180220BenjaminBrands b/doc/news/changes/minor/20180220BenjaminBrands
new file mode 100644 (file)
index 0000000..e5d0989
--- /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 9e62ea97a2e6b5648019db39ed1c5cae7e6d24e2..d0491cef2f07770b26e690fb97cb59f1ed0ee6cc 100644 (file)
--- 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 8e71fd8d33b2e95cf3806bb15a8cd4a449b37e52..cd23a761cd0dd59ca68e9e2511765c0a2f20154c 100644 (file)
--- 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);
 }