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