Implementation of CHOLMOD Solver

This class wrappes the dune-istl CHOLMOD method for dune-solvers module.

It is basically a copy of UMFPack solver.

There are problems with some METIS implementations (especially LibScotchMETIS)
causing random segfaults. Therefore the usage of METIS is disabled for matrices larger
than 3000 rows for now (whichs seems stable),
see https://gitlab.inria.fr/scotch/scotch/issues/8

dune/solvers/solvers/cholmodsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_CHOLMODSOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_CHOLMODSOLVER_HH
+
+/** \file
+    \brief A wrapper for the CHOLMOD sparse direct solver
+*/
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/solver.hh>
+#include <dune/istl/cholmod.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
+
+namespace Dune
+{
+
+namespace Solvers
+{
+
+/** \brief Wraps the CHOLMOD sparse symmetric direct solver
+*
+*  CHOLMOD uses the Cholesky decomposition A = LL^T for sparse
+*  symmetric matrices to solve linear systems Ax = b efficiently.
+*
+*  This class wraps the dune-istl CHOLMOD solver for dune-solvers to provide
+*  a Solvers::LinearSolver interface.
+*/
+template <class MatrixType, class VectorType>
+class CholmodSolver
+: public LinearSolver<MatrixType,VectorType>, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+  enum {blocksize = VectorType::block_type::dimension};
+
+  using field_type = typename VectorType::field_type;
+
+  using  MatrixBlock = Dune::FieldMatrix<field_type, blocksize, blocksize> ;
+  using  VectorBlock = Dune::FieldVector<field_type, blocksize>;
+
+public:
+
+  /** \brief Default constructor */
+  CholmodSolver ()
+  : LinearSolver<MatrixType,VectorType>(NumProc::FULL)
+  {}
+
+  /** \brief Constructor for a linear problem */
+  CholmodSolver (const MatrixType& matrix,
+                 VectorType& x,
+                 const VectorType& rhs,
+                 NumProc::VerbosityMode verbosity=NumProc::FULL)
+  : LinearSolver<MatrixType,VectorType>(verbosity),
+    matrix_(&matrix),
+    x_(&x),
+    rhs_(&rhs)
+  {}
+
+  /** \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.
+   */
+  void setProblem(const MatrixType& matrix,
+                  VectorType& x,
+                  const VectorType& rhs) override
+  {
+    matrix_ = &matrix;
+    x_ = &x;
+    rhs_ = &rhs;
+  }
+
+  virtual void solve()
+  {
+    // 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;
+
+    if (not this->hasIgnore())
+    {
+      // the apply() method doesn't like constant values
+      auto mutableRhs = *rhs_;
+      solver.setMatrix(*matrix_);
+      solver.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);
+
+      auto modifiedRhs = *rhs_;
+
+      // pull all ignored columns to the rhs
+      for(auto rowIt = matrix_->begin(); rowIt != matrix_->end(); rowIt++)
+      {
+        const auto row = rowIt.index();
+
+        for(auto ii=0; ii<blocksize; ii++)
+        {
+          if( ignore[row][ii] )
+            continue;
+
+          for(auto colIt = rowIt->begin(); colIt != rowIt->end(); colIt++)
+          {
+            const auto col = colIt.index();
+            for(auto jj=0; jj<blocksize; jj++)
+            {
+              if( ignore[col][jj] )
+              {
+                // we assume values only in the upper part
+                if ( row < col )
+                {
+                  // upper part -> use canonical indices
+                  modifiedRhs[row][ii] -= (*matrix_)[row][col][ii][jj] * (*x_)[col][jj];
+                }
+                else if (row == col and ii > jj )
+                {
+                  // diagonal block but there in the lower part -> swap ii <-> jj
+                  modifiedRhs[row][ii] -= (*matrix_)[row][col][jj][ii] * (*x_)[col][jj];
+                }
+                else
+                {
+                  // lower part -> swap col <-> row and ii <-> jj
+                  modifiedRhs[row][ii] -= (*matrix_)[col][row][jj][ii] * (*x_)[col][jj];
+                }
+              }
+            }
+          }
+        }
+      }
+
+      // Solve the modified system
+      solver.apply(*x_, modifiedRhs, statistics);
+    }
+  }
+
+  //! Pointer to the system matrix
+  const MatrixType* matrix_;
+
+  //! Pointer to the solution vector
+  VectorType* x_;
+
+  //! Pointer to the right hand side
+  const VectorType* rhs_;
+};
+
+}   // namespace Solvers
+}   // namespace Dune
+
+#endif

dune/solvers/test/CMakeLists.txt

@@ -22,3 +22,8 @@ if(SuiteSparse_UMFPACK_FOUND)
   dune_add_test(SOURCES umfpacksolvertest.cc)
   add_dune_suitesparse_flags(umfpacksolvertest)
 endif()
+
+if(SuiteSparse_CHOLMOD_FOUND)
+  dune_add_test(SOURCES cholmodsolvertest.cc)
+  add_dune_suitesparse_flags(cholmodsolvertest)
+endif()

dune/solvers/test/cholmodsolvertest.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
+#include <config.h>
+
+#include <cmath>
+#include <iostream>
+#include <sstream>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/cholmodsolver.hh>
+
+
+#include "common.hh"
+
+using namespace Dune;
+
+/**  \brief test for the CholmodSolver class
+  *
+  *  This test tests whether the Dirichlet problem for a Laplace operator is solved correctly for a "random" rhs by an CholmodSolver.
+  *  Implementation is more or less stolen from UMFPackSolverTest
+  */
+template <size_t blocksize>
+struct CHOLMODSolverTestSuite
+{
+  template <class GridType>
+  bool check(const GridType& grid)
+  {
+    double tol = 1e-12;
+
+    bool passed = true;
+
+    using Problem =
+        SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
+    Problem p(grid.leafGridView());
+
+    typedef Solvers::CholmodSolver<typename Problem::Matrix,
+                                   typename Problem::Vector> Solver;
+    Solver solver(p.A,p.u,p.rhs);
+    solver.setIgnore(p.ignore);
+
+    // solve problem
+    solver.preprocess();
+    solver.solve();
+
+    if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
+    {
+      std::cout << "The CholmodSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
+      std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
+      passed = false;
+    }
+
+    return passed;
+  }
+};
+
+
+
+int main(int argc, char** argv)
+{
+    Dune::MPIHelper::instance(argc, argv);
+
+    bool passed = true;
+
+    CHOLMODSolverTestSuite<1> testsuite1;
+    CHOLMODSolverTestSuite<2> testsuite2;
+    passed = checkWithStandardGrids(testsuite1);
+    passed = checkWithStandardGrids(testsuite2);
+
+    return passed ? 0 : 1;
+}