3 #include <dune/common/fmatrix.hh>
4 #include <dune/common/fvector.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; }
43 template<
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_);
86 DUNE_UNUSED_PARAMETER(reduction);
103 if (x.size() != b.size())
104 DUNE_THROW(Exception,
"Error in apply(): sizes of x and b do not match!");
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_;
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
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Define general, extensible interface for inverse operators.
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:7
DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator())
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:426
A vector of blocks with memory management.
Definition: bvector.hh:403
Definition: cholmod.hh:37
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(const Cholmod &)=delete
~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 setMatrix(const BCRSMatrix< FieldMatrix< T, k, k >> &matrix, const Ignore *ignore)
Set matrix and ignore nodes.
Definition: cholmod.hh:186
Cholmod & operator=(const Cholmod &)=delete
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
Definition: cholmod.hh:335
std::shared_ptr< Dune::InverseOperator< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL, const M &mat, const Dune::ParameterTree &, std::enable_if_t< isValidBlock< typename Dune::TypeListElement< 1, TL >::type::block_type >::value, int >=0) const
Definition: cholmod.hh:343
Definition: cholmod.hh:336
Definition: matrixutils.hh:26
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
Abstract base class for all solvers.
Definition: solver.hh:97
Category
Definition: solvercategory.hh:21
Definition: solverfactory.hh:124