Remove singular check for LU decomposition and add const qualifier to the CUDA kernels

include/micm/cuda/solver/cuda_lu_decomposition.cuh

@@ -14,8 +14,7 @@ namespace micm
         const CudaMatrixParam& A_param,
         CudaMatrixParam& L_param,
         CudaMatrixParam& U_param,
-        const LuDecomposeParam& devstruct,
-        bool& is_singular);
+        const LuDecomposeParam& devstruct);
 
     /// This is the function that will copy the constant data
     ///   members of class "CudaLuDecomposition" to the device;

include/micm/cuda/solver/cuda_lu_decomposition.hpp

@@ -85,42 +85,21 @@ namespace micm
     /// @param A is the sparse matrix to decompose
     /// @param L is the lower triangular matrix created by decomposition
     /// @param U is the upper triangular matrix created by decomposition
-    /// @param is_singular Flag that is set to true if A is singular; false otherwise
-    template<class SparseMatrixPolicy>
-    requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void Decompose(
-        const SparseMatrixPolicy& A,
-        SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
-
     template<class SparseMatrixPolicy>
     requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
         SparseMatrixPolicy& U) const;
   };
 
-  template<class SparseMatrixPolicy>
-  requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void CudaLuDecomposition::Decompose(
-      const SparseMatrixPolicy& A,
-      SparseMatrixPolicy& L,
-      SparseMatrixPolicy& U,
-      bool& is_singular) const
-  {
-    auto L_param = L.AsDeviceParam();  // we need to update lower matrix so it can't be constant and must be an lvalue
-    auto U_param = U.AsDeviceParam();  // we need to update upper matrix so it can't be constant and must be an lvalue
-    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular);
-  }
-
   template<class SparseMatrixPolicy>
   requires(CudaMatrix<SparseMatrixPolicy>&& VectorizableSparse<SparseMatrixPolicy>) void CudaLuDecomposition::Decompose(
       const SparseMatrixPolicy& A,
       SparseMatrixPolicy& L,
       SparseMatrixPolicy& U) const
   {
-    bool is_singular = false;
     auto L_param = L.AsDeviceParam();  // we need to update lower matrix so it can't be constant and must be an lvalue
     auto U_param = U.AsDeviceParam();  // we need to update upper matrix so it can't be constant and must be an lvalue
-    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular);
+    micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_);
   }
 }  // end of namespace micm

include/micm/cuda/util/cuda_param.hpp

@@ -33,7 +33,6 @@ struct ProcessSetParam
 /// This struct could be allocated on the host or device;
 struct LuDecomposeParam
 {
-  bool* is_singular = nullptr;
   std::pair<std::size_t, std::size_t>* niLU_ = nullptr;
   char* do_aik_ = nullptr;
   std::size_t* aik_ = nullptr;

include/micm/jit/solver/jit_linear_solver.hpp

@@ -42,12 +42,10 @@ namespace micm
 
     /// @brief Decompose the matrix into upper and lower triangular matrices and general JIT functions
     /// @param matrix Matrix that will be factored into lower and upper triangular matrices
-    /// @param is_singular Flag that will be set to true if matrix is singular; false otherwise
     void Factor(
         SparseMatrixPolicy& matrix,
         SparseMatrixPolicy& lower_matrix,
-        SparseMatrixPolicy& upper_matrix,
-        bool& is_singular) const;
+        SparseMatrixPolicy& upper_matrix) const;
 
     /// @brief Solve for x in Ax = b
     template<class MatrixPolicy>

include/micm/jit/solver/jit_linear_solver.inl

@@ -60,10 +60,9 @@ namespace micm
   inline void JitLinearSolver<L, SparseMatrixPolicy, LuDecompositionPolicy>::Factor(
       SparseMatrixPolicy &matrix,
       SparseMatrixPolicy &lower_matrix,
-      SparseMatrixPolicy &upper_matrix,
-      bool &is_singular) const
+      SparseMatrixPolicy &upper_matrix) const
   {
-    LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(matrix, lower_matrix, upper_matrix, is_singular);
+    LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(matrix, lower_matrix, upper_matrix);
   }
 
   template<std::size_t L, class SparseMatrixPolicy, class LuDecompositionPolicy>

include/micm/jit/solver/jit_lu_decomposition.hpp

@@ -39,9 +39,8 @@ namespace micm
     /// @param A Sparse matrix that will be decomposed
     /// @param lower The lower triangular matrix created by decomposition
     /// @param upper The upper triangular matrix created by decomposition
-    /// @param is_singular Flag that will be set to true if A is singular; false otherwise
     template<class SparseMatrixPolicy>
-    void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper, bool &is_singular)
+    void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper)
         const;
 
    private:

include/micm/jit/solver/jit_lu_decomposition.inl

