3 #ifndef DUNE_AMG_AMG_HH 4 #define DUNE_AMG_AMG_HH 7 #include <dune/common/exceptions.hh> 16 #include <dune/common/typetraits.hh> 17 #include <dune/common/exceptions.hh> 40 template<
class M,
class X,
class S,
class P,
class K,
class A>
54 class A=std::allocator<X> >
57 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
107 AMG(
const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
108 const SmootherArgs& smootherArgs,
const Parameters& parms);
122 AMG(
const Operator& fineOperator,
const C& criterion,
134 void pre(Domain& x, Range& b);
137 void apply(Domain& v,
const Range& d);
140 void post(Domain& x);
163 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
180 void createHierarchies(C& criterion, Operator& matrix,
206 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
210 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
234 void mgc(LevelContext& levelContext);
244 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
250 bool moveToCoarseLevel(LevelContext& levelContext);
256 void initIteratorsWithFineLevel(LevelContext& levelContext);
259 std::shared_ptr<OperatorHierarchy> matrices_;
261 SmootherArgs smootherArgs_;
263 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
265 std::shared_ptr<CoarseSolver> solver_;
275 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
276 typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
278 ScalarProductPointer scalarProduct_;
282 std::size_t preSteps_;
284 std::size_t postSteps_;
285 bool buildHierarchy_;
287 bool coarsesolverconverged;
288 std::shared_ptr<Smoother> coarseSmoother_;
290 std::size_t verbosity_;
293 template<
class M,
class X,
class S,
class PI,
class A>
295 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
296 smoothers_(amg.smoothers_), solver_(amg.solver_),
297 rhs_(), lhs_(), update_(),
298 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
299 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
300 buildHierarchy_(amg.buildHierarchy_),
301 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
302 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
312 template<
class M,
class X,
class S,
class PI,
class A>
314 const SmootherArgs& smootherArgs,
316 : matrices_(&matrices), smootherArgs_(smootherArgs),
317 smoothers_(new
Hierarchy<Smoother,A>), solver_(&coarseSolver),
318 rhs_(), lhs_(), update_(), scalarProduct_(0),
319 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
320 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
321 additive(parms.getAdditive()), coarsesolverconverged(true),
322 coarseSmoother_(), verbosity_(parms.debugLevel())
324 assert(matrices_->isBuilt());
327 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
330 template<
class M,
class X,
class S,
class PI,
class A>
334 const SmootherArgs& smootherArgs,
336 : smootherArgs_(smootherArgs),
337 smoothers_(new
Hierarchy<Smoother,A>), solver_(),
338 rhs_(), lhs_(), update_(), scalarProduct_(),
339 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
340 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
341 additive(criterion.getAdditive()), coarsesolverconverged(true),
342 coarseSmoother_(), verbosity_(criterion.debugLevel())
344 static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
345 "Matrix and Solver must match in terms of category!");
349 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
353 template<
class M,
class X,
class S,
class PI,
class A>
356 if(buildHierarchy_) {
360 coarseSmoother_.reset();
381 #if HAVE_SUITESPARSE_UMFPACK 389 template <
class M, SolverType>
393 static type*
create(
const M&
mat,
bool verbose,
bool reusevector )
395 DUNE_THROW(NotImplemented,
"DirectSolver not selected");
398 static std::string
name () {
return "None"; }
400 #if HAVE_SUITESPARSE_UMFPACK 402 struct Solver< M, umfpack >
405 static type* create(
const M&
mat,
bool verbose,
bool reusevector )
407 return new type(mat, verbose, reusevector );
409 static std::string name () {
return "UMFPack"; }
417 static type*
create(
const M&
mat,
bool verbose,
bool reusevector )
419 return new type(mat, verbose, reusevector );
421 static std::string
name () {
return "SuperLU"; }
428 static constexpr
bool isDirectSolver = solver != none;
429 static std::string
name() {
return SelectedSolver :: name (); }
432 return SelectedSolver :: create( mat, verbose, reusevector );
436 template<
class M,
class X,
class S,
class PI,
class A>
442 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
444 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
447 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
449 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels())
452 SmootherArgs sargs(smootherArgs_);
453 sargs.iterations = 1;
456 cargs.setArgs(sargs);
457 if(matrices_->redistributeInformation().back().isSetup()) {
459 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
460 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
462 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
463 cargs.setComm(*matrices_->parallelInformation().coarsest());
467 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
472 if( SolverSelector::isDirectSolver &&
473 (std::is_same<ParallelInformation,SequentialInformation>::value
474 || matrices_->parallelInformation().coarsest()->communicator().size()==1
475 || (matrices_->parallelInformation().coarsest().isRedistributed()
476 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
477 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
480 if(matrices_->parallelInformation().coarsest().isRedistributed())
482 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
485 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
492 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
494 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
495 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
499 if(matrices_->parallelInformation().coarsest().isRedistributed())
501 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
503 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
505 *coarseSmoother_, 1E-2, 1000, 0));
509 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
511 *coarseSmoother_, 1E-2, 1000, 0));
515 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
516 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels " 517 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
521 template<
class M,
class X,
class S,
class PI,
class A>
528 typedef typename M::matrix_type
Matrix;
535 const Matrix&
mat=matrices_->matrices().finest()->getmat();
536 for(RowIter row=mat.
begin(); row!=mat.
end(); ++row) {
537 bool isDirichlet =
true;
538 bool hasDiagonal =
false;
540 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
541 if(row.index()==
col.index()) {
549 if(isDirichlet && hasDiagonal)
550 diagonal.solve(x[row.index()], b[row.index()]);
553 if(smoothers_->levels()>0)
554 smoothers_->finest()->pre(x,b);
557 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
558 Range* copy =
new Range(b);
562 Domain* dcopy =
new Domain(x);
566 dcopy =
new Domain(x);
570 matrices_->coarsenVector(*rhs_);
571 matrices_->coarsenVector(*lhs_);
572 matrices_->coarsenVector(*update_);
578 Iterator coarsest = smoothers_->
coarsest();
579 Iterator smoother = smoothers_->finest();
580 RIterator rhs = rhs_->finest();
581 DIterator lhs = lhs_->finest();
582 if(smoothers_->levels()>0) {
584 assert(lhs_->levels()==rhs_->levels());
585 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
586 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
588 if(smoother!=coarsest)
589 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
590 smoother->pre(*lhs,*rhs);
591 smoother->pre(*lhs,*rhs);
601 template<
class M,
class X,
class S,
class PI,
class A>
604 return matrices_->
levels();
606 template<
class M,
class X,
class S,
class PI,
class A>
613 template<
class M,
class X,
class S,
class PI,
class A>
616 LevelContext levelContext;
624 initIteratorsWithFineLevel(levelContext);
627 *levelContext.lhs = v;
628 *levelContext.rhs = d;
629 *levelContext.update=0;
630 levelContext.level=0;
634 if(postSteps_==0||matrices_->maxlevels()==1)
635 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
637 v=*levelContext.update;
642 template<
class M,
class X,
class S,
class PI,
class A>
645 levelContext.smoother = smoothers_->finest();
646 levelContext.matrix = matrices_->matrices().finest();
647 levelContext.pinfo = matrices_->parallelInformation().finest();
648 levelContext.redist =
649 matrices_->redistributeInformation().begin();
650 levelContext.aggregates = matrices_->aggregatesMaps().begin();
651 levelContext.lhs = lhs_->finest();
652 levelContext.update = update_->finest();
653 levelContext.rhs = rhs_->finest();
656 template<
class M,
class X,
class S,
class PI,
class A>
661 bool processNextLevel=
true;
663 if(levelContext.redist->isSetup()) {
664 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
665 levelContext.rhs.getRedistributed());
666 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
667 if(processNextLevel) {
670 ++levelContext.pinfo;
673 static_cast<const Range&>(fineRhs.getRedistributed()),
674 *levelContext.pinfo);
679 ++levelContext.pinfo;
682 *levelContext.rhs, static_cast<const Range&>(*fineRhs),
683 *levelContext.pinfo);
686 if(processNextLevel) {
689 ++levelContext.update;
690 ++levelContext.matrix;
691 ++levelContext.level;
692 ++levelContext.redist;
694 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
696 ++levelContext.smoother;
697 ++levelContext.aggregates;
700 *levelContext.update=0;
702 return processNextLevel;
705 template<
class M,
class X,
class S,
class PI,
class A>
709 if(processNextLevel) {
710 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
712 --levelContext.smoother;
713 --levelContext.aggregates;
715 --levelContext.redist;
716 --levelContext.level;
718 --levelContext.matrix;
722 --levelContext.pinfo;
724 if(levelContext.redist->isSetup()) {
726 levelContext.lhs.getRedistributed()=0;
728 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
729 levelContext.lhs.getRedistributed(),
730 matrices_->getProlongationDampingFactor(),
731 *levelContext.pinfo, *levelContext.redist);
735 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
736 matrices_->getProlongationDampingFactor(),
737 *levelContext.pinfo);
741 if(processNextLevel) {
742 --levelContext.update;
746 *levelContext.update += *levelContext.lhs;
749 template<
class M,
class X,
class S,
class PI,
class A>
755 template<
class M,
class X,
class S,
class PI,
class A>
757 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
761 if(levelContext.redist->isSetup()) {
762 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
763 if(levelContext.rhs.getRedistributed().size()>0) {
765 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
766 levelContext.rhs.getRedistributed());
767 solver_->apply(levelContext.update.getRedistributed(),
768 levelContext.rhs.getRedistributed(), res);
770 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
771 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
773 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
774 solver_->apply(*levelContext.update, *levelContext.rhs, res);
778 coarsesolverconverged =
false;
783 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 784 bool processNextLevel = moveToCoarseLevel(levelContext);
786 if(processNextLevel) {
788 for(std::size_t i=0; i<gamma_; i++)
792 moveToFineLevel(levelContext, processNextLevel);
797 if(levelContext.matrix == matrices_->matrices().finest()) {
798 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
799 if(!coarsesolverconverged)
800 DUNE_THROW(MathError,
"Coarse solver did not converge");
808 template<
class M,
class X,
class S,
class PI,
class A>
815 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
820 ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
829 for(rhs=rhs_->
finest(); rhs != rhs_->
coarsest(); ++lhs, ++rhs, ++smoother) {
832 smoother->apply(*lhs, *rhs);
836 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 838 pinfo->copyOwnerToAll(*rhs, *rhs);
839 solver_->apply(*lhs, *rhs, res);
842 DUNE_THROW(MathError,
"Coarse solver did not converge");
858 template<
class M,
class X,
class S,
class PI,
class A>
861 DUNE_UNUSED_PARAMETER(x);
865 Iterator coarsest = smoothers_->
coarsest();
866 Iterator smoother = smoothers_->finest();
867 DIterator lhs = lhs_->finest();
868 if(smoothers_->levels()>0) {
869 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
870 smoother->post(*lhs);
871 if(smoother!=coarsest)
872 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
873 smoother->post(*lhs);
874 smoother->post(*lhs);
887 template<
class M,
class X,
class S,
class PI,
class A>
891 matrices_->getCoarsestAggregatesOnFinest(cont);
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:75
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:430
All parameters for AMG.
Definition: parameters.hh:390
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:417
Use the UMFPack package to directly solve linear systems – empty default class.
Definition: umfpack.hh:49
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:63
Col col
Definition: matrixmatrix.hh:349
X Range
The range type.
Definition: amg.hh:80
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:194
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:161
A generic dynamic dense matrix.
Definition: matrix.hh:554
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:91
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
Definition: example.cc:35
Definition: superlu.hh:59
Choose the approriate scalar product for a solver category.
Definition: scalarproducts.hh:76
Definition: transfer.hh:30
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:202
Implementations of the inverse operator interface.
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:55
Templates characterizing the type of a solver.
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:257
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:393
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:522
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
T block_type
Export the type representing the components.
Definition: matrix.hh:562
void post(Domain &x)
Clean up.
Definition: amg.hh:859
Prolongation and restriction for amg.
Two grid operator for AMG with Krylov cycle.
Definition: amg.hh:44
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:210
static std::string name()
Definition: amg.hh:429
static std::string name()
Definition: amg.hh:421
S Smoother
The type of the smoother.
Definition: amg.hh:88
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:260
Traits class for generically constructing non default constructable types.
Definition: novlpschwarz.hh:328
SelectedSolver ::type DirectSolver
Definition: amg.hh:427
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:218
Define base class for scalar product and norm.
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:574
~AMG()
Definition: amg.hh:354
Matrix ::field_type field_type
Definition: amg.hh:377
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:206
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:614
Matrix & mat
Definition: matrixmatrix.hh:345
Iterator finest()
Get an iterator positioned at the finest level.
Definition: hierarchy.hh:1379
Smoother SmootherType
Definition: amg.hh:190
Definition: basearray.hh:19
Solver< Matrix, solver > SelectedSolver
Definition: amg.hh:426
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
The solver category.
Definition: amg.hh:95
std::size_t level
The level index.
Definition: amg.hh:226
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:214
std::size_t maxlevels()
Definition: amg.hh:607
static std::string name()
Definition: amg.hh:398
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:198
Definition: solvertype.hh:13
Classes for the generic construction and application of the smoothers.
ConstIterator class for sequential access.
Definition: matrix.hh:397
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:889
SuperLU< M > type
Definition: amg.hh:416
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:408
Definition: umfpack.hh:55
SolverType
Definition: amg.hh:378
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:73
X Domain
The domain type.
Definition: amg.hh:78
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:222
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:82
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:313
M Operator
The matrix operator type.
Definition: amg.hh:64
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:430
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:750
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:316
Classes for using UMFPack with ISTL matrices.
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1385
an algebraic multigrid method using a Krylov-cycle.
Definition: amg.hh:41
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
Statistics about the application of an inverse operator.
Definition: solver.hh:31
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:26
Classes for using SuperLU with ISTL matrices.
std::size_t levels()
Definition: amg.hh:602
InverseOperator< Vector, Vector > type
Definition: amg.hh:392
Provides a classes representing the hierarchies in AMG.
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:71