8#include <dune/istl/solverfactory.hh>
29 const NoIgnore& operator[](std::size_t)
const {
return *
this; }
30 explicit operator bool()
const {
return false; }
31 std::size_t count()
const {
return 0; }
43template<
class T,
class A,
int k>
45 :
public InverseOperator<BlockVector<FieldVector<T,k>, A>, BlockVector<FieldVector<T,k>, A>>
73 cholmod_free_factor(&L_, &c_);
78 Cholmod(
const Cholmod&) =
delete;
79 Cholmod& operator=(
const Cholmod&) =
delete;
103 if (x.size() != b.size())
106 const auto& blocksize = k;
109 auto b2 = std::make_unique<double[]>(L_->n);
110 auto x2 = std::make_unique<double[]>(L_->n);
114 if (inverseSubIndices_.empty())
117 for (
const auto& bi : b)
118 for (
const auto& bii : bi)
124 for (
const auto& idx : inverseSubIndices_)
125 *bp++ = b[ idx / blocksize ][ idx % blocksize ];
129 auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
131 auto b4 =
static_cast<double*
>(b3->x);
132 std::copy(b2.get(), b2.get() + L_->n, b4);
135 auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
137 auto xp =
static_cast<double*
>(x3->x);
140 if (inverseSubIndices_.empty())
143 for (
int i=0, s=x.size(); i<s; i++)
144 for (
int ii=0, ss=x[i].size(); ii<ss; ii++)
150 for (
const auto& idx : inverseSubIndices_)
151 x[ idx / blocksize ][ idx % blocksize ] = *xp++;
167 const Impl::NoIgnore* noIgnore =
nullptr;
168 setMatrix(matrix, noIgnore);
185 template<
class Ignore>
189 const auto blocksize = k;
192 int N = blocksize * matrix.N();
194 N -= ignore->count();
204 const int nnz = blocksize * blocksize * matrix.nonzeroes();
206 const int nDiag = blocksize * blocksize * matrix.N();
208 const int nnzDiag = (blocksize * (blocksize+1)) / 2 * matrix.N();
210 const int nnzTri = (nnz - nDiag) / 2 + nnzDiag;
217 const auto deleter = [c = &this->c_](
auto* p) {
218 cholmod_free_sparse(&p, c);
220 auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
221 cholmod_allocate_sparse(N,
232 int* Ap =
static_cast<int*
>(M->p);
233 int* Ai =
static_cast<int*
>(M->i);
234 double* Ax =
static_cast<double*
>(M->x);
237 std::vector<std::size_t> subIndices;
242 inverseSubIndices_.resize(N);
243 subIndices.resize(matrix.M()*blocksize);
246 for( std::size_t block=0; block<matrix.N(); block++ )
248 for( std::size_t i=0; i<blocksize; i++ )
250 if( not (*ignore)[block][i] )
252 subIndices[ block*blocksize + i ] = j;
253 inverseSubIndices_[j] = block*blocksize + i;
262 for (
auto rowIt = matrix.begin(); rowIt != matrix.end(); rowIt++)
264 const auto row = rowIt.index();
265 for (std::size_t i=0; i<blocksize; i++)
267 if( ignore and (*ignore)[row][i] )
273 for (
auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
275 const auto col = colIt.index();
281 for (
auto j = (row == col ? i : 0); j<blocksize; j++)
283 if( ignore and (*ignore)[col][j] )
286 const auto jj = ignore ? subIndices[ blocksize*col + j ] : blocksize*col + j;
290 *Ax++ = (*colIt)[i][j];
301 L_ = cholmod_analyze(M.get(), &c_);
304 cholmod_factorize(M.get(), L_, &c_);
309 return SolverCategory::Category::sequential;
315 auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
317 const auto deleter = [c](
auto* p) {
318 cholmod_free_dense(&p, c);
320 return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
324 cholmod_factor* L_ =
nullptr;
327 bool nIsZero_ =
false;
331 std::vector<std::size_t> inverseSubIndices_;
335 struct CholmodCreator{
336 template<
class F>
struct isValidBlock : std::false_type{};
337 template<
int k>
struct isValidBlock<FieldVector<double,k>> : std::true_type{};
338 template<
int k>
struct isValidBlock<FieldVector<float,k>> : std::true_type{};
340 template<
class TL,
typename M>
341 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
342 typename Dune::TypeListElement<2, TL>::type>>
344 std::enable_if_t<isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
346 using D =
typename Dune::TypeListElement<1, TL>::type;
347 auto solver = std::make_shared<Dune::Cholmod<D>>();
348 solver->setMatrix(mat);
353 template<
typename TL,
typename M>
354 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
355 typename Dune::TypeListElement<2, TL>::type>>
357 std::enable_if_t<!isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
359 DUNE_THROW(UnsupportedType,
"Unsupported Type in Cholmod");
362 DUNE_REGISTER_DIRECT_SOLVER(
"cholmod", Dune::CholmodCreator());
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:426
A vector of blocks with memory management.
Definition: bvector.hh:403
void setMatrix(const BCRSMatrix< FieldMatrix< T, k, k > > &matrix, const Ignore *ignore)
Set matrix and ignore nodes.
Definition: cholmod.hh:186
void apply(X &x, B &b, double reduction, InverseOperatorResult &res)
simple forward to apply(X&, Y&, InverseOperatorResult&)
Definition: cholmod.hh:84
void setMatrix(const BCRSMatrix< FieldMatrix< T, k, k > > &matrix)
Set matrix without ignore nodes.
Definition: cholmod.hh:165
~Cholmod()
Destructor.
Definition: cholmod.hh:70
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: cholmod.hh:307
Cholmod()
Default constructor.
Definition: cholmod.hh:60
void apply(X &x, B &b, InverseOperatorResult &res)
solve the linear system Ax=b (possibly with respect to some ignore field)
Definition: cholmod.hh:95
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
A dense n x m matrix.
Definition: fmatrix.hh:69
vector space out of a tensor product of fields.
Definition: fvector.hh:96
Abstract base class for all solvers.
Definition: solver.hh:97
Hierarchical structure of string parameters.
Definition: parametertree.hh:35
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.
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:46
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:14
Define general, extensible interface for inverse operators.
Statistics about the application of an inverse operator.
Definition: solver.hh:46
int iterations
Number of iterations.
Definition: solver.hh:65
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Category
Definition: solvercategory.hh:21