@@ -209,22 +209,12 @@ namespace micm
   void JitLuDecomposition<L>::Decompose(
       const SparseMatrixPolicy &A,
       SparseMatrixPolicy &lower,
-      SparseMatrixPolicy &upper,
-      bool &is_singular) const
+      SparseMatrixPolicy &upper) const
   {
-    is_singular = false;
     decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data());
     for (size_t block = 0; block < A.NumberOfBlocks(); ++block)
     {
       auto diagonals = upper.DiagonalIndices(block);
-      for (const auto &diagonal : diagonals)
-      {
-        if (upper.AsVector()[diagonal] == 0)
-        {
-          is_singular = true;
-          return;
-        }
-      }
     }
   }
 }  // namespace micm

include/micm/solver/backward_euler.inl

@@ -66,8 +66,6 @@ namespace micm
     std::size_t n_successful_integrations = 0;
     std::size_t n_convergence_failures = 0;
 
-    bool singular = false;
-
     auto derived_class_temporary_variables =
         static_cast<BackwardEulerTemporaryVariables<MatrixPolicy>*>(state.temporary_variables_.get());
     auto& Yn = derived_class_temporary_variables->Yn_;
@@ -112,7 +110,7 @@ namespace micm
         // (y_{n+1} - y_n) / H = f(t_{n+1}, y_{n+1})
 
         // try to find the root by factoring and solving the linear system
-        linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular);
+        linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_);
         result.stats_.decompositions_++;
 
         // forcing_blk in camchem

include/micm/solver/linear_solver.hpp

@@ -78,12 +78,10 @@ namespace micm
 
     /// @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& matrix,
         SparseMatrixPolicy& lower_matrix,
-        SparseMatrixPolicy& upper_matrix,
-        bool& is_singular) const;
+        SparseMatrixPolicy& upper_matrix) const;
 
     /// @brief Solve for x in Ax = b. x should be a copy of b and after Solve finishes x will contain the result
     template<class MatrixPolicy>

include/micm/solver/linear_solver.inl

@@ -111,12 +111,11 @@ namespace micm
   inline void LinearSolver<SparseMatrixPolicy, LuDecompositionPolicy>::Factor(
       const SparseMatrixPolicy& matrix,
       SparseMatrixPolicy& lower_matrix,
-      SparseMatrixPolicy& upper_matrix,
-      bool& is_singular) const
+      SparseMatrixPolicy& upper_matrix) const
   {
     MICM_PROFILE_FUNCTION();
 
-    lu_decomp_.template Decompose<SparseMatrixPolicy>(matrix, lower_matrix, upper_matrix, is_singular);
+    lu_decomp_.template Decompose<SparseMatrixPolicy>(matrix, lower_matrix, upper_matrix);
   }
 
   template<class SparseMatrixPolicy, class LuDecompositionPolicy>

include/micm/solver/lu_decomposition.hpp

@@ -122,19 +122,16 @@ namespace micm
     /// @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<class SparseMatrixPolicy>
     requires(!VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
+        SparseMatrixPolicy& U) const;
     template<class SparseMatrixPolicy>
     requires(VectorizableSparse<SparseMatrixPolicy>) void Decompose(
         const SparseMatrixPolicy& A,
         SparseMatrixPolicy& L,
-        SparseMatrixPolicy& U,
-        bool& is_singular) const;
+        SparseMatrixPolicy& U) const;
 
    private:
     /// @brief Initialize arrays for the LU decomposition

include/micm/solver/lu_decomposition.inl

@@ -168,11 +168,9 @@ namespace micm
   requires(!VectorizableSparse<SparseMatrixPolicy>) inline void LuDecomposition::Decompose(
       const SparseMatrixPolicy& A,
       SparseMatrixPolicy& L,
-      SparseMatrixPolicy& U,
-      bool& is_singular) const
+      SparseMatrixPolicy& U) const
   {
     MICM_PROFILE_FUNCTION();
-    is_singular = false;
 
     // Loop over blocks
     for (std::size_t i_block = 0; i_block < A.NumberOfBlocks(); ++i_block)
@@ -226,30 +224,19 @@ 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;
-          }
           L_vector[lki_nkj->first] /= U_vector[*uii];
           ++lki_nkj;
           ++uii;
         }
       }
