Compare the 4x4 cofactor determinant against its own scale

opm/simulators/linalg/matrixblock.hh

@@ -34,6 +34,7 @@
 
 #include <opm/common/Exceptions.hpp>
 
+#include <algorithm>
 #include <limits>
 
 namespace Opm {
@@ -187,26 +188,46 @@ static inline K adjugateMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse)
            matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
 }
 
+//! max_i sum_j |a_ij * b_ji|, unchanged by a -> R*a*C for diagonal R and C.
+//! With b the adjugate of a this is the scale of the determinant; with b the
+//! inverse of a it is that scale divided by the determinant.
+template <template<class K> class Matrix, typename K>
+static inline K crossScale4(const Matrix<K>& a, const Matrix<K>& b)
+{
+    K rowSum[4];
+    for (int i = 0; i < 4; ++i) {
+        rowSum[i] = std::abs(a[i][0] * b[0][i]) + std::abs(a[i][1] * b[1][i])
+                  + std::abs(a[i][2] * b[2][i]) + std::abs(a[i][3] * b[3][i]);
+    }
+
+    return std::max(std::max(rowSum[0], rowSum[1]), std::max(rowSum[2], rowSum[3]));
+}
+
 //! invert 4x4 Matrix without changing the original matrix
 //!
 //! The cofactor expansion is used by default. It is branch free and roughly an
 //! order of magnitude faster than a pivoted LU at this size, which matters
 //! because every diagonal block of the ILU decomposition goes through here.
 //!
-//! Its determinant is however the product of the four column scales, so it
-//! underflows for blocks whose residuals are insensitive to some of the
-//! primary variables - a cell in which a phase is absent or immobile - even
-//! though those blocks are perfectly well invertible. Dividing the adjugate by
-//! such a determinant is unreliable, so fall back to Dune's invert(), which
-//! uses Gaussian elimination with partial pivoting and reports a matrix as
-//! singular only when a pivot is exactly zero. Block sizes 5 and above already
-//! use that routine.
+//! Its determinant is however the product of the four column scales, so an
+//! absolute threshold on it tests the units a block is written in rather than
+//! how close to singular it is. Compare it against its own scale instead, and
+//! fall back to Dune's invert() below that. Dune reports a matrix as singular
+//! only when a pivot is exactly zero, which a merely nearly dependent block
+//! does not produce, so the fallback result is checked by the same measure.
+//! Block sizes 5 and above already use that routine.
+//!
+//! The determinant is returned as the cofactors computed it, including after a
+//! fallback. No caller in tree reads it.
 template <template<class K> class Matrix, typename K>
 static inline K invertMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse)
 {
+    constexpr K singularLimit = K(1e-12);
+
     const K det = adjugateMatrix4<Matrix, K>(matrix, inverse);
 
-    if (std::abs(det) < 1e-40) {
+    // `inverse` still holds the adjugate. Negated, so a NaN determinant pivots.
+    if (!(std::abs(det) > singularLimit * crossScale4<Matrix, K>(matrix, inverse))) {
         inverse = matrix;
         try {
             inverse.invert();
@@ -215,6 +236,10 @@ static inline K invertMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse)
             inverse = std::numeric_limits<K>::quiet_NaN();
             DUNE_THROW(Dune::MatrixBlockError, "Singular matrix block");
         }
+        if (!(crossScale4<Matrix, K>(matrix, inverse) < K(1) / singularLimit)) {
+            inverse = std::numeric_limits<K>::quiet_NaN();
+            DUNE_THROW(Dune::MatrixBlockError, "Singular matrix block");
+        }
     }
     else {
         inverse *= 1.0 / det;

tests/test_invert.cpp

@@ -73,8 +73,8 @@ BOOST_AUTO_TEST_CASE(Invert4x4WithUnderflowingDeterminant)
     // A well conditioned matrix (condition number about 9.5) scaled such that its
     // determinant, which is homogeneous of degree one in every row, falls below the
     // 1e-40 threshold used to decide whether the cofactor determinant can be divided
-    // by. Scaling leaves the matrix just as invertible as it was, so this must be
-    // inverted through the pivoted fallback rather than reported as singular.
+    // by. Scaling leaves the matrix just as invertible as it was, so this must not
+    // be reported as singular.
     BaseType matrix;
     const double base[4][4] = {{2, 1, 0, 0},
                                {1, 2, 1, 0},
@@ -99,3 +99,73 @@ BOOST_AUTO_TEST_CASE(Invert4x4WithUnderflowingDeterminant)
         }
     }
 }
+
+BOOST_AUTO_TEST_CASE(Invert4x4RankThree)
+{
+    using BaseType = Dune::FieldMatrix<double, 4, 4>;
+
+    // Last row is the sum of the first two. The cofactor determinant comes out at
+    // roundoff rather than at zero, and elimination reaches no exactly zero pivot,
+    // so neither an absolute threshold nor the pivoting catches this.
+    const double base[4][4] = {{1,  2,  3,  4},
+                               {5,  6,  7,  8},
+                               {9, 10, 11, 13},
+                               {6,  8, 10, 12}};
+
+    for (const double scale : {1.0, 1e-7}) {
+        BaseType matrix;
+        for (int i = 0; i < 4; ++i) {
+            for (int j = 0; j < 4; ++j) {
+                matrix[i][j] = base[i][j] * scale;
+            }
+        }
+
+        BaseType inverse;
+        BOOST_CHECK_THROW(Opm::detail::invertMatrix4<Opm::detail::FMat4>(matrix, inverse),
+                          Dune::MatrixBlockError);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(Invert4x4DecisionIsScaleInvariant)
+{
+    using BaseType = Dune::FieldMatrix<double, 4, 4>;
+
+    // Scaling rows and columns leaves a block exactly as invertible as it was while
+    // moving its determinant through nearly two hundred orders of magnitude. Powers
+    // of two, so the scaling rounds nothing and the singular block keeps its zero
+    // pivot.
+    const double regular[4][4] = {{2, 1, 0, 0},
+                                  {1, 2, 1, 0},
+                                  {0, 1, 2, 1},
+                                  {0, 0, 1, 2}};
+    const double singular[4][4] = {{1, 5,  9, 13},
+                                   {2, 6, 10, 14},
+                                   {3, 7, 11, 15},
+                                   {4, 8, 12, 16}};
+
+    for (const int rowExp : {-40, -13, 0, 13, 40}) {
+        for (const int colExp : {-40, -13, 0, 13, 40}) {
+            BaseType matrix;
+            BaseType singularMatrix;
+            for (int i = 0; i < 4; ++i) {
+                for (int j = 0; j < 4; ++j) {
+                    const int exponent = rowExp + i + colExp - j;
+                    matrix[i][j] = std::ldexp(regular[i][j], exponent);
+                    singularMatrix[i][j] = std::ldexp(singular[i][j], exponent);
+                }
+            }
+
+            BaseType inverse;
+            BOOST_CHECK_NO_THROW(Opm::detail::invertMatrix4<Opm::detail::FMat4>(matrix, inverse));
+            const BaseType product = matrix.rightmultiply(inverse);
+            for (int i = 0; i < 4; ++i) {
+                for (int j = 0; j < 4; ++j) {
+                    BOOST_CHECK_SMALL(product[i][j] - (i == j ? 1.0 : 0.0), 1e-12);
+                }
+            }
+
+            BOOST_CHECK_THROW(Opm::detail::invertMatrix4<Opm::detail::FMat4>(singularMatrix, inverse),
+                              Dune::MatrixBlockError);
+        }
+    }
+}