From d34a21fe21abe340057fd15330f16626a042df80 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 22 Nov 2020 11:29:55 +0100
Subject: [PATCH] Implement addToDiagonal for scalars (interpreted as 1x1
 matrices)

This is required for dune-istl matrices than use scalars to end
the nesting recursion (and not FieldMatrix), like BCRSMatrix<double>.
---
 CHANGELOG.md                                     |  4 +++-
 dune/matrix-vector/addtodiagonal.hh              |  8 ++++++++
 dune/matrix-vector/test/staticmatrixtoolstest.cc | 10 ++++++++++
 3 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4a5d1c7..c766858 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,8 @@
 # Master (will become release 2.8)
 
-- ...
+- The method `addToDiagonal` can now also be called if the matrix is a scalar number type.
+  This is needed since nowadays scalar entries can end the nesting recursion of dune-istl
+  matrices.
 
 ## Deprecations and removals
 
diff --git a/dune/matrix-vector/addtodiagonal.hh b/dune/matrix-vector/addtodiagonal.hh
index 90ad152..56de2fa 100644
--- a/dune/matrix-vector/addtodiagonal.hh
+++ b/dune/matrix-vector/addtodiagonal.hh
@@ -11,6 +11,14 @@ namespace MatrixVector {
       x[i][i] += a;
   }
 
+  /** \brief Specialization for scalars (to be interpreted as 1x1 matrices) */
+  template <class Matrix>
+  static void addToDiagonal(Matrix& x, const Matrix& a)
+  {
+    static_assert(IsNumber<Matrix>::value, "Only scalars can be treated both as matrix and as number!");
+    x += a;
+  }
+
   /*
     the line
 
diff --git a/dune/matrix-vector/test/staticmatrixtoolstest.cc b/dune/matrix-vector/test/staticmatrixtoolstest.cc
index 368f523..b3566a5 100644
--- a/dune/matrix-vector/test/staticmatrixtoolstest.cc
+++ b/dune/matrix-vector/test/staticmatrixtoolstest.cc
@@ -35,6 +35,16 @@ private:
 
     // SQUARE MATRICES
 
+    // scalars
+    {
+      FieldType scalarMatrix_x = 1,
+                scalarMatrix_check = 4;
+
+      addToDiagonal(scalarMatrix_x, scalar_a);
+
+      passed = passed and
+               myDiff(scalarMatrix_x, scalarMatrix_check) < 1e-12;
+    }
     // case FM
     {
       Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_x = {{1, 2}, {3, 4}},
-- 
GitLab