Added templated blockinversion for C-style arrays

opm/simulators/linalg/MatrixBlock.hpp

@@ -537,6 +537,214 @@ namespace Detail
         }
     }
 
+    //! perform out of place matrix inversion on C-style arrays
+    //! must have a specified block_size
+    template <int block_size>
+    struct Inverter
+    {
+        template <typename K>
+        void operator()(const K *matrix, K *inverse)
+        {
+            throw std::logic_error("Not implemented");
+        }
+    };
+
+    //! perform out of place matrix inversion on C-style arrays
+    template <>
+    struct Inverter<4>
+    {
+        template <typename K>
+        void operator()(const K *matrix, K *inverse)
+        {
+            // based on Dune::FMatrixHelp::invertMatrix
+            inverse[0] = matrix[5] * matrix[10] * matrix[15] -
+                    matrix[5] * matrix[11] * matrix[14] -
+                    matrix[9] * matrix[6] * matrix[15] +
+                    matrix[9] * matrix[7] * matrix[14] +
+                    matrix[13] * matrix[6] * matrix[11] -
+                    matrix[13] * matrix[7] * matrix[10];
+
+            inverse[4] = -matrix[4] * matrix[10] * matrix[15] +
+                    matrix[4] * matrix[11] * matrix[14] +
+                    matrix[8] * matrix[6] * matrix[15] -
+                    matrix[8] * matrix[7] * matrix[14] -
+                    matrix[12] * matrix[6] * matrix[11] +
+                    matrix[12] * matrix[7] * matrix[10];
+
+            inverse[8] = matrix[4] * matrix[9] * matrix[15] -
+                    matrix[4] * matrix[11] * matrix[13] -
+                    matrix[8] * matrix[5] * matrix[15] +
+                    matrix[8] * matrix[7] * matrix[13] +
+                    matrix[12] * matrix[5] * matrix[11] -
+                    matrix[12] * matrix[7] * matrix[9];
+
+            inverse[12] = -matrix[4] * matrix[9] * matrix[14] +
+                    matrix[4] * matrix[10] * matrix[13] +
+                    matrix[8] * matrix[5] * matrix[14] -
+                    matrix[8] * matrix[6] * matrix[13] -
+                    matrix[12] * matrix[5] * matrix[10] +
+                    matrix[12] * matrix[6] * matrix[9];
+
+            inverse[1]= -matrix[1]  * matrix[10] * matrix[15] +
+                    matrix[1] * matrix[11] * matrix[14] +
+                    matrix[9] * matrix[2] * matrix[15] -
+                    matrix[9] * matrix[3] * matrix[14] -
+                    matrix[13] * matrix[2] * matrix[11] +
+                    matrix[13] * matrix[3] * matrix[10];
+
+            inverse[5] = matrix[0] * matrix[10] * matrix[15] -
+                    matrix[0] * matrix[11] * matrix[14] -
+                    matrix[8] * matrix[2] * matrix[15] +
+                    matrix[8] * matrix[3] * matrix[14] +
+                    matrix[12] * matrix[2] * matrix[11] -
+                    matrix[12] * matrix[3] * matrix[10];
+
+            inverse[9] = -matrix[0] * matrix[9] * matrix[15] +
+                    matrix[0] * matrix[11] * matrix[13] +
+                    matrix[8] * matrix[1] * matrix[15] -
+                    matrix[8] * matrix[3] * matrix[13] -
+                    matrix[12] * matrix[1] * matrix[11] +
+                    matrix[12] * matrix[3] * matrix[9];
+
+            inverse[13] = matrix[0] * matrix[9] * matrix[14] -
+                    matrix[0] * matrix[10] * matrix[13] -
+                    matrix[8] * matrix[1] * matrix[14] +
+                    matrix[8] * matrix[2] * matrix[13] +
+                    matrix[12] * matrix[1] * matrix[10] -
+                    matrix[12] * matrix[2] * matrix[9];
+
+            inverse[2] = matrix[1] * matrix[6] * matrix[15] -
+                    matrix[1] * matrix[7] * matrix[14] -
+                    matrix[5] * matrix[2] * matrix[15] +
+                    matrix[5] * matrix[3] * matrix[14] +
+                    matrix[13] * matrix[2] * matrix[7] -
+                    matrix[13] * matrix[3] * matrix[6];
+
+            inverse[6] = -matrix[0]  * matrix[6] * matrix[15] +
+                    matrix[0] * matrix[7] * matrix[14] +
+                    matrix[4] * matrix[2] * matrix[15] -
+                    matrix[4] * matrix[3] * matrix[14] -
+                    matrix[12] * matrix[2] * matrix[7] +
+                    matrix[12] * matrix[3] * matrix[6];
+
+            inverse[10] = matrix[0] * matrix[5] * matrix[15] -
+                    matrix[0] * matrix[7] * matrix[13] -
+                    matrix[4] * matrix[1] * matrix[15] +
+                    matrix[4] * matrix[3] * matrix[13] +
+                    matrix[12] * matrix[1] * matrix[7] -
+                    matrix[12] * matrix[3] * matrix[5];
+
+            inverse[14] = -matrix[0] * matrix[5] * matrix[14] +
+                    matrix[0] * matrix[6] * matrix[13] +
+                    matrix[4] * matrix[1] * matrix[14] -
+                    matrix[4] * matrix[2] * matrix[13] -
+                    matrix[12] * matrix[1] * matrix[6] +
+                    matrix[12] * matrix[2] * matrix[5];
+
+            inverse[3] = -matrix[1] * matrix[6] * matrix[11] +
+                    matrix[1] * matrix[7] * matrix[10] +
+                    matrix[5] * matrix[2] * matrix[11] -
+                    matrix[5] * matrix[3] * matrix[10] -
+                    matrix[9] * matrix[2] * matrix[7] +
+                    matrix[9] * matrix[3] * matrix[6];
+
+            inverse[7] = matrix[0] * matrix[6] * matrix[11] -
+                    matrix[0] * matrix[7] * matrix[10] -
+                    matrix[4] * matrix[2] * matrix[11] +
+                    matrix[4] * matrix[3] * matrix[10] +
+                    matrix[8] * matrix[2] * matrix[7] -
+                    matrix[8] * matrix[3] * matrix[6];
+
+            inverse[11] = -matrix[0] * matrix[5] * matrix[11] +
+                    matrix[0] * matrix[7] * matrix[9] +
+                    matrix[4] * matrix[1] * matrix[11] -
+                    matrix[4] * matrix[3] * matrix[9] -
+                    matrix[8] * matrix[1] * matrix[7] +
+                    matrix[8] * matrix[3] * matrix[5];
+
+            inverse[15] = matrix[0] * matrix[5] * matrix[10] -
+                    matrix[0] * matrix[6] * matrix[9] -
+                    matrix[4] * matrix[1] * matrix[10] +
+                    matrix[4] * matrix[2] * matrix[9] +
+                    matrix[8] * matrix[1] * matrix[6] -
+                    matrix[8] * matrix[2] * matrix[5];
+
+            K det = matrix[0] * inverse[0] + matrix[1] * inverse[4] +
+                    matrix[2] * inverse[8] + matrix[3] * inverse[12];
+
+            // return identity for singular or nearly singular matrices.
+            if (std::abs(det) < 1e-40) {
+                for (int i = 0; i < 4; ++i){
+                    inverse[4*i + i] = 1.0;
+                }
+            }
+            K inv_det = 1.0 / det;
+
+            for (unsigned int i = 0; i < 4 * 4; ++i) {
+                inverse[i] *= inv_det;
+            }
+        }
+    };
+
+    //! perform out of place matrix inversion on C-style arrays
+    template <>
+    struct Inverter<3>
+    {
+        template <typename K>
+        void operator()(const K *matrix, K *inverse)
+        {
+            // code generated by maple, copied from Dune::DenseMatrix
+            K t4  = matrix[0] * matrix[4];
+            K t6  = matrix[0] * matrix[5];
+            K t8  = matrix[1] * matrix[3];
+            K t10 = matrix[2] * matrix[3];
+            K t12 = matrix[1] * matrix[6];
+            K t14 = matrix[2] * matrix[6];
+
+            K det = (t4 * matrix[8] - t6 * matrix[7] - t8 * matrix[8] +
+                          t10 * matrix[7] + t12 * matrix[5] - t14 * matrix[4]);
+            K t17 = 1.0 / det;
+
+            inverse[0] =  (matrix[4] * matrix[8] - matrix[5] * matrix[7]) * t17;
+            inverse[1] = -(matrix[1] * matrix[8] - matrix[2] * matrix[7]) * t17;
+            inverse[2] =  (matrix[1] * matrix[5] - matrix[2] * matrix[4]) * t17;
+            inverse[3] = -(matrix[3] * matrix[8] - matrix[5] * matrix[6]) * t17;
+            inverse[4] =  (matrix[0] * matrix[8] - t14) * t17;
+            inverse[5] = -(t6 - t10) * t17;
+            inverse[6] =  (matrix[3] * matrix[7] - matrix[4] * matrix[6]) * t17;
+            inverse[7] = -(matrix[0] * matrix[7] - t12) * t17;
+            inverse[8] =  (t4 - t8) * t17;
+        }
+    };
+
+    //! perform out of place matrix inversion on C-style arrays
+    template <>
+    struct Inverter<2>
+    {
+        template <typename K>
+        void operator()(const K *matrix, K *inverse)
+        {
+            // code based on Dune::DenseMatrix
+            K detinv = matrix[0] * matrix[3] - matrix[1] * matrix[2];
+            detinv = 1 / detinv;
+            inverse[0] =  matrix[3] * detinv;
+            inverse[1] = -matrix[1] * detinv;
+            inverse[2] = -matrix[2] * detinv;
+            inverse[3] =  matrix[0] * detinv;
+        }
+    };
+
+    //! perform out of place matrix inversion on C-style arrays
+    template <>
+    struct Inverter<1>
+    {
+        template <typename K>
+        void operator()(const K *matrix, K *inverse)
+        {
+            inverse[0] = 1.0 / matrix[0];
+        }
+    };
+
 } // namespace Detail
 } // namespace Opm
 

