CholmodSolver: Allow to reuse a factorization

Computing Cholesky factorizations is expensive, and they absolutely
have to be reused if the same linear system is solved with several
different right hand sides.  The previous implementation of
CholmodSolver wouldn't allow this: The 'solve' method factorized
and solved together.

This patch makes the interface of the CholmodSolver class richer
to allow for factorization reuse.  It is now possible to set the
matrix, rhs and solution storage separately.  The matrix can be
factorized by calling a new 'factorize' method, and if 'solve'
is called afterwards, the matrix won't be factorized again.

On the other hand, when calling 'solve' by itself the 'solve'
method will do the factorization.  In other words, the patch
is fully backward-compatible.

CHANGELOG.md

@@ -6,6 +6,9 @@
 
 - CholmodSolver: `CholmodSolver` can be used with any blocked matrix/vector type supported by `flatMatrixForEach` and `flatVectorForEach` in dune-istl.
 
+- CholmodSolver: There are now various new methods that allow to factorize the matrix once,
+  and use the factorization to solve linear systems with several right hand sides.
+
 ## Deprecations and removals
 
 - The file `blockgsstep.hh` has been removed after four years of deprecation.  Use the replacement

dune/solvers/solvers/cholmodsolver.hh

@@ -46,7 +46,14 @@ public:
   /** \brief Default constructor */
   CholmodSolver ()
   : LinearSolver<MatrixType,VectorType>(NumProc::FULL)
-  {}
+  {
+    // Bug in LibScotchMetis:
+    // metis_dswitch: see cholmod.h: LibScotchMetis throws segmentation faults for large problems.
+    // This caused a hard to understand CI problems and therefore, METIS is disabled for problems
+    // larger than defined in metis_nswitch for now (default 3000).
+    // Bug report: https://gitlab.inria.fr/scotch/scotch/issues/8
+    cholmodCommonObject().metis_dswitch = 0.0;
+  }
 
   /** \brief Constructor for a linear problem */
   CholmodSolver (const MatrixType& matrix,
@@ -57,56 +64,117 @@ public:
     matrix_(&matrix),
     x_(&x),
     rhs_(&rhs)
-  {}
+  {
+    // Bug in LibScotchMetis:
+    // metis_dswitch: see cholmod.h: LibScotchMetis throws segmentation faults for large problems.
+    // This caused a hard to understand CI problems and therefore, METIS is disabled for problems
+    // larger than defined in metis_nswitch for now (default 3000).
+    // Bug report: https://gitlab.inria.fr/scotch/scotch/issues/8
+    cholmodCommonObject().metis_dswitch = 0.0;
+  }
 
   /** \brief Set the linear problem to solve
    *
-   * Warning! The matrix is assumed to be symmetric!
-   * CHOLMOD uses only the upper part of the matrix, the lower part is ignored and can be empty for optimal memory usage.
+   * The matrix is assumed to be symmetric! CHOLMOD uses only the
+   * upper part of the matrix, the lower part is ignored and can be
+   * empty for optimal memory use.
    */
   void setProblem(const MatrixType& matrix,
                   VectorType& x,
                   const VectorType& rhs) override
+  {
+    setMatrix(matrix);
+    setSolutionVector(x);
+    setRhs(rhs);
+  }
+
+  /** \brief Set the matrix for the linear linear problem to solve
+   *
+   * The matrix is assumed to be symmetric! CHOLMOD uses only the
+   * upper part of the matrix, the lower part is ignored and can be
+   * empty for optimal memory usage.
+   *
+   * \warning Unlike the setMatrix method of the dune-istl Cholmod class,
+   * the method here does not factorize the matrix!   Call the factorize
+   * method for that.
+   */
+  void setMatrix(const MatrixType& matrix)
   {
     matrix_ = &matrix;
-    x_ = &x;
+    isFactorized_ = false;
+  }
+
+  /** \brief Set the rhs vector for the linear problem
+   */
+  void setRhs(const VectorType& rhs)
+  {
     rhs_ = &rhs;
   }
 
-  virtual void solve()
+  /** \brief Set the vector where to write the solution of the linear problem
+   */
+  void setSolutionVector(VectorType& x)
+  {
+    x_ = &x;
+  }
+
+  /** \brief Compute the Cholesky decomposition of the matrix
+   *
+   * You must have set a matrix to be able to call this,
+   * either by calling setProblem or the appropriate constructor.
+   *
+   * The degrees of freedom to ignore must have been set before
+   * calling the `factorize` method.
+   */
+  void factorize()
+  {
+    if (not this->hasIgnore())
+    {
+      // setMatrix will decompose the matrix
+      istlCholmodSolver_.setMatrix(*matrix_);
+
+      // get the error code from Cholmod in case something happened
+      errorCode_ = cholmodCommonObject().status;
+    }
+    else
+    {
+      const auto& ignore = this->ignore();
+
+      // The setMatrix method stores only the selected entries of the matrix (determined by the ignore field)
+      istlCholmodSolver_.setMatrix(*matrix_,&ignore);
+      // get the error code from Cholmod in case something happened
+      errorCode_ = cholmodCommonObject().status;
+    }
+
+    if (errorCode_ >= 0)  // okay or warning, but no error
+      isFactorized_ = true;
+  }
+
+  /** \brief Solve the linear system
+   *
+   * This factorizes the matrix if it hasn't been factorized earlier
+   */
+  virtual void solve() override
   {
     // prepare the solver
     Dune::InverseOperatorResult statistics;
-    Dune::Cholmod<VectorType> solver;
 
-    // Bug in LibScotchMetis:
-    // metis_dswitch: see cholmod.h: LibScotchMetis throws segmentation faults for large problems.
-    // This caused a hard to understand CI problems and therefore, METIS is disabled for problems
-    // larger than defined in metis_nswitch for now (default 3000).
-    // Bug report: https://gitlab.inria.fr/scotch/scotch/issues/8
-    solver.cholmodCommonObject().metis_dswitch = 0.0;
+    // Factorize the matrix now, if the caller hasn't done so explicitly earlier.
+    if (!isFactorized_)
+      factorize();
 
     if (not this->hasIgnore())
     {
       // the apply() method doesn't like constant values
       auto mutableRhs = *rhs_;
 
-      // setMatrix will decompose the matrix
-      solver.setMatrix(*matrix_);
-      // get the error code from Cholmod in case something happened
-      errorCode_ = solver.cholmodCommonObject().status;
-
-      solver.apply(*x_, mutableRhs, statistics);
+      // The back-substitution
+      istlCholmodSolver_.apply(*x_, mutableRhs, statistics);
     }
     else
     {
       const auto& ignore = this->ignore();
 
-      // The setMatrix method stores only the selected entries of the matrix (determined by the ignore field)
-      solver.setMatrix(*matrix_,&ignore);
-      // get the error code from Cholmod in case something happened
-      errorCode_ = solver.cholmodCommonObject().status;
-
       // create flat vectors
       using T = typename VectorType::field_type;
       std::size_t flatSize = flatVectorForEach(ignore, [](...){});
@@ -148,10 +216,15 @@ public:
       });
 
       // Solve the modified system
-      solver.apply(*x_, modifiedRhs, statistics);
+      istlCholmodSolver_.apply(*x_, modifiedRhs, statistics);
     }
   }
 
