3#ifndef DUNE_ISTL_FASTAMG_HH
4#define DUNE_ISTL_FASTAMG_HH
20#include "fastamgsmoother.hh"
56 template<
class M,
class X,
class PI=SequentialInformation,
class A=std::allocator<X> >
136 std::size_t levels();
138 std::size_t maxlevels();
167 void createHierarchies(C& criterion,
168 const std::shared_ptr<const Operator>& matrixptr,
190 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
194 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
214 void mgc(LevelContext& levelContext,
Domain& x,
const Range& b);
222 void presmooth(LevelContext& levelContext,
Domain& x,
const Range& b);
230 void postsmooth(LevelContext& levelContext,
Domain& x,
const Range& b);
238 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel,
245 bool moveToCoarseLevel(LevelContext& levelContext);
251 void initIteratorsWithFineLevel(LevelContext& levelContext);
254 std::shared_ptr<OperatorHierarchy> matrices_;
256 std::shared_ptr<CoarseSolver> solver_;
258 std::shared_ptr<Hierarchy<Range,A>> rhs_;
260 std::shared_ptr<Hierarchy<Domain,A>> lhs_;
262 std::shared_ptr<Hierarchy<Domain,A>> residual_;
267 std::shared_ptr<ScalarProduct> scalarProduct_;
271 std::size_t preSteps_;
273 std::size_t postSteps_;
275 bool buildHierarchy_;
277 bool coarsesolverconverged;
279 typedef std::shared_ptr<Smoother> SmootherPointer;
280 SmootherPointer coarseSmoother_;
282 std::size_t verbosity_;
285 template<
class M,
class X,
class PI,
class A>
287 : matrices_(amg.matrices_), solver_(amg.solver_),
288 rhs_(), lhs_(), residual_(), scalarProduct_(amg.scalarProduct_),
289 gamma_(amg.gamma_), preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
290 symmetric(amg.symmetric), coarsesolverconverged(amg.coarsesolverconverged),
291 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
294 template<
class M,
class X,
class PI,
class A>
298 rhs_(), lhs_(), residual_(), scalarProduct_(),
299 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
300 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
301 symmetric(symmetric_), coarsesolverconverged(true),
302 coarseSmoother_(), verbosity_(parms.debugLevel())
304 if(preSteps_>1||postSteps_>1)
306 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
307 preSteps_=postSteps_=0;
309 assert(matrices_->isBuilt());
310 static_assert(std::is_same<PI,SequentialInformation>::value,
311 "Currently only sequential runs are supported");
313 template<
class M,
class X,
class PI,
class A>
320 : solver_(), rhs_(), lhs_(), residual_(), scalarProduct_(), gamma_(parms.getGamma()),
321 preSteps_(parms.getNoPreSmoothSteps()), postSteps_(parms.getNoPostSmoothSteps()),
322 buildHierarchy_(true),
323 symmetric(symmetric_), coarsesolverconverged(true),
324 coarseSmoother_(), verbosity_(criterion.debugLevel())
326 if(preSteps_>1||postSteps_>1)
328 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
329 preSteps_=postSteps_=1;
331 static_assert(std::is_same<PI,SequentialInformation>::value,
332 "Currently only sequential runs are supported");
337 createHierarchies(criterion, matrixptr, pinfo);
340 template<
class M,
class X,
class PI,
class A>
343 const std::shared_ptr<const Operator>& matrixptr,
347 matrices_ = std::make_shared<OperatorHierarchy>(
348 std::const_pointer_cast<Operator>(matrixptr),
351 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
353 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
354 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.
elapsed()<<
" seconds."<<std::endl;
356 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
360 sargs.iterations = 1;
363 cargs.setArgs(sargs);
364 if(matrices_->redistributeInformation().back().isSetup()) {
366 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
367 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
369 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
370 cargs.setComm(*matrices_->parallelInformation().coarsest());
374 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
376#if HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK
377#if HAVE_SUITESPARSE_UMFPACK
378#define DIRECTSOLVER UMFPack
380#define DIRECTSOLVER SuperLU
383 if(std::is_same<ParallelInformation,SequentialInformation>::value
384 || matrices_->parallelInformation().coarsest()->communicator().size()==1
385 || (matrices_->parallelInformation().coarsest().isRedistributed()
386 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
387 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
388 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
389 std::cout<<
"Using superlu"<<std::endl;
390 if(matrices_->parallelInformation().coarsest().isRedistributed())
392 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
394 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
398 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
403 if(matrices_->parallelInformation().coarsest().isRedistributed())
405 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
407 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
409 *coarseSmoother_, 1E-2, 1000, 0));
413 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
415 *coarseSmoother_, 1E-2, 1000, 0));
419 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
420 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.
elapsed()<<
" seconds."<<std::endl;
424 template<
class M,
class X,
class PI,
class A>
432 typedef typename M::matrix_type
Matrix;
439 const Matrix& mat=matrices_->matrices().finest()->getmat();
440 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
441 bool isDirichlet =
true;
442 bool hasDiagonal =
false;
444 for(ColIter col=row->begin(); col!=row->end(); ++col) {
445 if(row.index()==col.index()) {
447 hasDiagonal = (*col != zero);
453 if(isDirichlet && hasDiagonal)
454 diag->solve(x[row.index()], b[row.index()]);
457 std::cout<<
" Preprocessing Dirichlet took "<<watch1.
elapsed()<<std::endl;
460 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
461 rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
462 lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
463 residual_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
464 matrices_->coarsenVector(*rhs_);
465 matrices_->coarsenVector(*lhs_);
466 matrices_->coarsenVector(*residual_);
473 template<
class M,
class X,
class PI,
class A>
476 return matrices_->levels();
478 template<
class M,
class X,
class PI,
class A>
479 std::size_t FastAMG<M,X,PI,A>::maxlevels()
481 return matrices_->maxlevels();
485 template<
class M,
class X,
class PI,
class A>
488 LevelContext levelContext;
490 initIteratorsWithFineLevel(levelContext);
492 assert(v.two_norm()==0);
495 if(matrices_->maxlevels()==1){
498 mgc(levelContext, v, b);
500 mgc(levelContext, v, d);
501 if(postSteps_==0||matrices_->maxlevels()==1)
502 levelContext.pinfo->copyOwnerToAll(v, v);
505 template<
class M,
class X,
class PI,
class A>
508 levelContext.matrix = matrices_->matrices().finest();
509 levelContext.pinfo = matrices_->parallelInformation().finest();
510 levelContext.redist =
511 matrices_->redistributeInformation().begin();
512 levelContext.aggregates = matrices_->aggregatesMaps().begin();
513 levelContext.lhs = lhs_->finest();
514 levelContext.residual = residual_->finest();
515 levelContext.rhs = rhs_->finest();
516 levelContext.level=0;
519 template<
class M,
class X,
class PI,
class A>
520 bool FastAMG<M,X,PI,A>
521 ::moveToCoarseLevel(LevelContext& levelContext)
523 bool processNextLevel=
true;
525 if(levelContext.redist->isSetup()) {
527 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.residual),
528 levelContext.residual.getRedistributed());
529 processNextLevel = levelContext.residual.getRedistributed().size()>0;
530 if(processNextLevel) {
532 ++levelContext.pinfo;
533 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
534 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
535 static_cast<const Range&
>(levelContext.residual.getRedistributed()),
536 *levelContext.pinfo);
541 ++levelContext.pinfo;
542 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
543 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
544 static_cast<const Range&
>(*levelContext.residual), *levelContext.pinfo);
547 if(processNextLevel) {
549 ++levelContext.residual;
551 ++levelContext.matrix;
552 ++levelContext.level;
553 ++levelContext.redist;
555 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
557 ++levelContext.aggregates;
561 *levelContext.residual=0;
563 return processNextLevel;
566 template<
class M,
class X,
class PI,
class A>
567 void FastAMG<M,X,PI,A>
568 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel, Domain& x)
570 if(processNextLevel) {
571 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
573 --levelContext.aggregates;
575 --levelContext.redist;
576 --levelContext.level;
578 --levelContext.matrix;
579 --levelContext.residual;
584 if(levelContext.redist->isSetup()) {
587 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
588 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
589 levelContext.lhs.getRedistributed(),
590 matrices_->getProlongationDampingFactor(),
591 *levelContext.pinfo, *levelContext.redist);
593 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
594 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
595 matrices_->getProlongationDampingFactor(), *levelContext.pinfo);
601 if(processNextLevel) {
608 template<
class M,
class X,
class PI,
class A>
609 void FastAMG<M,X,PI,A>
610 ::presmooth(LevelContext& levelContext, Domain& x,
const Range& b)
612 constexpr auto bl = blockLevel<typename M::matrix_type>();
613 GaussSeidelPresmoothDefect<bl>::apply(levelContext.matrix->getmat(),
615 *levelContext.residual,
619 template<
class M,
class X,
class PI,
class A>
620 void FastAMG<M,X,PI,A>
621 ::postsmooth(LevelContext& levelContext, Domain& x,
const Range& b)
623 constexpr auto bl = blockLevel<typename M::matrix_type>();
624 GaussSeidelPostsmoothDefect<bl>
625 ::apply(levelContext.matrix->getmat(), x, *levelContext.residual, b);
629 template<
class M,
class X,
class PI,
class A>
632 return IsDirectSolver< CoarseSolver>::value;
635 template<
class M,
class X,
class PI,
class A>
638 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
642 if(levelContext.redist->isSetup()) {
643 levelContext.redist->redistribute(b, levelContext.rhs.getRedistributed());
644 if(levelContext.rhs.getRedistributed().size()>0) {
646 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
647 levelContext.rhs.getRedistributed());
648 solver_->apply(levelContext.lhs.getRedistributed(), levelContext.rhs.getRedistributed(), res);
650 levelContext.redist->redistributeBackward(v, levelContext.lhs.getRedistributed());
651 levelContext.pinfo->copyOwnerToAll(v, v);
653 levelContext.pinfo->copyOwnerToAll(b, b);
654 solver_->apply(v,
const_cast<Range&
>(b), res);
660 coarsesolverconverged =
false;
666#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
667 bool processNextLevel = moveToCoarseLevel(levelContext);
669 if(processNextLevel) {
671 for(std::size_t i=0; i<gamma_; i++)
672 mgc(levelContext, *levelContext.lhs, *levelContext.rhs);
675 moveToFineLevel(levelContext, processNextLevel, v);
680 if(levelContext.matrix == matrices_->matrices().finest()) {
681 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
682 if(!coarsesolverconverged)
683 DUNE_THROW(MathError,
"Coarse solver did not converge");
692 template<
class M,
class X,
class PI,
class A>
700 template<
class M,
class X,
class PI,
class A>
704 matrices_->getCoarsestAggregatesOnFinest(cont);
A fast (sequential) algebraic multigrid based on agglomeration that saves memory bandwidth.
Definition: fastamg.hh:58
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:214
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:217
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:59
All parameters for AMG.
Definition: parameters.hh:391
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
The negation of a set. An item is contained in the set if and only if it is not contained in the nega...
Definition: enumset.hh:94
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:50
Sequential SSOR preconditioner.
Definition: preconditioners.hh:139
A simple stop watch.
Definition: timer.hh:41
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:55
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:75
A few common exception classes.
Traits for type conversions and type information.
Some generic functions for pretty printing vectors and matrices.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: fastamg.hh:182
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: fastamg.hh:702
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: fastamg.hh:186
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: fastamg.hh:148
void post(Domain &x)
Clean up.
Definition: fastamg.hh:693
X Domain
The domain type.
Definition: fastamg.hh:75
Hierarchy< Domain, A >::Iterator residual
The iterator over the residuals.
Definition: fastamg.hh:202
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: fastamg.hh:121
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: fastamg.hh:70
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: fastamg.hh:190
X Range
The range type.
Definition: fastamg.hh:77
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: fastamg.hh:68
M Operator
The matrix operator type.
Definition: fastamg.hh:61
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: fastamg.hh:79
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: fastamg.hh:630
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: fastamg.hh:198
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: fastamg.hh:206
std::size_t level
The level index.
Definition: fastamg.hh:210
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: fastamg.hh:486
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: fastamg.hh:72
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: fastamg.hh:425
FastAMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const Parameters &parms, bool symmetric=true)
Construct a new amg with a specific coarse solver.
Definition: fastamg.hh:295
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: fastamg.hh:194
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:42
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:50
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
Provides a classes representing the hierarchies in AMG.
Dune namespace.
Definition: alignedallocator.hh:11
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:70
Define general preconditioner interface.
Define base class for scalar product and norm.
Classes for the generic construction and application of the smoothers.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
Statistics about the application of an inverse operator.
Definition: solver.hh:46
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Classes for using SuperLU with ISTL matrices.
Prolongation and restriction for amg.
Classes for using UMFPack with ISTL matrices.