opm/simulators/linalg/bda/BILU0.cpp

@@ -20,6 +20,7 @@
 #include <config.h>
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/simulators/linalg/MatrixBlock.hpp>
 #include <dune/common/timer.hh>
 
 #include <opm/simulators/linalg/bda/BdaSolver.hpp>
@@ -170,6 +171,7 @@ namespace bda
         double pivot[bs * bs];
 
         int LSize = 0;
+        Opm::Detail::Inverter<bs> inverter;   // reuse inverter to invert blocks
 
         Timer t_decomposition;
 
@@ -219,7 +221,7 @@ namespace bda
                 }
             }
             // store the inverse in the diagonal!
-            blockInvert3x3(LUmat->nnzValues + ij * bs * bs, invDiagVals + i * bs * bs);
+            inverter(LUmat->nnzValues + ij * bs * bs, invDiagVals + i * bs * bs);
 
             memcpy(LUmat->nnzValues + ij * bs * bs, invDiagVals + i * bs * bs, sizeof(double) * bs * bs);
         }

opm/simulators/linalg/bda/BlockedMatrix.cpp

@@ -128,33 +128,6 @@ void blockMult(double *mat1, double *mat2, double *resMat) {
 }
 
 
-/* Calculate the inverse of a block. This function is specific for only 3x3 block size.*/
-
-void blockInvert3x3(double *mat, double *res) {
-    // code generated by maple, copied from DUNE
-    double t4  = mat[0] * mat[4];
-    double t6  = mat[0] * mat[5];
-    double t8  = mat[1] * mat[3];
-    double t10 = mat[2] * mat[3];
-    double t12 = mat[1] * mat[6];
-    double t14 = mat[2] * mat[6];
-
-    double det = (t4 * mat[8] - t6 * mat[7] - t8 * mat[8] +
-                  t10 * mat[7] + t12 * mat[5] - t14 * mat[4]);
-    double t17 = 1.0 / det;
-
-    res[0] =  (mat[4] * mat[8] - mat[5] * mat[7]) * t17;
-    res[1] = -(mat[1] * mat[8] - mat[2] * mat[7]) * t17;
-    res[2] =  (mat[1] * mat[5] - mat[2] * mat[4]) * t17;
-    res[3] = -(mat[3] * mat[8] - mat[5] * mat[6]) * t17;
-    res[4] =  (mat[0] * mat[8] - t14) * t17;
-    res[5] = -(t6 - t10) * t17;
-    res[6] =  (mat[3] * mat[7] - mat[4] * mat[6]) * t17;
-    res[7] = -(mat[0] * mat[7] - t12) * t17;
-    res[8] =  (t4 - t8) * t17;
-}
-
-
 #define INSTANTIATE_BDA_FUNCTIONS(n)                                 \
 template BlockedMatrix *allocateBlockedMatrix<n>(int Nb, int nnzbs); \
 template void sortBlockedRow<n>(int *, double *, int, int);          \

opm/simulators/linalg/bda/BlockedMatrix.hpp

@@ -70,11 +70,6 @@ void blockMultSub(double *a, double *b, double *c);
 template <unsigned int block_size>
 void blockMult(double *mat1, double *mat2, double *resMat);
 
-/// Calculate the inverse of a block. This function is specific for only 3x3 block size.
-/// \param[in] mat               input block
-/// \param[inout] res            output block
-void blockInvert3x3(double *mat, double *res);
-
 } // end namespace bda
 
 #endif