add singularity check to LU decomposition

NCAR · Oct 20, 2023 · e8b3f7f · e8b3f7f
1 parent 11822d9
commit e8b3f7f
Show file tree

Hide file tree

Showing 9 changed files with 134 additions and 7 deletions.
diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp
@@ -71,7 +71,14 @@ namespace micm
         const std::function<LuDecompositionPolicy(const SparseMatrixPolicy<T>&)> create_lu_decomp);
 
     /// @brief Decompose the matrix into upper and lower triangular matrices
+    /// @param matrix Matrix to decompose into lower and upper triangular matrices
+    /// @param is_singular Flag that is set to true if matrix is singular; false otherwise
     void Factor(const SparseMatrixPolicy<T>& matrix);
+
+    /// @brief Decompose the matrix into upper and lower triangular matrices
+    /// @param matrix Matrix to decompose into lower and upper triangular matrices
+    /// @param is_singular Flag that is set to true if matrix is singular; false otherwise
+    void Factor(const SparseMatrixPolicy<T>& matrix, bool& is_singular);
 
     /// @brief Solve for x in Ax = b
     template<template<class> class MatrixPolicy>

diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl
@@ -110,6 +110,12 @@ namespace micm
     lu_decomp_.template Decompose<T, SparseMatrixPolicy>(matrix, lower_matrix_, upper_matrix_);
   }
 
+  template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
+  inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Factor(const SparseMatrixPolicy<T>& matrix, bool& is_singular)
+  {
+    lu_decomp_.template Decompose<T, SparseMatrixPolicy>(matrix, lower_matrix_, upper_matrix_, is_singular);
+  }
+
   template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
   template<template<class> class MatrixPolicy>
     requires(!VectorizableDense<MatrixPolicy<T>> || !VectorizableSparse<SparseMatrixPolicy<T>>)

diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp
@@ -91,15 +91,28 @@ namespace micm
     /// @param L The lower triangular matrix created by decomposition
     /// @param U The upper triangular matrix created by decomposition
     template<typename T, template<class> class SparseMatrixPolicy>
-    requires(!VectorizableSparse<SparseMatrixPolicy<T>>) void Decompose(
+    void Decompose(
         const SparseMatrixPolicy<T>& A,
         SparseMatrixPolicy<T>& L,
         SparseMatrixPolicy<T>& U) const;
+
+    /// @brief Perform an LU decomposition on a given A matrix
+    /// @param A Sparse matrix to decompose
+    /// @param L The lower triangular matrix created by decomposition
+    /// @param U The upper triangular matrix created by decomposition
+    /// @param is_singular Flag that is set to true if A is singular; false otherwise
+    template<typename T, template<class> class SparseMatrixPolicy>
+    requires(!VectorizableSparse<SparseMatrixPolicy<T>>) void Decompose(
+        const SparseMatrixPolicy<T>& A,
+        SparseMatrixPolicy<T>& L,
+        SparseMatrixPolicy<T>& U,
+        bool& is_singular) const;
     template<typename T, template<class> class SparseMatrixPolicy>
     requires(VectorizableSparse<SparseMatrixPolicy<T>>) void Decompose(
         const SparseMatrixPolicy<T>& A,
         SparseMatrixPolicy<T>& L,
-        SparseMatrixPolicy<T>& U) const;
+        SparseMatrixPolicy<T>& U,
+        bool& is_singular) const;
   };
 
 }  // namespace micm

diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl
@@ -145,9 +145,17 @@ namespace micm
   }
 
   template<typename T, template<class> class SparseMatrixPolicy>
-    requires(!VectorizableSparse<SparseMatrixPolicy<T>>)
   inline void LuDecomposition::Decompose(const SparseMatrixPolicy<T>& A, SparseMatrixPolicy<T>& L, SparseMatrixPolicy<T>& U)
       const
