From 9f20ba1de75aa79d3e285d9a1f6aee5b4c64460c Mon Sep 17 00:00:00 2001
From: Niklas Fehn <fehn@lnm.mw.tum.de>
Date: Tue, 7 Mar 2017 11:41:42 +0100
Subject: [PATCH] Add operators *= and /= to LAPACKFullMatrix

---
 include/deal.II/lac/lapack_full_matrix.h | 12 +++++
 source/lac/lapack_full_matrix.cc         | 35 ++++++++++++
 tests/lapack/full_matrix_13.cc           | 69 ++++++++++++++++++++++++
 tests/lapack/full_matrix_13.output       |  2 +
 4 files changed, 118 insertions(+)
 create mode 100644 tests/lapack/full_matrix_13.cc
 create mode 100644 tests/lapack/full_matrix_13.output

diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h
index ecf25ee85f..e31861687e 100644
--- a/include/deal.II/lac/lapack_full_matrix.h
+++ b/include/deal.II/lac/lapack_full_matrix.h
@@ -121,6 +121,18 @@ public:
   LAPACKFullMatrix<number> &
   operator = (const double d);
 
+  /**
+   * This operator multiplies all entries by a fixed factor.
+   */
+  LAPACKFullMatrix<number> &
+  operator*= (const number factor);
+
+  /**
+   * This operator divides all entries by a fixed factor.
+   */
+  LAPACKFullMatrix<number> &
+  operator/= (const number factor);
+
   /**
    * Assignment from different matrix classes, performing the usual conversion
    * to the transposed format expected by LAPACK. This assignment operator
diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 29029c4ddc..7b46ef45a8 100644
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -130,6 +130,41 @@ LAPACKFullMatrix<number>::operator = (const double d)
 }
 
 
+template <typename number>
+LAPACKFullMatrix<number> &
+LAPACKFullMatrix<number>::operator*= (const number factor)
+{
+  Assert(state == LAPACKSupport::matrix ||
+         state == LAPACKSupport::inverse_matrix,
+         ExcState(state));
+
+  for (unsigned int column = 0; column<this->n(); ++column)
+    for (unsigned int row = 0; row<this->m(); ++row)
+      (*this)(row,column) *= factor;
+
+  return *this;
+}
+
+
+template <typename number>
+LAPACKFullMatrix<number> &
+LAPACKFullMatrix<number>::operator/= (const number factor)
+{
+  Assert(state == LAPACKSupport::matrix ||
+         state == LAPACKSupport::inverse_matrix,
+         ExcState(state));
+
+  AssertIsFinite(factor);
+  Assert (factor != number(0.), ExcZero() );
+
+  for (unsigned int column = 0; column<this->n(); ++column)
+    for (unsigned int row = 0; row<this->m(); ++row)
+      (*this)(row,column) /= factor;
+
+  return *this;
+}
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::vmult (
diff --git a/tests/lapack/full_matrix_13.cc b/tests/lapack/full_matrix_13.cc
new file mode 100644
index 0000000000..cd649d1449
--- /dev/null
+++ b/tests/lapack/full_matrix_13.cc
@@ -0,0 +1,69 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::operator*= and operator/=
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+  const unsigned int m=7;
+  const unsigned int n=11;
+  LAPACKFullMatrix<double> A(m,n);
+  for (unsigned int i=0; i<m; ++i)
+    for (unsigned int j=0; j<n; ++j)
+      A(i,j) = (double)rand()/RAND_MAX;
+
+  LAPACKFullMatrix<double> A_check(A);
+
+  const double factor = 2.345;
+
+  // multiply by factor
+  A *= factor;
+
+  // test multiplication
+  for (unsigned int i=0; i<m; ++i)
+    for (unsigned int j=0; j<n; ++j)
+      AssertThrow(std::abs(A(i,j)-factor*A_check(i,j))<1.e-12,ExcInternalError());
+
+  // divide by factor
+  A /= factor;
+
+  // test division
+  for (unsigned int i=0; i<m; ++i)
+    for (unsigned int j=0; j<n; ++j)
+      AssertThrow(std::abs(A(i,j)-A_check(i,j))<1.e-12,ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test();
+}
diff --git a/tests/lapack/full_matrix_13.output b/tests/lapack/full_matrix_13.output
new file mode 100644
index 0000000000..0fd8fc12f0
--- /dev/null
+++ b/tests/lapack/full_matrix_13.output
@@ -0,0 +1,2 @@
+
+DEAL::OK
-- 
2.39.5