7#if HAVE_SUITESPARSE_CHOLMOD || defined DOXYGEN
14#include <dune/istl/solverregistry.hh>
15#include <dune/istl/foreach.hh>
36 const NoIgnore& operator[](std::size_t)
const {
return *
this; }
37 explicit operator bool()
const {
return false; }
38 static constexpr std::size_t
size() {
return 0; }
43 template<
class BlockedVector,
class FlatVector>
44 void copyToFlatVector(
const BlockedVector& blockedVector, FlatVector& flatVector)
48 flatVector.resize(len);
51 flatVector[offset] = entry;
56 template<
class FlatVector>
57 void copyToFlatVector(
const NoIgnore&, FlatVector&)
63 template<
class FlatVector,
class BlockedVector>
64 void copyToBlockedVector(
const FlatVector& flatVector, BlockedVector& blockedVector)
67 entry = flatVector[offset];
73 template <
class Index>
74 struct CholmodMethodChooser;
78 struct CholmodMethodChooser<int>
81 static cholmod_dense* allocate_dense(
size_t nrow,
size_t ncol,
size_t d,
int xtype, cholmod_common *c)
83 return ::cholmod_allocate_dense(nrow,ncol,d,xtype,c);
87 static cholmod_sparse* allocate_sparse(
size_t nrow,
size_t ncol,
size_t nzmax,
int sorted,
int packed,
int stype,
int xtype, cholmod_common *c)
89 return ::cholmod_allocate_sparse(nrow,ncol,nzmax,
sorted,packed,stype,xtype,c);
93 static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
95 return ::cholmod_analyze(A,c);
98 static int defaults(cholmod_common *c)
100 return ::cholmod_defaults(c);
103 static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
105 return ::cholmod_factorize(A,L,c);
108 static int finish(cholmod_common *c)
110 return ::cholmod_finish(c);
113 static int free_dense (cholmod_dense **X, cholmod_common *c)
115 return ::cholmod_free_dense(X,c);
118 static int free_factor(cholmod_factor **L, cholmod_common *c)
120 return ::cholmod_free_factor(L,c);
123 static int free_sparse(cholmod_sparse **A, cholmod_common *c)
125 return ::cholmod_free_sparse(A,c);
129 static cholmod_dense* solve(
int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
131 return ::cholmod_solve(sys,L,B,c);
134 static int start(cholmod_common *c)
136 return ::cholmod_start(c);
142 struct CholmodMethodChooser<SuiteSparse_long>
145 static cholmod_dense* allocate_dense(
size_t nrow,
size_t ncol,
size_t d,
int xtype, cholmod_common *c)
147 return ::cholmod_l_allocate_dense(nrow,ncol,d,xtype,c);
151 static cholmod_sparse* allocate_sparse(
size_t nrow,
size_t ncol,
size_t nzmax,
int sorted,
int packed,
int stype,
int xtype, cholmod_common *c)
153 return ::cholmod_l_allocate_sparse(nrow,ncol,nzmax,
sorted,packed,stype,xtype,c);
157 static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
159 return ::cholmod_l_analyze(A,c);
162 static int defaults(cholmod_common *c)
164 return ::cholmod_l_defaults(c);
167 static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
169 return ::cholmod_l_factorize(A,L,c);
172 static int finish(cholmod_common *c)
174 return ::cholmod_l_finish(c);
177 static int free_dense (cholmod_dense **X, cholmod_common *c)
179 return ::cholmod_l_free_dense(X,c);
182 static int free_factor (cholmod_factor **L, cholmod_common *c)
184 return ::cholmod_l_free_factor(L,c);
187 static int free_sparse(cholmod_sparse **A, cholmod_common *c)
189 return ::cholmod_l_free_sparse(A,c);
193 static cholmod_dense* solve(
int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
195 return ::cholmod_l_solve(sys,L,B,c);
198 static int start(cholmod_common *c)
200 return ::cholmod_l_start(c);
214template<
class Vector,
class Index=
int>
217 static_assert(std::is_same_v<Index,int> || std::is_same_v<Index,SuiteSparse_long>,
218 "Index type must be either 'int' or 'SuiteSparse_long'!");
220 using CholmodMethod = Impl::CholmodMethodChooser<Index>;
231 CholmodMethod::start(&c_);
242 CholmodMethod::free_factor(&L_, &c_);
243 CholmodMethod::finish(&c_);
271 if (x.size() != b.size())
275 auto b2 = std::make_unique<double[]>(L_->n);
276 auto x2 = std::make_unique<double[]>(L_->n);
282 if ( subIndices_.empty() )
283 bp[ flatIndex ] = entry;
286 bp[ subIndices_[ flatIndex ] ] = entry;
290 auto b3 = make_cholmod_dense(CholmodMethod::allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
293 auto b4 =
static_cast<double*
>(b3->x);
294 std::copy(b2.get(), b2.get() + L_->n, b4);
297 auto x3 = make_cholmod_dense(CholmodMethod::solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
299 auto xp =
static_cast<double*
>(x3->x);
303 if ( subIndices_.empty() )
304 entry = xp[ flatIndex ];
307 entry = xp[ subIndices_[ flatIndex ] ];
321 template<
class Matrix>
324 const Impl::NoIgnore* noIgnore =
nullptr;
342 template<
class Matrix,
class Ignore>
347 size_t numberOfIgnoredDofs = 0;
350 auto [flatRows,flatCols] =
flatMatrixForEach( matrix, [&](
auto&& ,
auto&& flatRowIndex,
auto&& flatColIndex){
351 if( flatRowIndex <= flatColIndex )
355 std::vector<bool> flatIgnore;
359 Impl::copyToFlatVector(*ignore,flatIgnore);
360 numberOfIgnoredDofs = std::count(flatIgnore.begin(),flatIgnore.end(),
true);
363 nIsZero_ = (size_t(flatRows) <= numberOfIgnoredDofs);
371 size_t N = flatRows - numberOfIgnoredDofs;
378 const auto deleter = [c = &this->c_](
auto* p) {
379 CholmodMethod::free_sparse(&p, c);
381 auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
382 CholmodMethod::allocate_sparse(N,
393 Index* Ap =
static_cast<Index*
>(M->p);
394 Index* Ai =
static_cast<Index*
>(M->i);
395 double* Ax =
static_cast<double*
>(M->x);
403 std::size_t subIndexCounter = 0;
405 for ( std::size_t i=0; i<flatRows; i++ )
407 if ( not flatIgnore[ i ] )
409 subIndices_[ i ] = subIndexCounter++;
416 flatMatrixForEach(matrix, [&](
auto&& ,
auto&& flatRowIndex,
auto&& flatColIndex){
419 if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
423 if ( flatRowIndex > flatColIndex )
427 auto idx = ignore ? subIndices_[flatRowIndex] : flatRowIndex;
434 for (
size_t i=0; i<N; i++ )
440 std::vector<std::size_t> rowPosition(N,0);
443 flatMatrixForEach(matrix, [&](
auto&& entry,
auto&& flatRowIndex,
auto&& flatColIndex){
446 if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
450 if ( flatRowIndex > flatColIndex )
454 auto rowIdx = ignore ? subIndices_[flatRowIndex] : flatRowIndex;
455 auto colIdx = ignore ? subIndices_[flatColIndex] : flatColIndex;
456 auto rowStart = Ap[rowIdx];
457 auto rowPos = rowPosition[rowIdx];
458 Ai[ rowStart + rowPos ] = colIdx;
459 Ax[ rowStart + rowPos ] = entry;
460 rowPosition[rowIdx]++;
465 L_ = CholmodMethod::analyze(M.get(), &c_);
468 CholmodMethod::factorize(M.get(), L_, &c_);
473 return SolverCategory::Category::sequential;
509 auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
511 const auto deleter = [c](
auto* p) {
512 CholmodMethod::free_dense(&p, c);
514 return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
518 cholmod_factor* L_ =
nullptr;
521 bool nIsZero_ =
false;
524 std::vector<std::size_t> subIndices_;
527 DUNE_REGISTER_SOLVER(
"cholmod",
529 -> std::shared_ptr<
typename decltype(opTraits)::solver_type>
531 using OpTraits =
decltype(opTraits);
532 using M =
typename OpTraits::matrix_type;
533 using D =
typename OpTraits::domain_type;
535 if constexpr (OpTraits::isParallel){
536 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
539 if constexpr (OpTraits::isAssembled &&
541 (std::is_same_v<typename FieldTraits<D>::field_type,
double> ||
542 std::is_same_v<typename FieldTraits<D>::field_type,
float>)){
543 const auto& A = opTraits.getAssembledOpOrThrow(op);
544 const M& mat = A->getmat();
545 auto solver = std::make_shared<Dune::Cholmod<D>>();
546 solver->setMatrix(mat);
550 "Unsupported Type in Cholmod (only double and float supported)");
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Dune wrapper for SuiteSparse/CHOLMOD solver.
Definition: cholmod.hh:216
cholmod_common & cholmodCommonObject()
return a reference to the CHOLMOD common object for advanced option settings
Definition: cholmod.hh:481
void setMatrix(const Matrix &matrix)
Set matrix without ignore nodes.
Definition: cholmod.hh:322
cholmod_factor & cholmodFactor()
The CHOLMOD data structure that stores the factorization.
Definition: cholmod.hh:491
void apply(Vector &x, Vector &b, InverseOperatorResult &res) override
solve the linear system Ax=b (possibly with respect to some ignore field)
Definition: cholmod.hh:263
~Cholmod()
Destructor.
Definition: cholmod.hh:239
const cholmod_factor & cholmodFactor() const
The CHOLMOD data structure that stores the factorization.
Definition: cholmod.hh:501
Cholmod()
Default constructor.
Definition: cholmod.hh:229
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category)
Definition: cholmod.hh:471
void apply(Vector &x, Vector &b, double reduction, InverseOperatorResult &res) override
simple forward to apply(X&, Y&, InverseOperatorResult&)
Definition: cholmod.hh:253
void setMatrix(const Matrix &matrix, const Ignore *ignore)
Set matrix and ignore nodes.
Definition: cholmod.hh:343
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Abstract base class for all solvers.
Definition: solver.hh:101
A generic dynamic dense matrix.
Definition: matrix.hh:561
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune namespace.
Definition: alignedallocator.hh:13
constexpr auto sorted(std::integer_sequence< T, II... > seq, Compare comp)
Sort a given sequence by the comparator comp.
Definition: integersequence.hh:119
std::pair< std::size_t, std::size_t > flatMatrixForEach(Matrix &&matrix, F &&f, std::size_t rowOffset=0, std::size_t colOffset=0)
Traverse a blocked matrix and call a functor at each scalar entry.
Definition: foreach.hh:132
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
std::size_t flatVectorForEach(Vector &&vector, F &&f, std::size_t offset=0)
Traverse a blocked vector and call a functor at each scalar entry.
Definition: foreach.hh:95
Define general, extensible interface for inverse operators.
Statistics about the application of an inverse operator.
Definition: solver.hh:50
int iterations
Number of iterations.
Definition: solver.hh:69
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Category
Definition: solvercategory.hh:23