5#ifndef DUNE_ISTL_FASTAMG_HH
6#define DUNE_ISTL_FASTAMG_HH
22#include "fastamgsmoother.hh"
58 template<
class M,
class X,
class PI=SequentialInformation,
class A=std::allocator<X> >
106 FastAMG(std::shared_ptr<const Operator> fineOperator,
131 criterion, parms, symmetric, pinfo)
161 std::size_t levels();
163 std::size_t maxlevels();
192 void createHierarchies(C& criterion,
193 std::shared_ptr<const Operator> fineOperator,
215 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
219 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
239 void mgc(LevelContext& levelContext,
Domain& x,
const Range& b);
247 void presmooth(LevelContext& levelContext,
Domain& x,
const Range& b);
255 void postsmooth(LevelContext& levelContext,
Domain& x,
const Range& b);
263 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel,
270 bool moveToCoarseLevel(LevelContext& levelContext);
276 void initIteratorsWithFineLevel(LevelContext& levelContext);
279 std::shared_ptr<OperatorHierarchy> matrices_;
281 std::shared_ptr<CoarseSolver> solver_;
283 std::shared_ptr<Hierarchy<Range,A>> rhs_;
285 std::shared_ptr<Hierarchy<Domain,A>> lhs_;
287 std::shared_ptr<Hierarchy<Domain,A>> residual_;
292 std::shared_ptr<ScalarProduct> scalarProduct_;
296 std::size_t preSteps_;
298 std::size_t postSteps_;
300 bool buildHierarchy_;
302 bool coarsesolverconverged;
304 typedef std::shared_ptr<Smoother> SmootherPointer;
305 SmootherPointer coarseSmoother_;
307 std::size_t verbosity_;
310 template<
class M,
class X,
class PI,
class A>
312 : matrices_(amg.matrices_), solver_(amg.solver_),
313 rhs_(), lhs_(), residual_(), scalarProduct_(amg.scalarProduct_),
314 gamma_(amg.gamma_), preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
315 symmetric(amg.symmetric), coarsesolverconverged(amg.coarsesolverconverged),
316 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
319 template<
class M,
class X,
class PI,
class A>
323 rhs_(), lhs_(), residual_(), scalarProduct_(),
324 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
325 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
326 symmetric(symmetric_), coarsesolverconverged(true),
327 coarseSmoother_(), verbosity_(parms.debugLevel())
329 if(preSteps_>1||postSteps_>1)
331 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
332 preSteps_=postSteps_=0;
334 assert(matrices_->isBuilt());
335 static_assert(std::is_same<PI,SequentialInformation>::value,
336 "Currently only sequential runs are supported");
338 template<
class M,
class X,
class PI,
class A>
345 : solver_(), rhs_(), lhs_(), residual_(), scalarProduct_(), gamma_(parms.getGamma()),
346 preSteps_(parms.getNoPreSmoothSteps()), postSteps_(parms.getNoPostSmoothSteps()),
347 buildHierarchy_(true),
348 symmetric(symmetric_), coarsesolverconverged(true),
349 coarseSmoother_(), verbosity_(criterion.debugLevel())
351 if(preSteps_>1||postSteps_>1)
353 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
354 preSteps_=postSteps_=1;
356 static_assert(std::is_same<PI,SequentialInformation>::value,
357 "Currently only sequential runs are supported");
361 createHierarchies(criterion, std::move(fineOperator), pinfo);
364 template<
class M,
class X,
class PI,
class A>
367 std::shared_ptr<const Operator> fineOperator,
371 matrices_ = std::make_shared<OperatorHierarchy>(
372 std::const_pointer_cast<Operator>(std::move(fineOperator)),
375 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
377 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
378 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.
elapsed()<<
" seconds."<<std::endl;
380 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
384 sargs.iterations = 1;
387 cargs.setArgs(sargs);
388 if(matrices_->redistributeInformation().back().isSetup()) {
390 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
391 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
393 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
394 cargs.setComm(*matrices_->parallelInformation().coarsest());
398 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
400#if HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK
401#if HAVE_SUITESPARSE_UMFPACK
402#define DIRECTSOLVER UMFPack
404#define DIRECTSOLVER SuperLU
407 if(std::is_same<ParallelInformation,SequentialInformation>::value
408 || matrices_->parallelInformation().coarsest()->communicator().size()==1
409 || (matrices_->parallelInformation().coarsest().isRedistributed()
410 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
411 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
412 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
413 std::cout<<
"Using superlu"<<std::endl;
414 if(matrices_->parallelInformation().coarsest().isRedistributed())
416 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
418 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
422 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
427 if(matrices_->parallelInformation().coarsest().isRedistributed())
429 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
431 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
433 *coarseSmoother_, 1E-2, 1000, 0));
437 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
439 *coarseSmoother_, 1E-2, 1000, 0));
443 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
444 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.
elapsed()<<
" seconds."<<std::endl;
448 template<
class M,
class X,
class PI,
class A>
456 typedef typename M::matrix_type
Matrix;
463 const Matrix& mat=matrices_->matrices().finest()->getmat();
464 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
465 bool isDirichlet =
true;
466 bool hasDiagonal =
false;
468 for(ColIter col=row->begin(); col!=row->end(); ++col) {
469 if(row.index()==col.index()) {
471 hasDiagonal = (*col != zero);
477 if(isDirichlet && hasDiagonal)
480 x[row.index()] = b[row.index()]/(*diag);
482 diag->solve(x[row.index()], b[row.index()]);
486 std::cout<<
" Preprocessing Dirichlet took "<<watch1.
elapsed()<<std::endl;
489 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
490 rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
491 lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
492 residual_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
493 matrices_->coarsenVector(*rhs_);
494 matrices_->coarsenVector(*lhs_);
495 matrices_->coarsenVector(*residual_);
502 template<
class M,
class X,
class PI,
class A>
505 return matrices_->levels();
507 template<
class M,
class X,
class PI,
class A>
508 std::size_t FastAMG<M,X,PI,A>::maxlevels()
510 return matrices_->maxlevels();
514 template<
class M,
class X,
class PI,
class A>
517 LevelContext levelContext;
519 initIteratorsWithFineLevel(levelContext);
521 assert(v.two_norm()==0);
524 if(matrices_->maxlevels()==1){
527 mgc(levelContext, v, b);
529 mgc(levelContext, v, d);
530 if(postSteps_==0||matrices_->maxlevels()==1)
531 levelContext.pinfo->copyOwnerToAll(v, v);
534 template<
class M,
class X,
class PI,
class A>
537 levelContext.matrix = matrices_->matrices().finest();
538 levelContext.pinfo = matrices_->parallelInformation().finest();
539 levelContext.redist =
540 matrices_->redistributeInformation().begin();
541 levelContext.aggregates = matrices_->aggregatesMaps().begin();
542 levelContext.lhs = lhs_->finest();
543 levelContext.residual = residual_->finest();
544 levelContext.rhs = rhs_->finest();
545 levelContext.level=0;
548 template<
class M,
class X,
class PI,
class A>
549 bool FastAMG<M,X,PI,A>
550 ::moveToCoarseLevel(LevelContext& levelContext)
552 bool processNextLevel=
true;
554 if(levelContext.redist->isSetup()) {
556 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.residual),
557 levelContext.residual.getRedistributed());
558 processNextLevel = levelContext.residual.getRedistributed().size()>0;
559 if(processNextLevel) {
561 ++levelContext.pinfo;
562 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
563 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
564 static_cast<const Range&
>(levelContext.residual.getRedistributed()),
565 *levelContext.pinfo);
570 ++levelContext.pinfo;
571 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
572 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
573 static_cast<const Range&
>(*levelContext.residual), *levelContext.pinfo);
576 if(processNextLevel) {
578 ++levelContext.residual;
580 ++levelContext.matrix;
581 ++levelContext.level;
582 ++levelContext.redist;
584 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
586 ++levelContext.aggregates;
590 *levelContext.residual=0;
592 return processNextLevel;
595 template<
class M,
class X,
class PI,
class A>
596 void FastAMG<M,X,PI,A>
597 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel, Domain& x)
599 if(processNextLevel) {
600 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
602 --levelContext.aggregates;
604 --levelContext.redist;
605 --levelContext.level;
607 --levelContext.matrix;
608 --levelContext.residual;
613 if(levelContext.redist->isSetup()) {
616 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
617 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
618 levelContext.lhs.getRedistributed(),
619 matrices_->getProlongationDampingFactor(),
620 *levelContext.pinfo, *levelContext.redist);
622 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
623 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
624 matrices_->getProlongationDampingFactor(), *levelContext.pinfo);
630 if(processNextLevel) {
637 template<
class M,
class X,
class PI,
class A>
638 void FastAMG<M,X,PI,A>
639 ::presmooth(LevelContext& levelContext, Domain& x,
const Range& b)
641 constexpr auto bl = blockLevel<typename M::matrix_type>();
642 GaussSeidelPresmoothDefect<bl>::apply(levelContext.matrix->getmat(),
644 *levelContext.residual,
648 template<
class M,
class X,
class PI,
class A>
649 void FastAMG<M,X,PI,A>
650 ::postsmooth(LevelContext& levelContext, Domain& x,
const Range& b)
652 constexpr auto bl = blockLevel<typename M::matrix_type>();
653 GaussSeidelPostsmoothDefect<bl>
654 ::apply(levelContext.matrix->getmat(), x, *levelContext.residual, b);
658 template<
class M,
class X,
class PI,
class A>
661 return IsDirectSolver< CoarseSolver>::value;
664 template<
class M,
class X,
class PI,
class A>
667 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
671 if(levelContext.redist->isSetup()) {
672 levelContext.redist->redistribute(b, levelContext.rhs.getRedistributed());
673 if(levelContext.rhs.getRedistributed().size()>0) {
675 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
676 levelContext.rhs.getRedistributed());
677 solver_->apply(levelContext.lhs.getRedistributed(), levelContext.rhs.getRedistributed(), res);
679 levelContext.redist->redistributeBackward(v, levelContext.lhs.getRedistributed());
680 levelContext.pinfo->copyOwnerToAll(v, v);
682 levelContext.pinfo->copyOwnerToAll(b, b);
683 solver_->apply(v,
const_cast<Range&
>(b), res);
689 coarsesolverconverged =
false;
695#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
696 bool processNextLevel = moveToCoarseLevel(levelContext);
698 if(processNextLevel) {
700 for(std::size_t i=0; i<gamma_; i++)
701 mgc(levelContext, *levelContext.lhs, *levelContext.rhs);
704 moveToFineLevel(levelContext, processNextLevel, v);
709 if(levelContext.matrix == matrices_->matrices().finest()) {
710 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
711 if(!coarsesolverconverged)
712 DUNE_THROW(MathError,
"Coarse solver did not converge");
721 template<
class M,
class X,
class PI,
class A>
729 template<
class M,
class X,
class PI,
class A>
733 matrices_->getCoarsestAggregatesOnFinest(cont);
A fast (sequential) algebraic multigrid based on agglomeration that saves memory bandwidth.
Definition: fastamg.hh:60
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:220
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:223
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:61
All parameters for AMG.
Definition: parameters.hh:416
ConstIterator class for sequential access.
Definition: matrix.hh:404
A generic dynamic dense matrix.
Definition: matrix.hh:561
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
T block_type
Export the type representing the components.
Definition: matrix.hh:568
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:96
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:33
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Sequential SSOR preconditioner.
Definition: preconditioners.hh:142
A simple stop watch.
Definition: timer.hh:31
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:47
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:67
A few common exception classes.
Some generic functions for pretty printing vectors and matrices.
Define base class for scalar product and norm.
Classes for using UMFPack with ISTL matrices.
#define DUNE_THROW(E,...)
Definition: exceptions.hh:312
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: fastamg.hh:207
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: fastamg.hh:731
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: fastamg.hh:211
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: fastamg.hh:173
void post(Domain &x)
Clean up.
Definition: fastamg.hh:722
X Domain
The domain type.
Definition: fastamg.hh:77
Hierarchy< Domain, A >::Iterator residual
The iterator over the residuals.
Definition: fastamg.hh:227
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: fastamg.hh:146
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: fastamg.hh:72
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: fastamg.hh:215
X Range
The range type.
Definition: fastamg.hh:79
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: fastamg.hh:70
M Operator
The matrix operator type.
Definition: fastamg.hh:63
FastAMG(const Operator &fineOperator, const C &criterion, const Parameters &parms=Parameters(), bool symmetric=true, const ParallelInformation &pinfo=ParallelInformation())
Construct an AMG with an inexact coarse solver based on the smoother.
Definition: fastamg.hh:125
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: fastamg.hh:81
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: fastamg.hh:659
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: fastamg.hh:223
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: fastamg.hh:231
std::size_t level
The level index.
Definition: fastamg.hh:235
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: fastamg.hh:515
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: fastamg.hh:74
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: fastamg.hh:449
FastAMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const Parameters &parms, bool symmetric=true)
Construct a new amg with a specific coarse solver.
Definition: fastamg.hh:320
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: fastamg.hh:219
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:406
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:428
Provides a classes representing the hierarchies in AMG.
Dune namespace.
Definition: alignedallocator.hh:13
std::shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:72
Define general preconditioner interface.
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:66
Statistics about the application of an inverse operator.
Definition: solver.hh:50
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Category
Definition: solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:25
Classes for using SuperLU with ISTL matrices.
Prolongation and restriction for amg.
Traits for type conversions and type information.