From 401441dc10ad1d126816f5f30ee0c375b1c85e16 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Wed, 24 Aug 2022 20:33:08 +0200
Subject: [PATCH] Revert "changed: unify matrixblock.hh with downstream
 MatrixBlock.hpp"

---
 .../common/fvbasediscretization.hh            |   1 +
 opm/simulators/linalg/matrixblock.hh          | 306 +++++++++---------
 2 files changed, 146 insertions(+), 161 deletions(-)

diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh
index 2aab24c0e..eee4115a0 100644
--- a/opm/models/discretization/common/fvbasediscretization.hh
+++ b/opm/models/discretization/common/fvbasediscretization.hh
@@ -55,6 +55,7 @@
 #include <opm/models/utils/timer.hh>
 #include <opm/models/utils/timerguard.hh>
 #include <opm/models/io/vtkprimaryvarsmodule.hh>
+#include <opm/simulators/linalg/matrixblock.hh>
 
 #include <opm/material/common/MathToolbox.hpp>
 #include <opm/material/common/Valgrind.hpp>
diff --git a/opm/simulators/linalg/matrixblock.hh b/opm/simulators/linalg/matrixblock.hh
index 4672bcb4a..1f95bf851 100644
--- a/opm/simulators/linalg/matrixblock.hh
+++ b/opm/simulators/linalg/matrixblock.hh
@@ -23,206 +23,190 @@
 #ifndef EWOMS_MATRIX_BLOCK_HH
 #define EWOMS_MATRIX_BLOCK_HH
 
-#include <dune/common/dynmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/version.hh>
-#include <dune/common/typetraits.hh>
-
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/owneroverlapcopy.hh>
 #include <dune/istl/superlu.hh>
 #include <dune/istl/umfpack.hh>
-#include <dune/istl/matrixutils.hh>
 
