Merge branch 'feature/cholmod-reuse-factorization' into 'master'

CholmodSolver: Allow to reuse a factorization

See merge request agnumpde/dune-solvers!67

CHANGELOG.md

@@ -10,6 +10,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)
 {
@@ -69,8 +142,10 @@ int main(int argc, char** argv)
 
     CHOLMODSolverTestSuite<1> testsuite1;
     CHOLMODSolverTestSuite<2> testsuite2;
-    passed = checkWithStandardGrids(testsuite1);
-    passed = checkWithStandardGrids(testsuite2);
+    passed &= checkWithStandardGrids(testsuite1);
+    passed &= checkWithStandardGrids(testsuite2);
+
+    passed &= checkMultipleSolves();
 
     return passed ? 0 : 1;
 }