4 #ifndef DUNE_AMG_AMG_HH
5 #define DUNE_AMG_AMG_HH
8 #include <dune/common/exceptions.hh>
17 #include <dune/common/typetraits.hh>
18 #include <dune/common/exceptions.hh>
41 template<
class M,
class X,
class S,
class P,
class K,
class A>
55 class A=std::allocator<X> >
58 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
117 std::size_t preSmoothingSteps,
118 std::size_t postSmoothingSteps,
119 bool additive=
false) DUNE_DEPRECATED;
153 AMG(const
Operator& fineOperator, const C& criterion,
155 std::
size_t preSmoothingSteps,
156 std::
size_t postSmoothingSteps,
172 AMG(const
Operator& fineOperator, const C& criterion,
213 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
230 void createHierarchies(C& criterion,
Operator& matrix,
256 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
260 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
284 void mgc(LevelContext& levelContext);
294 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
300 bool moveToCoarseLevel(LevelContext& levelContext);
306 void initIteratorsWithFineLevel(LevelContext& levelContext);
309 Dune::shared_ptr<OperatorHierarchy> matrices_;
313 Dune::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
315 Dune::shared_ptr<CoarseSolver> solver_;
325 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
326 typedef Dune::shared_ptr<ScalarProduct> ScalarProductPointer;
328 ScalarProductPointer scalarProduct_;
332 std::size_t preSteps_;
334 std::size_t postSteps_;
335 bool buildHierarchy_;
337 bool coarsesolverconverged;
338 Dune::shared_ptr<Smoother> coarseSmoother_;
340 std::size_t verbosity_;
343 template<
class M,
class X,
class S,
class PI,
class A>
345 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
346 smoothers_(amg.smoothers_), solver_(amg.solver_),
347 rhs_(), lhs_(), update_(),
348 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
349 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
350 buildHierarchy_(amg.buildHierarchy_),
351 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
352 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
362 template<
class M,
class X,
class S,
class PI,
class A>
365 std::size_t gamma, std::size_t preSmoothingSteps,
366 std::size_t postSmoothingSteps,
bool additive_)
367 : matrices_(&matrices), smootherArgs_(smootherArgs),
369 rhs_(), lhs_(), update_(), scalarProduct_(0),
370 gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
371 additive(additive_), coarsesolverconverged(true),
372 coarseSmoother_(), verbosity_(2)
374 assert(matrices_->isBuilt());
377 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
380 template<
class M,
class X,
class S,
class PI,
class A>
384 : matrices_(&matrices), smootherArgs_(smootherArgs),
386 rhs_(), lhs_(), update_(), scalarProduct_(0),
387 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
388 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
389 additive(parms.getAdditive()), coarsesolverconverged(true),
390 coarseSmoother_(), verbosity_(parms.debugLevel())
392 assert(matrices_->isBuilt());
395 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
398 template<
class M,
class X,
class S,
class PI,
class A>
403 std::size_t gamma, std::size_t preSmoothingSteps,
404 std::size_t postSmoothingSteps,
407 : smootherArgs_(smootherArgs),
409 rhs_(), lhs_(), update_(), scalarProduct_(), gamma_(gamma),
410 preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
411 additive(additive_), coarsesolverconverged(true),
412 coarseSmoother_(), verbosity_(criterion.debugLevel())
414 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
415 "Matrix and Solver must match in terms of category!");
419 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
422 template<
class M,
class X,
class S,
class PI,
class A>
428 : smootherArgs_(smootherArgs),
430 rhs_(), lhs_(), update_(), scalarProduct_(),
431 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
432 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
433 additive(criterion.getAdditive()), coarsesolverconverged(true),
434 coarseSmoother_(), verbosity_(criterion.debugLevel())
436 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
437 "Matrix and Solver must match in terms of category!");
441 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
445 template<
class M,
class X,
class S,
class PI,
class A>
448 if(buildHierarchy_) {
452 coarseSmoother_.reset();
465 template<
class M,
class X,
class S,
class PI,
class A>
471 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
473 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
476 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
478 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
480 SmootherArgs sargs(smootherArgs_);
481 sargs.iterations = 1;
484 cargs.setArgs(sargs);
485 if(matrices_->redistributeInformation().back().isSetup()) {
487 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
488 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
490 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
491 cargs.setComm(*matrices_->parallelInformation().coarsest());
494 coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
495 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
497 #if HAVE_SUPERLU || HAVE_UMFPACK
499 #define DIRECTSOLVER UMFPack
501 #define DIRECTSOLVER SuperLU
504 if(is_same<ParallelInformation,SequentialInformation>::value
505 || matrices_->parallelInformation().coarsest()->communicator().size()==1
506 || (matrices_->parallelInformation().coarsest().isRedistributed()
507 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
508 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
509 if(matrices_->parallelInformation().coarsest().isRedistributed())
511 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
513 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
517 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
518 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
519 std::cout<<
"Using a direct coarse solver (" <<
static_cast< DIRECTSOLVER<typename M::matrix_type>*
>(solver_.get())->name() <<
")" << std::endl;
524 if(matrices_->parallelInformation().coarsest().isRedistributed())
526 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
528 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
530 *coarseSmoother_, 1E-2, 1000, 0));
534 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
536 *coarseSmoother_, 1E-2, 1000, 0));
540 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
541 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
542 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
546 template<
class M,
class X,
class S,
class PI,
class A>
553 typedef typename M::matrix_type
Matrix;
560 const Matrix&
mat=matrices_->matrices().finest()->getmat();
562 bool isDirichlet =
true;
563 bool hasDiagonal =
false;
566 if(
row.index()==
col.index()) {
574 if(isDirichlet && hasDiagonal)
575 diagonal.solve(x[
row.index()], b[
row.index()]);
578 if(smoothers_->levels()>0)
579 smoothers_->finest()->pre(x,b);
582 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
595 matrices_->coarsenVector(*rhs_);
596 matrices_->coarsenVector(*lhs_);
597 matrices_->coarsenVector(*update_);
603 Iterator coarsest = smoothers_->
coarsest();
604 Iterator smoother = smoothers_->finest();
605 RIterator rhs = rhs_->finest();
606 DIterator lhs = lhs_->finest();
607 if(smoothers_->levels()>0) {
609 assert(lhs_->levels()==rhs_->levels());
610 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
611 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
613 if(smoother!=coarsest)
614 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
615 smoother->pre(*lhs,*rhs);
616 smoother->pre(*lhs,*rhs);
626 template<
class M,
class X,
class S,
class PI,
class A>
629 return matrices_->levels();
631 template<
class M,
class X,
class S,
class PI,
class A>
634 return matrices_->maxlevels();
638 template<
class M,
class X,
class S,
class PI,
class A>
641 LevelContext levelContext;
649 initIteratorsWithFineLevel(levelContext);
652 *levelContext.lhs = v;
653 *levelContext.rhs = d;
654 *levelContext.update=0;
655 levelContext.level=0;
659 if(postSteps_==0||matrices_->maxlevels()==1)
660 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
662 v=*levelContext.update;
667 template<
class M,
class X,
class S,
class PI,
class A>
670 levelContext.smoother = smoothers_->finest();
671 levelContext.matrix = matrices_->matrices().finest();
672 levelContext.pinfo = matrices_->parallelInformation().finest();
673 levelContext.redist =
674 matrices_->redistributeInformation().begin();
675 levelContext.aggregates = matrices_->aggregatesMaps().begin();
676 levelContext.lhs = lhs_->finest();
677 levelContext.update = update_->finest();
678 levelContext.rhs = rhs_->finest();
681 template<
class M,
class X,
class S,
class PI,
class A>
683 ::moveToCoarseLevel(LevelContext& levelContext)
686 bool processNextLevel=
true;
688 if(levelContext.redist->isSetup()) {
689 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
690 levelContext.rhs.getRedistributed());
691 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
692 if(processNextLevel) {
695 ++levelContext.pinfo;
698 static_cast<const Range&>(fineRhs.getRedistributed()),
699 *levelContext.pinfo);
704 ++levelContext.pinfo;
707 *levelContext.rhs, static_cast<const Range&>(*fineRhs),
708 *levelContext.pinfo);
711 if(processNextLevel) {
714 ++levelContext.update;
715 ++levelContext.matrix;
716 ++levelContext.level;
717 ++levelContext.redist;
719 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
721 ++levelContext.smoother;
722 ++levelContext.aggregates;
725 *levelContext.update=0;
727 return processNextLevel;
730 template<
class M,
class X,
class S,
class PI,
class A>
732 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
734 if(processNextLevel) {
735 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
737 --levelContext.smoother;
738 --levelContext.aggregates;
740 --levelContext.redist;
741 --levelContext.level;
743 --levelContext.matrix;
747 --levelContext.pinfo;
749 if(levelContext.redist->isSetup()) {
751 levelContext.lhs.getRedistributed()=0;
753 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
754 levelContext.lhs.getRedistributed(),
755 matrices_->getProlongationDampingFactor(),
756 *levelContext.pinfo, *levelContext.redist);
760 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
761 matrices_->getProlongationDampingFactor(),
762 *levelContext.pinfo);
766 if(processNextLevel) {
767 --levelContext.update;
771 *levelContext.update += *levelContext.lhs;
774 template<
class M,
class X,
class S,
class PI,
class A>
780 template<
class M,
class X,
class S,
class PI,
class A>
782 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
786 if(levelContext.redist->isSetup()) {
787 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
788 if(levelContext.rhs.getRedistributed().size()>0) {
790 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
791 levelContext.rhs.getRedistributed());
792 solver_->apply(levelContext.update.getRedistributed(),
793 levelContext.rhs.getRedistributed(), res);
795 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
796 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
798 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
799 solver_->apply(*levelContext.update, *levelContext.rhs, res);
803 coarsesolverconverged =
false;
808 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
809 bool processNextLevel = moveToCoarseLevel(levelContext);
811 if(processNextLevel) {
813 for(std::size_t i=0; i<gamma_; i++)
817 moveToFineLevel(levelContext, processNextLevel);
822 if(levelContext.matrix == matrices_->matrices().finest()) {
823 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
824 if(!coarsesolverconverged)
825 DUNE_THROW(MathError,
"Coarse solver did not converge");
833 template<
class M,
class X,
class S,
class PI,
class A>
834 void AMG<M,X,S,PI,A>::additiveMgc(){
837 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
840 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
845 ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
851 lhs = lhs_->finest();
854 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
857 smoother->apply(*lhs, *rhs);
861 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
862 InverseOperatorResult res;
863 pinfo->copyOwnerToAll(*rhs, *rhs);
864 solver_->apply(*lhs, *rhs, res);
867 DUNE_THROW(MathError,
"Coarse solver did not converge");
883 template<
class M,
class X,
class S,
class PI,
class A>
886 DUNE_UNUSED_PARAMETER(x);
890 Iterator coarsest = smoothers_->
coarsest();
891 Iterator smoother = smoothers_->finest();
892 DIterator lhs = lhs_->finest();
893 if(smoothers_->levels()>0) {
894 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
895 smoother->post(*lhs);
896 if(smoother!=coarsest)
897 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
898 smoother->post(*lhs);
899 smoother->post(*lhs);
912 template<
class M,
class X,
class S,
class PI,
class A>
916 matrices_->getCoarsestAggregatesOnFinest(cont);
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:914
The solver category.
Definition: amg.hh:96
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:92
void post(Domain &x)
Clean up.
Definition: amg.hh:884
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:26
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:252
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:83
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:72
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:76
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:261
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:547
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:639
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1386
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Two grid operator for AMG with Krylov cycle.
Definition: amg.hh:45
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Smoother SmootherType
Definition: amg.hh:240
Definition: solvertype.hh:13
S Smoother
The type of the smoother.
Definition: amg.hh:89
Templates characterizing the type of a solver.
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:775
A generic dynamic dense matrix.
Definition: matrix.hh:24
std::size_t level
The level index.
Definition: amg.hh:276
Prolongation and restriction for amg.
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:256
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:85
Col col
Definition: matrixmatrix.hh:347
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:45
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:260
an algebraic multigrid method using a Krylov-cycle.
Definition: amg.hh:42
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:74
Implementations of the inverse operator interface.
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:317
M Operator
The matrix operator type.
Definition: amg.hh:65
Matrix & A
Definition: matrixmatrix.hh:216
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:268
~AMG()
Definition: amg.hh:446
VariableBlockVector< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition: matrix.hh:50
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:63
Provides a classes representing the hierarchies in AMG.
Statistics about the application of an inverse operator.
Definition: solver.hh:31
Row row
Definition: matrixmatrix.hh:345
X Domain
The domain type.
Definition: amg.hh:79
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:272
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
T block_type
Export the type representing the components.
Definition: matrix.hh:32
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:264
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:258
Choose the approriate scalar product for a solver category.
Definition: scalarproducts.hh:75
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:56
X Range
The range type.
Definition: amg.hh:81
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:244
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:211
Classes for using UMFPack with ISTL matrices.
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:408
Define base class for scalar product and norm.
std::size_t levels()
Definition: amg.hh:627
Matrix & mat
Definition: matrixmatrix.hh:343
Classes for the generic construction and application of the smoothers.
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:430
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, bool additive=false)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:363
std::size_t maxlevels()
Definition: amg.hh:632
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:248
Classes for using SuperLU with ISTL matrices.
All parameters for AMG.
Definition: parameters.hh:390
Definition: example.cc:34