diff --git a/dune/solvers/solvers/cholmodsolver.hh b/dune/solvers/solvers/cholmodsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..50114c64e15cc2daa39e1540c5a2f9211822c9cf
--- /dev/null
+++ b/dune/solvers/solvers/cholmodsolver.hh
@@ -0,0 +1,165 @@
+// -*- 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
diff --git a/dune/solvers/test/CMakeLists.txt b/dune/solvers/test/CMakeLists.txt
index 30d1946e9479d5f34d42864aec42d70f63f1e945..314ad6ed239fa0103fca5211e0d5fc39fa4f70de 100644
--- a/dune/solvers/test/CMakeLists.txt
+++ b/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()
diff --git a/dune/solvers/test/cholmodsolvertest.cc b/dune/solvers/test/cholmodsolvertest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3411fa91ea6475724ba9bce4d38a88935f6e4d43
--- /dev/null
+++ b/dune/solvers/test/cholmodsolvertest.cc
@@ -0,0 +1,76 @@
+// -*- 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;
+}