+  {
+    bool is_singular;
+    Decompose<T, SparseMatrixPolicy>(A, L, U, is_singular);
+  }
+
+  template<typename T, template<class> class SparseMatrixPolicy>
+    requires(!VectorizableSparse<SparseMatrixPolicy<T>>)
+  inline void LuDecomposition::Decompose(const SparseMatrixPolicy<T>& A, SparseMatrixPolicy<T>& L, SparseMatrixPolicy<T>& U, bool& is_singular)
+      const
   {
     // Loop over blocks
     for (std::size_t i_block = 0; i_block < A.size(); ++i_block)
@@ -164,6 +172,7 @@ namespace micm
       auto lki_nkj = lki_nkj_.begin();
       auto lkj_uji = lkj_uji_.begin();
       auto uii = uii_.begin();
+      is_singular = false;
       for (auto& inLU : niLU_)
       {
         // Upper trianglur matrix
@@ -189,6 +198,11 @@ namespace micm
             L_vector[lki_nkj->first] -= L_vector[lkj_uji->first] * U_vector[lkj_uji->second];
             ++lkj_uji;
           }
+          if( U_vector[*uii] == 0.0 )
+          {
+            is_singular = true;
+            return;
+          }
           L_vector[lki_nkj->first] /= U_vector[*uii];
           ++lki_nkj;
           ++uii;
@@ -199,10 +213,9 @@ namespace micm
 
   template<typename T, template<class> class SparseMatrixPolicy>
     requires(VectorizableSparse<SparseMatrixPolicy<T>>)
-  inline void LuDecomposition::Decompose(const SparseMatrixPolicy<T>& A, SparseMatrixPolicy<T>& L, SparseMatrixPolicy<T>& U)
+  inline void LuDecomposition::Decompose(const SparseMatrixPolicy<T>& A, SparseMatrixPolicy<T>& L, SparseMatrixPolicy<T>& U, bool& is_singular)
       const
   {
-    const std::size_t n_cells = A.GroupVectorSize();
     // Loop over groups of blocks
     for (std::size_t i_group = 0; i_group < A.NumberOfGroups(A.size()); ++i_group)
     {
@@ -218,6 +231,8 @@ namespace micm
       auto lki_nkj = lki_nkj_.begin();
       auto lkj_uji = lkj_uji_.begin();
       auto uii = uii_.begin();
+      is_singular = false;
+      const std::size_t n_cells = std::min(A.GroupVectorSize(), A.size() - i_group * A.GroupVectorSize());
       for (auto& inLU : niLU_)
       {
         // Upper trianglur matrix
@@ -256,7 +271,13 @@ namespace micm
             ++lkj_uji;
           }
           for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell)
+          {
+            if (U_vector[*uii + i_cell] == 0.0) {
+              is_singular = true;
+              return;
+            }
             L_vector[lki_nkj->first + i_cell] /= U_vector[*uii + i_cell];
+          }
           ++lki_nkj;
           ++uii;
         }

diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp
@@ -84,6 +84,7 @@ namespace micm
 
     size_t number_of_grid_cells_{ 1 };  // Number of grid cells to solve simultaneously
     bool reorder_state_{ true };        // Reorder state during solver construction to minimize LU fill-in
+    bool check_singularity_{ false };   // Check for singular A matrix in linear solve of A x = b
 
     // Print RosenbrockSolverParameters to console
     void print() const;

diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl
@@ -757,12 +757,18 @@ namespace micm
     // From my understanding the fortran do loop would only ever do one iteration and is equivalent to what's below
     SparseMatrixPolicy<double> jacobian = jacobian_;
     uint64_t n_consecutive = 0;
-    singular = true;
+    singular = false;
     while (true)
     {
       double alpha = 1 / (H * gamma);
       AlphaMinusJacobian(jacobian, alpha);
-      linear_solver_.Factor(jacobian);
+      if (parameters_.check_singularity_)
+      {
+        linear_solver_.Factor(jacobian, singular);
+      } else {
+        singular = false;
+        linear_solver_.Factor(jacobian);
+      }
       singular = false;  // TODO This should be evaluated in some way
       stats_.decompositions += 1;
       if (!singular)

diff --git a/test/unit/solver/test_jit_lu_decomposition.cpp b/test/unit/solver/test_jit_lu_decomposition.cpp
@@ -41,6 +41,32 @@ TEST(JitLuDecomposition, DenseMatrixVectorOrdering)
       });
 }
 
