Merge branch 'feature/cholmod-solver' into 'master'

Implementation of CHOLMOD Solver

See merge request agnumpde/dune-solvers!36

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;
+}