-#include <limits>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/typetraits.hh>
 
 namespace Opm {
-namespace detail {
+namespace MatrixBlockHelp {
 
 template <typename K, int m, int n>
-static inline void invertMatrix(Dune::FieldMatrix<K,m,n>& matrix)
+static inline void invertMatrix(Dune::FieldMatrix<K, m, n>& matrix)
+{ matrix.invert(); }
+
+template <typename K>
+static inline void invertMatrix(Dune::FieldMatrix<K, 1, 1>& matrix)
 {
-    matrix.invert();
+    matrix[0][0] = 1.0/matrix[0][0];
 }
 
 template <typename K>
-static inline void invertMatrix(Dune::FieldMatrix<K,1,1>& matrix)
+static inline void invertMatrix(Dune::FieldMatrix<K, 2, 2>& matrix)
 {
-    Dune::FieldMatrix<K,1,1> tmp(matrix);
-    Dune::FMatrixHelp::invertMatrix(tmp,matrix);
+    Dune::FieldMatrix<K, 2, 2> tmp(matrix);
+    Dune::FMatrixHelp::invertMatrix(tmp, matrix);
 }
 
 template <typename K>
-static inline void invertMatrix(Dune::FieldMatrix<K,2,2>& matrix)
+static inline void invertMatrix(Dune::FieldMatrix<K, 3, 3>& matrix)
 {
-    Dune::FieldMatrix<K,2,2> tmp(matrix);
-    Dune::FMatrixHelp::invertMatrix(tmp,matrix);
+    Dune::FieldMatrix<K, 3, 3> tmp(matrix);
+    Dune::FMatrixHelp::invertMatrix(tmp, matrix);
 }
 
 template <typename K>
-static inline void invertMatrix(Dune::FieldMatrix<K,3,3>& matrix)
+static inline void invertMatrix(Dune::FieldMatrix<K, 4, 4>& matrix)
 {
-    Dune::FieldMatrix<K,3,3> tmp(matrix);
-    Dune::FMatrixHelp::invertMatrix(tmp,matrix);
-}
+    Dune::FieldMatrix<K, 4, 4> tmp(matrix);
 
-//! invert 4x4 Matrix without changing the original matrix
-template <template<class K> class Matrix, typename K>
-static inline K invertMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse)
-{
-    inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
-            matrix[1][1] * matrix[2][3] * matrix[3][2] -
-            matrix[2][1] * matrix[1][2] * matrix[3][3] +
-            matrix[2][1] * matrix[1][3] * matrix[3][2] +
-            matrix[3][1] * matrix[1][2] * matrix[2][3] -
-            matrix[3][1] * matrix[1][3] * matrix[2][2];
+    matrix[0][0] =
+        tmp[1][1]*tmp[2][2]*tmp[3][3] -
+        tmp[1][1]*tmp[2][3]*tmp[3][2] -
+        tmp[2][1]*tmp[1][2]*tmp[3][3] +
+        tmp[2][1]*tmp[1][3]*tmp[3][2] +
+        tmp[3][1]*tmp[1][2]*tmp[2][3] -
+        tmp[3][1]*tmp[1][3]*tmp[2][2];
 
-    inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
-            matrix[1][0] * matrix[2][3] * matrix[3][2] +
-            matrix[2][0] * matrix[1][2] * matrix[3][3] -
-            matrix[2][0] * matrix[1][3] * matrix[3][2] -
-            matrix[3][0] * matrix[1][2] * matrix[2][3] +
-            matrix[3][0] * matrix[1][3] * matrix[2][2];
+    matrix[1][0] =
+       -tmp[1][0]*tmp[2][2]*tmp[3][3] +
+        tmp[1][0]*tmp[2][3]*tmp[3][2] +
+        tmp[2][0]*tmp[1][2]*tmp[3][3] -
+        tmp[2][0]*tmp[1][3]*tmp[3][2] -
+        tmp[3][0]*tmp[1][2]*tmp[2][3] +
+        tmp[3][0]*tmp[1][3]*tmp[2][2];
 
-    inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
-            matrix[1][0] * matrix[2][3] * matrix[3][1] -
-            matrix[2][0] * matrix[1][1] * matrix[3][3] +
-            matrix[2][0] * matrix[1][3] * matrix[3][1] +
-            matrix[3][0] * matrix[1][1] * matrix[2][3] -
-            matrix[3][0] * matrix[1][3] * matrix[2][1];
+    matrix[2][0] =
+        tmp[1][0]*tmp[2][1]*tmp[3][3] -
+        tmp[1][0]*tmp[2][3]*tmp[3][1] -
+        tmp[2][0]*tmp[1][1]*tmp[3][3] +
+        tmp[2][0]*tmp[1][3]*tmp[3][1] +
+        tmp[3][0]*tmp[1][1]*tmp[2][3] -
+        tmp[3][0]*tmp[1][3]*tmp[2][1];
 
-    inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
-            matrix[1][0] * matrix[2][2] * matrix[3][1] +
-            matrix[2][0] * matrix[1][1] * matrix[3][2] -
-            matrix[2][0] * matrix[1][2] * matrix[3][1] -
-            matrix[3][0] * matrix[1][1] * matrix[2][2] +
-            matrix[3][0] * matrix[1][2] * matrix[2][1];
+    matrix[3][0] =
+       -tmp[1][0]*tmp[2][1]*tmp[3][2] +
+        tmp[1][0]*tmp[2][2]*tmp[3][1] +
+        tmp[2][0]*tmp[1][1]*tmp[3][2] -
+        tmp[2][0]*tmp[1][2]*tmp[3][1] -
+        tmp[3][0]*tmp[1][1]*tmp[2][2] +
+        tmp[3][0]*tmp[1][2]*tmp[2][1];
 
-    inverse[0][1]= -matrix[0][1]  * matrix[2][2] * matrix[3][3] +
-            matrix[0][1] * matrix[2][3] * matrix[3][2] +
-            matrix[2][1] * matrix[0][2] * matrix[3][3] -
-            matrix[2][1] * matrix[0][3] * matrix[3][2] -
-            matrix[3][1] * matrix[0][2] * matrix[2][3] +
-            matrix[3][1] * matrix[0][3] * matrix[2][2];
+    matrix[0][1] =
+       -tmp[0][1]*tmp[2][2]*tmp[3][3] +
+        tmp[0][1]*tmp[2][3]*tmp[3][2] +
+        tmp[2][1]*tmp[0][2]*tmp[3][3] -
+        tmp[2][1]*tmp[0][3]*tmp[3][2] -
+        tmp[3][1]*tmp[0][2]*tmp[2][3] +
+        tmp[3][1]*tmp[0][3]*tmp[2][2];
 
-    inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
-            matrix[0][0] * matrix[2][3] * matrix[3][2] -
-            matrix[2][0] * matrix[0][2] * matrix[3][3] +
-            matrix[2][0] * matrix[0][3] * matrix[3][2] +
-            matrix[3][0] * matrix[0][2] * matrix[2][3] -
-            matrix[3][0] * matrix[0][3] * matrix[2][2];
+    matrix[1][1] =
+        tmp[0][0]*tmp[2][2]*tmp[3][3] -
+        tmp[0][0]*tmp[2][3]*tmp[3][2] -
+        tmp[2][0]*tmp[0][2]*tmp[3][3] +
+        tmp[2][0]*tmp[0][3]*tmp[3][2] +
+        tmp[3][0]*tmp[0][2]*tmp[2][3] -
+        tmp[3][0]*tmp[0][3]*tmp[2][2];
 
-    inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
-            matrix[0][0] * matrix[2][3] * matrix[3][1] +
-            matrix[2][0] * matrix[0][1] * matrix[3][3] -
-            matrix[2][0] * matrix[0][3] * matrix[3][1] -
-            matrix[3][0] * matrix[0][1] * matrix[2][3] +
-            matrix[3][0] * matrix[0][3] * matrix[2][1];
+    matrix[2][1] =
+       -tmp[0][0]*tmp[2][1]*tmp[3][3] +
+        tmp[0][0]*tmp[2][3]*tmp[3][1] +
+        tmp[2][0]*tmp[0][1]*tmp[3][3] -
+        tmp[2][0]*tmp[0][3]*tmp[3][1] -
+        tmp[3][0]*tmp[0][1]*tmp[2][3] +
+        tmp[3][0]*tmp[0][3]*tmp[2][1];
 
-    inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
-            matrix[0][0] * matrix[2][2] * matrix[3][1] -
-            matrix[2][0] * matrix[0][1] * matrix[3][2] +
-            matrix[2][0] * matrix[0][2] * matrix[3][1] +
-            matrix[3][0] * matrix[0][1] * matrix[2][2] -
-            matrix[3][0] * matrix[0][2] * matrix[2][1];
+    matrix[3][1] =
+        tmp[0][0]*tmp[2][1]*tmp[3][2] -
+        tmp[0][0]*tmp[2][2]*tmp[3][1] -
+        tmp[2][0]*tmp[0][1]*tmp[3][2] +
+        tmp[2][0]*tmp[0][2]*tmp[3][1] +
+        tmp[3][0]*tmp[0][1]*tmp[2][2] -
+        tmp[3][0]*tmp[0][2]*tmp[2][1];
 
-    inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
-            matrix[0][1] * matrix[1][3] * matrix[3][2] -
-            matrix[1][1] * matrix[0][2] * matrix[3][3] +
-            matrix[1][1] * matrix[0][3] * matrix[3][2] +
-            matrix[3][1] * matrix[0][2] * matrix[1][3] -
-            matrix[3][1] * matrix[0][3] * matrix[1][2];
+    matrix[0][2] =
+        tmp[0][1]*tmp[1][2]*tmp[3][3] -
+        tmp[0][1]*tmp[1][3]*tmp[3][2] -
+        tmp[1][1]*tmp[0][2]*tmp[3][3] +
+        tmp[1][1]*tmp[0][3]*tmp[3][2] +
+        tmp[3][1]*tmp[0][2]*tmp[1][3] -
+        tmp[3][1]*tmp[0][3]*tmp[1][2];
 
-    inverse[1][2] = -matrix[0][0]  * matrix[1][2] * matrix[3][3] +
-            matrix[0][0] * matrix[1][3] * matrix[3][2] +
-            matrix[1][0] * matrix[0][2] * matrix[3][3] -
-            matrix[1][0] * matrix[0][3] * matrix[3][2] -
-            matrix[3][0] * matrix[0][2] * matrix[1][3] +
-            matrix[3][0] * matrix[0][3] * matrix[1][2];
+    matrix[1][2] =
+       -tmp[0][0] *tmp[1][2]*tmp[3][3] +
+        tmp[0][0]*tmp[1][3]*tmp[3][2] +
+        tmp[1][0]*tmp[0][2]*tmp[3][3] -
+        tmp[1][0]*tmp[0][3]*tmp[3][2] -
+        tmp[3][0]*tmp[0][2]*tmp[1][3] +
+        tmp[3][0]*tmp[0][3]*tmp[1][2];
 
-    inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
-            matrix[0][0] * matrix[1][3] * matrix[3][1] -
-            matrix[1][0] * matrix[0][1] * matrix[3][3] +
-            matrix[1][0] * matrix[0][3] * matrix[3][1] +
-            matrix[3][0] * matrix[0][1] * matrix[1][3] -
-            matrix[3][0] * matrix[0][3] * matrix[1][1];
+    matrix[2][2] =
+        tmp[0][0]*tmp[1][1]*tmp[3][3] -
+        tmp[0][0]*tmp[1][3]*tmp[3][1] -
+        tmp[1][0]*tmp[0][1]*tmp[3][3] +
+        tmp[1][0]*tmp[0][3]*tmp[3][1] +
+        tmp[3][0]*tmp[0][1]*tmp[1][3] -
+        tmp[3][0]*tmp[0][3]*tmp[1][1];
 
-    inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
-            matrix[0][0] * matrix[1][2] * matrix[3][1] +
-            matrix[1][0] * matrix[0][1] * matrix[3][2] -
-            matrix[1][0] * matrix[0][2] * matrix[3][1] -
-            matrix[3][0] * matrix[0][1] * matrix[1][2] +
-            matrix[3][0] * matrix[0][2] * matrix[1][1];
+    matrix[3][2] =
+       -tmp[0][0]*tmp[1][1]*tmp[3][2] +
+        tmp[0][0]*tmp[1][2]*tmp[3][1] +
+        tmp[1][0]*tmp[0][1]*tmp[3][2] -
+        tmp[1][0]*tmp[0][2]*tmp[3][1] -
+        tmp[3][0]*tmp[0][1]*tmp[1][2] +
+        tmp[3][0]*tmp[0][2]*tmp[1][1];
 
-    inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
-            matrix[0][1] * matrix[1][3] * matrix[2][2] +
-            matrix[1][1] * matrix[0][2] * matrix[2][3] -
-            matrix[1][1] * matrix[0][3] * matrix[2][2] -
-            matrix[2][1] * matrix[0][2] * matrix[1][3] +
-            matrix[2][1] * matrix[0][3] * matrix[1][2];
+    matrix[0][3] =
+       -tmp[0][1]*tmp[1][2]*tmp[2][3] +
+        tmp[0][1]*tmp[1][3]*tmp[2][2] +
+        tmp[1][1]*tmp[0][2]*tmp[2][3] -
+        tmp[1][1]*tmp[0][3]*tmp[2][2] -
+        tmp[2][1]*tmp[0][2]*tmp[1][3] +
+        tmp[2][1]*tmp[0][3]*tmp[1][2];
 
-    inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
-            matrix[0][0] * matrix[1][3] * matrix[2][2] -
-            matrix[1][0] * matrix[0][2] * matrix[2][3] +
-            matrix[1][0] * matrix[0][3] * matrix[2][2] +
-            matrix[2][0] * matrix[0][2] * matrix[1][3] -
-            matrix[2][0] * matrix[0][3] * matrix[1][2];
+    matrix[1][3] =
+        tmp[0][0]*tmp[1][2]*tmp[2][3] -
+        tmp[0][0]*tmp[1][3]*tmp[2][2] -
+        tmp[1][0]*tmp[0][2]*tmp[2][3] +
+        tmp[1][0]*tmp[0][3]*tmp[2][2] +
+        tmp[2][0]*tmp[0][2]*tmp[1][3] -
+        tmp[2][0]*tmp[0][3]*tmp[1][2];
 
-    inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
-            matrix[0][0] * matrix[1][3] * matrix[2][1] +
-            matrix[1][0] * matrix[0][1] * matrix[2][3] -
-            matrix[1][0] * matrix[0][3] * matrix[2][1] -
-            matrix[2][0] * matrix[0][1] * matrix[1][3] +
-            matrix[2][0] * matrix[0][3] * matrix[1][1];
+    matrix[2][3] =
+       -tmp[0][0]*tmp[1][1]*tmp[2][3] +
+        tmp[0][0]*tmp[1][3]*tmp[2][1] +
+        tmp[1][0]*tmp[0][1]*tmp[2][3] -
+        tmp[1][0]*tmp[0][3]*tmp[2][1] -
+        tmp[2][0]*tmp[0][1]*tmp[1][3] +
+        tmp[2][0]*tmp[0][3]*tmp[1][1];
 
-    inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
-            matrix[0][0] * matrix[1][2] * matrix[2][1] -
-            matrix[1][0] * matrix[0][1] * matrix[2][2] +
-            matrix[1][0] * matrix[0][2] * matrix[2][1] +
-            matrix[2][0] * matrix[0][1] * matrix[1][2] -
-            matrix[2][0] * matrix[0][2] * matrix[1][1];
+    matrix[3][3] =
+        tmp[0][0]*tmp[1][1]*tmp[2][2] -
+        tmp[0][0]*tmp[1][2]*tmp[2][1] -
+        tmp[1][0]*tmp[0][1]*tmp[2][2] +
+        tmp[1][0]*tmp[0][2]*tmp[2][1] +
+        tmp[2][0]*tmp[0][1]*tmp[1][2] -
+        tmp[2][0]*tmp[0][2]*tmp[1][1];
 
-    K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] +
-            matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
+    K det =
+        tmp[0][0]*matrix[0][0] +
+        tmp[0][1]*matrix[1][0] +
+        tmp[0][2]*matrix[2][0] +
+        tmp[0][3]*matrix[3][0];
 
     // return identity for singular or nearly singular matrices.
-    if (std::abs(det) < 1e-40) {
-        inverse = std::numeric_limits<K>::quiet_NaN();
-        det = 0.0;
-    } else
-      inverse *= 1.0 / det;
-
-    return det;
+    if (std::abs(det) < 1e-40)
+        matrix = std::numeric_limits<K>::quiet_NaN();
+    else
+        matrix *= 1.0/det;
 }
-
-template<class K> using FMat4 = Dune::FieldMatrix<K,4,4>;
-
-template <typename K>
-static inline void invertMatrix(Dune::FieldMatrix<K,4,4>& matrix)
-{
-    FMat4<K> tmp(matrix);
-    invertMatrix4<FMat4>(tmp, matrix);
-}
-
-template <typename K>
-static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
-{
-    // this function is only for 4 X 4 matrix
-    // for 4 X 4 matrix, using the invertMatrix() function above
-    // it is for temporary usage, mainly to reduce the huge burden of testing
-    // what algorithm should be used to invert 4 X 4 matrix will be handled
-    // as a seperate issue
-    if (matrix.rows() == 4) {
-         Dune::DynamicMatrix<K> A = matrix;
-         invertMatrix4(A, matrix);
-         return;
-     }
-
-#if ! DUNE_VERSION_NEWER( DUNE_COMMON, 2, 7 )
-     Dune::FMatrixPrecision<K>::set_singular_limit(1.e-30);
-#endif
-     matrix.invert();
-}
-
-} // namespace detail
+} // namespace MatrixBlockHelp
 
 template <class Scalar, int n, int m>
 class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
@@ -243,7 +227,7 @@ public:
     {}
 
     void invert()
-    { detail::invertMatrix(asBase()); }
+    { Opm::MatrixBlockHelp::invertMatrix(asBase()); }
 
     const BaseType& asBase() const
     { return static_cast<const BaseType&>(*this); }