+  cholmod_common& cholmodCommonObject()
+  {
+    return istlCholmodSolver_.cholmodCommonObject();
+  }
+
   /** \brief return the error code of the Cholmod factorize call
    *
    * In setMatrix() the matrix factorization takes place and Cholmod is
@@ -165,6 +238,9 @@ public:
     return errorCode_;
   }
 
+  //! dune-istl Cholmod solver object
+  Dune::Cholmod<VectorType> istlCholmodSolver_;
+
   //! Pointer to the system matrix
   const MatrixType* matrix_;
 
@@ -176,6 +252,13 @@ public:
 
   //! error code of Cholmod factorize call
   int errorCode_ = 0;
+
+  /** \brief Whether istlCholmodSolver_ currently holds a factorized matrix.
+   *
+   * isFactorized==true is equivalent to istlCholmodSolver_->L_ != nullptr,
+   * but the latter information is private to the dune-istl Cholmod class.
+   */
+  bool isFactorized_ = false;
 };
 
 }   // namespace Solvers

dune/solvers/test/cholmodsolvertest.cc

@@ -59,7 +59,80 @@ struct CHOLMODSolverTestSuite
   }
 };
 
+// Check whether we can reuse a factorization for several right hand sides
+bool checkMultipleSolves()
+{
+  bool passed = true;
+
+  using Matrix = BCRSMatrix<double>;
+  using Vector = BlockVector<double>;
+
+  // Construct an example matrix
+  Matrix matrix(4, 4,
+                3,    // Expected average number of nonzeros
+                      // per row
+                0.2,  // Size of the compression buffer,
+                      // as a fraction of the expected total number of nonzeros
+                BCRSMatrix<double>::implicit);
+
+  matrix.entry(0,0) = 2;
+  matrix.entry(0,1) = 1;
+  matrix.entry(1,0) = 1;
+  matrix.entry(1,1) = 2;
+  matrix.entry(1,2) = 1;
+  matrix.entry(2,1) = 1;
+  matrix.entry(2,2) = 2;
+  matrix.entry(2,3) = 1;
+  matrix.entry(3,2) = 1;
+  matrix.entry(3,3) = 2;
+
+  matrix.compress();
+
+  // Set up the solver
+  Solvers::CholmodSolver<Matrix,Vector> cholmodSolver;
+
+  // Where to write the solution
+  Vector x(4);
+  cholmodSolver.setSolutionVector(x);
+
+  // Factorize the matrix
+  cholmodSolver.setMatrix(matrix);
+  cholmodSolver.factorize();
+
+  // Test with a particular right-hand side
+  Vector rhs = {3,4,4,3};
+
+  // The corresponding solution
+  Vector reference = {1,1,1,1};
+
+  cholmodSolver.setRhs(rhs);
+  cholmodSolver.solve();
+
+  // Check whether result is correct
+  auto diff = reference;
+  diff -= x;
+  if (diff.infinity_norm() > 1e-8)
+  {
+    std::cerr << "CholmodSolver didn't produce the correct result!" << std::endl;
+    passed = false;
+  }
 
+  // Test with a different right-hand side
+  // This reuses the original factorization.
+  rhs = {4,6,6,5};
+  reference = {1,2,1,2};
+  cholmodSolver.solve();
+
+  diff = reference;
+  diff -= x;
+  if (diff.infinity_norm() > 1e-8)
+  {
+    std::cerr << "CholmodSolver didn't produce the correct result!" << std::endl;
+    passed = false;
+  }
+
+  return passed;
+}
 
 int main(int argc, char** argv)
 {
@@ -72,5 +145,7 @@ int main(int argc, char** argv)
     passed &= checkWithStandardGrids(testsuite1);
     passed &= checkWithStandardGrids(testsuite2);
 
+    passed &= checkMultipleSolves();
+
     return passed ? 0 : 1;
 }