From 93cafc4ae11e910701215e9455acef526d03bd10 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Wed, 25 Oct 2017 13:32:02 +0200
Subject: [PATCH] add LAPACKFullMatrix::add()

---
 doc/news/changes/minor/20171025DenisDavydov |  3 +
 include/deal.II/lac/lapack_full_matrix.h    |  6 ++
 source/lac/lapack_full_matrix.cc            | 22 +++++++
 tests/lapack/full_matrix_21.cc              | 69 +++++++++++++++++++++
 tests/lapack/full_matrix_21.output          |  5 ++
 5 files changed, 105 insertions(+)
 create mode 100644 doc/news/changes/minor/20171025DenisDavydov
 create mode 100644 tests/lapack/full_matrix_21.cc
 create mode 100644 tests/lapack/full_matrix_21.output

diff --git a/doc/news/changes/minor/20171025DenisDavydov b/doc/news/changes/minor/20171025DenisDavydov
new file mode 100644
index 0000000000..736c471604
--- /dev/null
+++ b/doc/news/changes/minor/20171025DenisDavydov
@@ -0,0 +1,3 @@
+New: add LAPACKFullMatrix<number>::add (const number a, const LAPACKFullMatrix<number> &A)
+<br>
+(Denis Davydov, 2017/10/25)
diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h
index ede5e3e48c..8228c23426 100644
--- a/include/deal.II/lac/lapack_full_matrix.h
+++ b/include/deal.II/lac/lapack_full_matrix.h
@@ -134,6 +134,12 @@ public:
   LAPACKFullMatrix<number> &
   operator/= (const number factor);
 
+  /**
+   * Simple addition of a scaled matrix, i.e. <tt>*this += a*A</tt>.
+   */
+  void add (const number                    a,
+            const LAPACKFullMatrix<number> &A);
+
   /**
    * 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 879c906561..c6918d2a3b 100644
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -171,6 +171,28 @@ LAPACKFullMatrix<number>::operator/= (const number factor)
 }
 
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::add (const number              a,
+                               const LAPACKFullMatrix<number> &A)
+{
+  Assert(state == LAPACKSupport::matrix ||
+         state == LAPACKSupport::inverse_matrix,
+         ExcState(state));
+
+  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
+  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
+
+  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);
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::vmult (
diff --git a/tests/lapack/full_matrix_21.cc b/tests/lapack/full_matrix_21.cc
new file mode 100644
index 0000000000..7a2fd4038b
--- /dev/null
+++ b/tests/lapack/full_matrix_21.cc
@@ -0,0 +1,69 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::add() by comparing to FullMatrix
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Full matrix:
+  FullMatrix<NumberType> A(size), B(size);
+  create_random(A);
+  create_random(B);
+
+  // Lapack:
+  LAPACKFullMatrix<NumberType> C(size), D(size);
+  C = A;
+  D = B;
+
+  // do addition
+  NumberType val = 0.1234;
+  A.add(val,B);
+  C.add(val,D);
+
+  FullMatrix<NumberType> diff(size);
+  diff = C;
+  diff.add(-1., A);
+
+  const NumberType error = diff.frobenius_norm();
+  deallog << "difference: " << error << std::endl;
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<unsigned int> sizes = {{17,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_21.output b/tests/lapack/full_matrix_21.output
new file mode 100644
index 0000000000..20782d8779
--- /dev/null
+++ b/tests/lapack/full_matrix_21.output
@@ -0,0 +1,5 @@
+
+DEAL::size=17
+DEAL::difference: 0.00000
+DEAL::size=391
+DEAL::difference: 0.00000
-- 
2.39.5