-      // check the bottom right corner of the matrix
-      if (U_vector[*uii] == 0.0)
-      {
-        is_singular = true;
-      }
     }
   }
 
   template<class SparseMatrixPolicy>
   requires(VectorizableSparse<SparseMatrixPolicy>) inline void LuDecomposition::Decompose(
       const SparseMatrixPolicy& A,
       SparseMatrixPolicy& L,
-      SparseMatrixPolicy& U,
-      bool& is_singular) const
+      SparseMatrixPolicy& U) const
   {
     MICM_PROFILE_FUNCTION();
 
@@ -258,7 +245,6 @@ namespace micm
     const std::size_t A_GroupSizeOfFlatBlockSize = A.GroupSize(A.FlatBlockSize());
     const std::size_t L_GroupSizeOfFlatBlockSize = L.GroupSize(L.FlatBlockSize());
     const std::size_t U_GroupSizeOfFlatBlockSize = U.GroupSize(U.FlatBlockSize());
-    is_singular = false;
 
     // Loop over groups of blocks
     for (std::size_t i_group = 0; i_group < A.NumberOfGroups(A_BlockSize); ++i_group)
@@ -328,29 +314,12 @@ namespace micm
           const std::size_t uii_deref = *uii;
           for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell)
           {
-            if (U_vector[uii_deref + i_cell] == 0.0)
-            {
-              is_singular = true;
-            }
             L_vector[lki_nkj_first + i_cell] /= U_vector[uii_deref + i_cell];
           }
           ++lki_nkj;
           ++uii;
         }
       }
-      const std::size_t uii_deref = *uii;
-      if (n_cells != A_GroupVectorSize)
-      {
-        // check the bottom right corner of the matrix
-        for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell)
-        {
-          if (U_vector[uii_deref + i_cell] == 0.0)
-          {
-            is_singular = true;
-            break;
-          }
-        }
-      }
     }
   }
 

include/micm/solver/rosenbrock.hpp

@@ -105,14 +105,12 @@ namespace micm
     /// @brief Perform the LU decomposition of the matrix
     /// @param H The time step
     /// @param gamma The gamma value
-    /// @param singular A flag to indicate if the matrix is singular
     /// @param number_densities The number densities
     /// @param stats The solver stats
     /// @param state The state
     void LinearFactor(
         double& H,
         const double gamma,
-        bool& singular,
         const auto& number_densities,
         SolverStats& stats,
         auto& state) const;

include/micm/solver/rosenbrock.inl

@@ -69,14 +69,8 @@ namespace micm
       //  Repeat step calculation until current step accepted
       while (!accepted)
       {
-        bool is_singular{ false };
         // Form and factor the rosenbrock ode jacobian
-        LinearFactor(H, parameters_.gamma_[0], is_singular, Y, stats, state);
-        if (is_singular)
-        {
-          result.state_ = SolverState::RepeatedlySingularMatrix;
-          break;
-        }
+        LinearFactor(H, parameters_.gamma_[0], Y, stats, state);
 
         // Compute the stages
         for (uint64_t stage = 0; stage < parameters_.stages_; ++stage)
@@ -230,37 +224,17 @@ namespace micm
   inline void AbstractRosenbrockSolver<RatesPolicy, LinearSolverPolicy, Derived>::LinearFactor(
       double& H,
       const double gamma,
-      bool& singular,
       const auto& number_densities,
       SolverStats& stats,
       auto& state) const
   {
     MICM_PROFILE_FUNCTION();
 
-    uint64_t n_consecutive = 0;
-    singular = false;
-    while (true)
-    {
-      double alpha = 1 / (H * gamma);
-      static_cast<const Derived*>(this)->AlphaMinusJacobian(state.jacobian_, alpha);
-
-      linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular);
-      stats.decompositions_ += 1;
+    double alpha = 1 / (H * gamma);
+    static_cast<const Derived*>(this)->AlphaMinusJacobian(state.jacobian_, alpha);
 
-      // if we are checking for singularity and the matrix is not singular, we can break the loop
-      // if we are not checking for singularity, we always break the loop
-      if (!singular || !parameters_.check_singularity_)
-        break;
-
-      stats.singular_ += 1;
-      if (++n_consecutive > 5)
-        break;
-      H /= 2;
-      // Reconstruct the Jacobian matrix if a substepping is performed here
-      state.jacobian_.Fill(0);
-      rates_.SubtractJacobianTerms(state.rate_constants_, number_densities, state.jacobian_);
-      stats.jacobian_updates_ += 1;
-    }
+    linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_);
+    stats.decompositions_ += 1;
   }
 
   template<class RatesPolicy, class LinearSolverPolicy, class Derived>

include/micm/solver/rosenbrock_solver_parameters.hpp

@@ -63,8 +63,6 @@ namespace micm
     std::vector<double> absolute_tolerance_;
     double relative_tolerance_{ 1e-6 };
 
-    bool check_singularity_{ false };  // Check for singular A matrix in linear solve of A x = b
-
     // Print RosenbrockSolverParameters to console
     void Print() const;
 

include/micm/solver/solver_result.hpp

@@ -45,9 +45,6 @@ namespace micm
     uint64_t decompositions_{};
     /// @brief The number of linear solves
     uint64_t solves_{};
-    /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix
-    /// is never singular
-    uint64_t singular_{};
 
     /// @brief Set all member variables to zero
     void Reset();
@@ -62,7 +59,6 @@ namespace micm
     rejected_ = 0;
     decompositions_ = 0;
     solves_ = 0;
-    singular_ = 0;
   }
 
   inline std::string SolverStateToString(const SolverState& state)