DUNE-FEM (unstable)

cholmod.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#pragma once
6
7#if HAVE_SUITESPARSE_CHOLMOD || defined DOXYGEN
8
12#include <dune/istl/bvector.hh>
13#include<dune/istl/solver.hh>
14#include <dune/istl/solverregistry.hh>
15#include <dune/istl/foreach.hh>
16
17#include <vector>
18#include <memory>
19
20#include <cholmod.h>
21
22namespace Dune {
23
24namespace Impl{
25
34 struct NoIgnore
35 {
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; }
39
40 };
41
42
43 template<class BlockedVector, class FlatVector>
44 void copyToFlatVector(const BlockedVector& blockedVector, FlatVector& flatVector)
45 {
46 // traverse the vector once just to compute the size
47 std::size_t len = flatVectorForEach(blockedVector, [&](auto&&, auto...){});
48 flatVector.resize(len);
49
50 flatVectorForEach(blockedVector, [&](auto&& entry, auto offset){
51 flatVector[offset] = entry;
52 });
53 }
54
55 // special (dummy) case for NoIgnore
56 template<class FlatVector>
57 void copyToFlatVector(const NoIgnore&, FlatVector&)
58 {
59 // just do nothing
60 return;
61 }
62
63 template<class FlatVector, class BlockedVector>
64 void copyToBlockedVector(const FlatVector& flatVector, BlockedVector& blockedVector)
65 {
66 flatVectorForEach(blockedVector, [&](auto& entry, auto offset){
67 entry = flatVector[offset];
68 });
69 }
70
71 // wrapper class for C function calls to CHOLMOD itself.
72 // The CHOLMOD API has different functions for different index types.
73 template <class Index>
74 struct CholmodMethodChooser;
75
76 // specialization using 'int' to store indices
77 template <>
78 struct CholmodMethodChooser<int>
79 {
80 [[nodiscard]]
81 static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
82 {
83 return ::cholmod_allocate_dense(nrow,ncol,d,xtype,c);
84 }
85
86 [[nodiscard]]
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)
88 {
89 return ::cholmod_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
90 }
91
92 [[nodiscard]]
93 static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
94 {
95 return ::cholmod_analyze(A,c);
96 }
97
98 static int defaults(cholmod_common *c)
99 {
100 return ::cholmod_defaults(c);
101 }
102
103 static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
104 {
105 return ::cholmod_factorize(A,L,c);
106 }
107
108 static int finish(cholmod_common *c)
109 {
110 return ::cholmod_finish(c);
111 }
112
113 static int free_dense (cholmod_dense **X, cholmod_common *c)
114 {
115 return ::cholmod_free_dense(X,c);
116 }
117
118 static int free_factor(cholmod_factor **L, cholmod_common *c)
119 {
120 return ::cholmod_free_factor(L,c);
121 }
122
123 static int free_sparse(cholmod_sparse **A, cholmod_common *c)
124 {
125 return ::cholmod_free_sparse(A,c);
126 }
127
128 [[nodiscard]]
129 static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
130 {
131 return ::cholmod_solve(sys,L,B,c);
132 }
133
134 static int start(cholmod_common *c)
135 {
136 return ::cholmod_start(c);
137 }
138 };
139
140 // specialization using 'SuiteSparse_long' to store indices
141 template <>
142 struct CholmodMethodChooser<SuiteSparse_long>
143 {
144 [[nodiscard]]
145 static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
146 {
147 return ::cholmod_l_allocate_dense(nrow,ncol,d,xtype,c);
148 }
149
150 [[nodiscard]]
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)
152 {
153 return ::cholmod_l_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
154 }
155
156 [[nodiscard]]
157 static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
158 {
159 return ::cholmod_l_analyze(A,c);
160 }
161
162 static int defaults(cholmod_common *c)
163 {
164 return ::cholmod_l_defaults(c);
165 }
166
167 static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
168 {
169 return ::cholmod_l_factorize(A,L,c);
170 }
171
172 static int finish(cholmod_common *c)
173 {
174 return ::cholmod_l_finish(c);
175 }
176
177 static int free_dense (cholmod_dense **X, cholmod_common *c)
178 {
179 return ::cholmod_l_free_dense(X,c);
180 }
181
182 static int free_factor (cholmod_factor **L, cholmod_common *c)
183 {
184 return ::cholmod_l_free_factor(L,c);
185 }
186
187 static int free_sparse(cholmod_sparse **A, cholmod_common *c)
188 {
189 return ::cholmod_l_free_sparse(A,c);
190 }
191
192 [[nodiscard]]
193 static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
194 {
195 return ::cholmod_l_solve(sys,L,B,c);
196 }
197
198 static int start(cholmod_common *c)
199 {
200 return ::cholmod_l_start(c);
201 }
202 };
203
204} //namespace Impl
205
214template<class Vector, class Index=int>
215class Cholmod : public InverseOperator<Vector, Vector>
216{
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'!");
219
220 using CholmodMethod = Impl::CholmodMethodChooser<Index>;
221
222public:
223
230 {
231 CholmodMethod::start(&c_);
232 }
233
240 {
241 if (L_)
242 CholmodMethod::free_factor(&L_, &c_);
243 CholmodMethod::finish(&c_);
244 }
245
246 // forbid copying to avoid freeing memory twice
247 Cholmod(const Cholmod&) = delete;
248 Cholmod& operator=(const Cholmod&) = delete;
249
250
253 void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
254 {
255 apply(x,b,res);
256 }
257
263 void apply(Vector& x, Vector& b, InverseOperatorResult& res) override
264 {
265 // do nothing if N=0
266 if ( nIsZero_ )
267 {
268 return;
269 }
270
271 if (x.size() != b.size())
272 DUNE_THROW(Exception, "Error in apply(): sizes of x and b do not match!");
273
274 // cast to double array
275 auto b2 = std::make_unique<double[]>(L_->n);
276 auto x2 = std::make_unique<double[]>(L_->n);
277
278 // copy to cholmod
279 auto bp = b2.get();
280
281 flatVectorForEach(b, [&](auto&& entry, auto&& flatIndex){
282 if ( subIndices_.empty() )
283 bp[ flatIndex ] = entry;
284 else
285 if( subIndices_[ flatIndex ] != std::numeric_limits<std::size_t>::max() )
286 bp[ subIndices_[ flatIndex ] ] = entry;
287 });
288
289 // create a cholmod dense object
290 auto b3 = make_cholmod_dense(CholmodMethod::allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
291
292 // cast because void-ptr
293 auto b4 = static_cast<double*>(b3->x);
294 std::copy(b2.get(), b2.get() + L_->n, b4);
295
296 // solve for a cholmod x object
297 auto x3 = make_cholmod_dense(CholmodMethod::solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
298 // cast because void-ptr
299 auto xp = static_cast<double*>(x3->x);
300
301 // copy into x
302 flatVectorForEach(x, [&](auto&& entry, auto&& flatIndex){
303 if ( subIndices_.empty() )
304 entry = xp[ flatIndex ];
305 else
306 if( subIndices_[ flatIndex ] != std::numeric_limits<std::size_t>::max() )
307 entry = xp[ subIndices_[ flatIndex ] ];
308 });
309
310 // statistics for a direct solver
311 res.iterations = 1;
312 res.converged = true;
313 }
314
315
321 template<class Matrix>
322 void setMatrix(const Matrix& matrix)
323 {
324 const Impl::NoIgnore* noIgnore = nullptr;
325 setMatrix(matrix, noIgnore);
326 }
327
342 template<class Matrix, class Ignore>
343 void setMatrix(const Matrix& matrix, const Ignore* ignore)
344 {
345 // count the number of entries and diagonal entries
346 size_t nonZeros = 0;
347 size_t numberOfIgnoredDofs = 0;
348
349
350 auto [flatRows,flatCols] = flatMatrixForEach( matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
351 if( flatRowIndex <= flatColIndex )
352 nonZeros++;
353 });
354
355 std::vector<bool> flatIgnore;
356
357 if ( ignore )
358 {
359 Impl::copyToFlatVector(*ignore,flatIgnore);
360 numberOfIgnoredDofs = std::count(flatIgnore.begin(),flatIgnore.end(),true);
361 }
362
363 nIsZero_ = (size_t(flatRows) <= numberOfIgnoredDofs);
364
365 if ( nIsZero_ )
366 {
367 return;
368 }
369
370 // Total number of rows
371 size_t N = flatRows - numberOfIgnoredDofs;
372
373 /*
374 * CHOLMOD uses compressed-column sparse matrices, but for symmetric
375 * matrices this is the same as the compressed-row sparse matrix used
376 * by DUNE. So we can just store Mᵀ instead of M (as M = Mᵀ).
377 */
378 const auto deleter = [c = &this->c_](auto* p) {
379 CholmodMethod::free_sparse(&p, c);
380 };
381 auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
382 CholmodMethod::allocate_sparse(N, // # rows
383 N, // # cols
384 nonZeros, // # of nonzeroes
385 1, // indices are sorted ( 1 = true)
386 1, // matrix is "packed" ( 1 = true)
387 -1, // stype of matrix ( -1 = consider the lower part only )
388 CHOLMOD_REAL, // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
389 &c_ // cholmod_common ptr
390 ), deleter);
391
392 // copy the data of BCRS matrix to Cholmod Sparse matrix
393 Index* Ap = static_cast<Index*>(M->p);
394 Index* Ai = static_cast<Index*>(M->i);
395 double* Ax = static_cast<double*>(M->x);
396
397
398 if ( ignore )
399 {
400 // init the mapping
401 subIndices_.resize(flatRows,std::numeric_limits<std::size_t>::max());
402
403 std::size_t subIndexCounter = 0;
404
405 for ( std::size_t i=0; i<flatRows; i++ )
406 {
407 if ( not flatIgnore[ i ] )
408 {
409 subIndices_[ i ] = subIndexCounter++;
410 }
411 }
412 }
413
414 // at first, we need to compute the row starts "Ap"
415 // therefore, we count all (not ignored) entries in each row and in the end we accumulate everything
416 flatMatrixForEach(matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
417
418 // stop if ignored
419 if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
420 return;
421
422 // stop if in lower half
423 if ( flatRowIndex > flatColIndex )
424 return;
425
426 // ok, count the entry
427 auto idx = ignore ? subIndices_[flatRowIndex] : flatRowIndex;
428 Ap[idx+1]++;
429
430 });
431
432 // now accumulate
433 Ap[0] = 0;
434 for ( size_t i=0; i<N; i++ )
435 {
436 Ap[i+1] += Ap[i];
437 }
438
439 // we need a compressed row position counter
440 std::vector<std::size_t> rowPosition(N,0);
441
442 // now we can set the entries
443 flatMatrixForEach(matrix, [&](auto&& entry, auto&& flatRowIndex, auto&& flatColIndex){
444
445 // stop if ignored
446 if ( ignore and ( flatIgnore[flatRowIndex] or flatIgnore[flatColIndex] ) )
447 return;
448
449 // stop if in lower half
450 if ( flatRowIndex > flatColIndex )
451 return;
452
453 // ok, set the entry
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]++;
461
462 });
463
464 // Now analyse the pattern and optimal row order
465 L_ = CholmodMethod::analyze(M.get(), &c_);
466
467 // Do the factorization (this may take some time)
468 CholmodMethod::factorize(M.get(), L_, &c_);
469 }
470
472 {
473 return SolverCategory::Category::sequential;
474 }
475
481 cholmod_common& cholmodCommonObject()
482 {
483 return c_;
484 }
485
491 cholmod_factor& cholmodFactor()
492 {
493 return *L_;
494 }
495
501 const cholmod_factor& cholmodFactor() const
502 {
503 return *L_;
504 }
505private:
506
507 // create a std::unique_ptr to a cholmod_dense object with a deleter
508 // that calls the appropriate cholmod cleanup routine
509 auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
510 {
511 const auto deleter = [c](auto* p) {
512 CholmodMethod::free_dense(&p, c);
513 };
514 return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
515 }
516
517 cholmod_common c_;
518 cholmod_factor* L_ = nullptr;
519
520 // indicator for a 0x0 problem (due to ignore dof's)
521 bool nIsZero_ = false;
522
523 // vector mapping all indices in flat order to the not ignored indices
524 std::vector<std::size_t> subIndices_;
525};
526
527 DUNE_REGISTER_SOLVER("cholmod",
528 [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
529 -> std::shared_ptr<typename decltype(opTraits)::solver_type>
530 {
531 using OpTraits = decltype(opTraits);
532 using M = typename OpTraits::matrix_type;
533 using D = typename OpTraits::domain_type;
534 // works only for sequential operators
535 if constexpr (OpTraits::isParallel){
536 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
537 DUNE_THROW(Dune::InvalidStateException, "CholMod works only for sequential operators.");
538 }
539 if constexpr (OpTraits::isAssembled &&
540 // check whether the Matrix field_type is double or float
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);
547 return solver;
548 }
549 DUNE_THROW(UnsupportedType,
550 "Unsupported Type in Cholmod (only double and float supported)");
551 });
552
553} /* namespace Dune */
554
555#endif // HAVE_SUITESPARSE_CHOLMOD
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)