+TEST(JitLuDecomposition, SingularMatrixVectorOrdering)
+{
+  auto jit{ micm::JitCompiler::create() };
+  if (auto err = jit.takeError())
+  {
+    llvm::logAllUnhandledErrors(std::move(err), llvm::errs(), "[JIT Error]");
+    EXPECT_TRUE(false);
+  }
+  testSingularMatrix<Group1SparseVectorMatrix, micm::JitLuDecomposition<1>>(
+      [&](const Group1SparseVectorMatrix<double>& matrix) -> micm::JitLuDecomposition<1> {
+        return micm::JitLuDecomposition<1>{ jit.get(), matrix };
+      });
+  testSingularMatrix<Group2SparseVectorMatrix, micm::JitLuDecomposition<2>>(
+      [&](const Group2SparseVectorMatrix<double>& matrix) -> micm::JitLuDecomposition<2> {
+        return micm::JitLuDecomposition<2>{ jit.get(), matrix };
+      });
+  testSingularMatrix<Group3SparseVectorMatrix, micm::JitLuDecomposition<3>>(
+      [&](const Group3SparseVectorMatrix<double>& matrix) -> micm::JitLuDecomposition<3> {
+        return micm::JitLuDecomposition<3>{ jit.get(), matrix };
+      });
+  testSingularMatrix<Group4SparseVectorMatrix, micm::JitLuDecomposition<4>>(
+      [&](const Group4SparseVectorMatrix<double>& matrix) -> micm::JitLuDecomposition<4> {
+        return micm::JitLuDecomposition<4>{ jit.get(), matrix };
+      });
+}
+
 TEST(JitLuDecomposition, RandomMatrixVectorOrdering)
 {
   auto jit{ micm::JitCompiler::create() };

diff --git a/test/unit/solver/test_lu_decomposition.cpp b/test/unit/solver/test_lu_decomposition.cpp
@@ -24,6 +24,12 @@ TEST(LuDecomposition, DenseMatrixStandardOrdering)
       [](const SparseMatrixTest<double>& matrix) -> micm::LuDecomposition { return micm::LuDecomposition{ matrix }; });
 }
 
+TEST(LuDecomposition, SingularMatrixStandardOrdering)
+{
+    testSingularMatrix<SparseMatrixTest, micm::LuDecomposition>(
+        [](const SparseMatrixTest<double>& matrix) -> micm::LuDecomposition { return micm::LuDecomposition{ matrix }; });
+}
+
 TEST(LuDecomposition, RandomMatrixStandardOrdering)
 {
   testRandomMatrix<SparseMatrixTest, micm::LuDecomposition>(
@@ -52,6 +58,22 @@ TEST(LuDecomposition, DenseMatrixVectorOrdering)
       { return micm::LuDecomposition{ matrix }; });
 }
 
+TEST(LuDecomposition, SingluarMatrixVectorOrdering)
+{
+  testSingularMatrix<Group1SparseVectorMatrix, micm::LuDecomposition>(
+      [](const Group1SparseVectorMatrix<double>& matrix) -> micm::LuDecomposition
+      { return micm::LuDecomposition{ matrix }; });
+  testSingularMatrix<Group2SparseVectorMatrix, micm::LuDecomposition>(
+      [](const Group2SparseVectorMatrix<double>& matrix) -> micm::LuDecomposition
+      { return micm::LuDecomposition{ matrix }; });
+  testSingularMatrix<Group3SparseVectorMatrix, micm::LuDecomposition>(
+      [](const Group3SparseVectorMatrix<double>& matrix) -> micm::LuDecomposition
+      { return micm::LuDecomposition{ matrix }; });
+  testSingularMatrix<Group4SparseVectorMatrix, micm::LuDecomposition>(
+      [](const Group4SparseVectorMatrix<double>& matrix) -> micm::LuDecomposition
+      { return micm::LuDecomposition{ matrix }; });
+}
+
 TEST(LuDecomposition, RandomMatrixVectorOrdering)
 {
   testRandomMatrix<Group1SparseVectorMatrix, micm::LuDecomposition>(

diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp
@@ -104,6 +104,31 @@ void testDenseMatrix(const std::function<LuDecompositionPolicy(const SparseMatri
       A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-5); });
 }
 
+template<template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
+void testSingularMatrix(const std::function<LuDecompositionPolicy(const SparseMatrixPolicy<double>&)> create_lu_decomp)
+{
+  SparseMatrixPolicy<double> A = SparseMatrixPolicy<double>(SparseMatrixPolicy<double>::create(2)
+                                                                .initial_value(1.0e-30)
+                                                                .with_element(0, 0)
+                                                                .with_element(0, 1)
+                                                                .with_element(1, 0)
+                                                                .with_element(1, 1));
+
+  A[0][0][0] = 0;
+  A[0][0][1] = 1;
+  A[0][1][0] = 1;
+  A[0][1][1] = 1;
+
+  LuDecompositionPolicy lud = create_lu_decomp(A);
+  auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0E-30);
+  bool is_singular{ false };
+  lud.template Decompose<double, SparseMatrixPolicy>(A, LU.first, LU.second, is_singular);
+  EXPECT_TRUE(is_singular);
+  A[0][0][0] = 12;
+  lud.template Decompose<double, SparseMatrixPolicy>(A, LU.first, LU.second, is_singular);
+  EXPECT_FALSE(is_singular);
+}
+
 template<template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
 void testRandomMatrix(
     const std::function<LuDecompositionPolicy(const SparseMatrixPolicy<double>&)> create_lu_decomp,