40 template<
class M,
class X,
class S,
class P,
class K,
class A>
56 template<
class M,
class X,
class S,
class PI=SequentialInformation,
57 class A=std::allocator<X> >
60 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
153 std::size_t levels();
155 std::size_t maxlevels();
184 void createHierarchies(C& criterion,
Operator& matrix,
210 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
214 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
238 void mgc(LevelContext& levelContext);
248 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
254 bool moveToCoarseLevel(LevelContext& levelContext);
260 void initIteratorsWithFineLevel(LevelContext& levelContext);
263 std::shared_ptr<OperatorHierarchy> matrices_;
267 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
269 std::shared_ptr<CoarseSolver> solver_;
279 std::shared_ptr<ScalarProduct> scalarProduct_;
283 std::size_t preSteps_;
285 std::size_t postSteps_;
286 bool buildHierarchy_;
288 bool coarsesolverconverged;
289 std::shared_ptr<Smoother> coarseSmoother_;
293 std::size_t verbosity_;
296 template<
class M,
class X,
class S,
class PI,
class A>
298 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
299 smoothers_(amg.smoothers_), solver_(amg.solver_),
300 rhs_(), lhs_(), update_(),
301 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
302 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
303 buildHierarchy_(amg.buildHierarchy_),
304 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
305 coarseSmoother_(amg.coarseSmoother_),
306 category_(amg.category_),
307 verbosity_(amg.verbosity_)
317 template<
class M,
class X,
class S,
class PI,
class A>
323 rhs_(), lhs_(), update_(), scalarProduct_(0),
324 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
325 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
326 additive(parms.getAdditive()), coarsesolverconverged(true),
330 verbosity_(parms.debugLevel())
332 assert(matrices_->isBuilt());
335 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
338 template<
class M,
class X,
class S,
class PI,
class A>
344 : smootherArgs_(smootherArgs),
346 rhs_(), lhs_(), update_(), scalarProduct_(),
347 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
348 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
349 additive(criterion.getAdditive()), coarsesolverconverged(true),
352 verbosity_(criterion.debugLevel())
355 DUNE_THROW(InvalidSolverCategory,
"Matrix and Communication must have the same SolverCategory!");
359 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
363 template<
class M,
class X,
class S,
class PI,
class A>
366 if(buildHierarchy_) {
370 coarseSmoother_.reset();
383 template <
class Matrix,
385 struct DirectSolverSelector
388 enum SolverType { umfpack, superlu,
none };
390 static constexpr SolverType solver =
391#if DISABLE_AMG_DIRECTSOLVER
393#elif HAVE_SUITESPARSE_UMFPACK
394 UMFPackMethodChooser< field_type > :: valid ? umfpack :
none ;
401 template <
class M, SolverType>
404 typedef InverseOperator<Vector,Vector> type;
405 static type* create(
const M& mat,
bool verbose,
bool reusevector )
407 DUNE_THROW(NotImplemented,
"DirectSolver not selected");
410 static std::string name () {
return "None"; }
412#if HAVE_SUITESPARSE_UMFPACK
414 struct Solver< M, umfpack >
416 typedef UMFPack< M > type;
417 static type* create(
const M& mat,
bool verbose,
bool reusevector )
419 return new type(mat, verbose, reusevector );
421 static std::string name () {
return "UMFPack"; }
426 struct Solver< M, superlu >
428 typedef SuperLU< M > type;
429 static type* create(
const M& mat,
bool verbose,
bool reusevector )
431 return new type(mat, verbose, reusevector );
433 static std::string name () {
return "SuperLU"; }
438 typedef Solver< Matrix, solver > SelectedSolver ;
439 typedef typename SelectedSolver :: type DirectSolver;
440 static constexpr bool isDirectSolver = solver !=
none;
441 static std::string name() {
return SelectedSolver :: name (); }
442 static DirectSolver* create(
const Matrix& mat,
bool verbose,
bool reusevector )
444 return SelectedSolver :: create( mat, verbose, reusevector );
448 template<
class M,
class X,
class S,
class PI,
class A>
450 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
454 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
456 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
459 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
465 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
466 && ( ! matrices_->redistributeInformation().back().isSetup() ||
467 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
470 SmootherArgs sargs(smootherArgs_);
471 sargs.iterations = 1;
474 cargs.setArgs(sargs);
475 if(matrices_->redistributeInformation().back().isSetup()) {
477 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
478 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
480 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
481 cargs.setComm(*matrices_->parallelInformation().coarsest());
485 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
487 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
490 if( SolverSelector::isDirectSolver &&
491 (std::is_same<ParallelInformation,SequentialInformation>::value
492 || matrices_->parallelInformation().coarsest()->communicator().size()==1
493 || (matrices_->parallelInformation().coarsest().isRedistributed()
494 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
495 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
498 if(matrices_->parallelInformation().coarsest().isRedistributed())
500 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
503 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
510 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
512 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
513 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
517 if(matrices_->parallelInformation().coarsest().isRedistributed())
519 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
524 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
526 *coarseSmoother_, 1E-2, 1000, 0));
531 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
533 *coarseSmoother_, 1E-2, 1000, 0));
552 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
553 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
554 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
558 template<
class M,
class X,
class S,
class PI,
class A>
565 typedef typename M::matrix_type
Matrix;
572 const Matrix& mat=matrices_->matrices().finest()->getmat();
573 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
574 bool isDirichlet =
true;
575 bool hasDiagonal =
false;
577 for(ColIter col=row->begin(); col!=row->end(); ++col) {
578 if(row.index()==col.index()) {
586 if(isDirichlet && hasDiagonal)
587 diagonal.solve(x[row.index()], b[row.index()]);
590 if(smoothers_->levels()>0)
591 smoothers_->finest()->pre(x,b);
594 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
607 matrices_->coarsenVector(*rhs_);
608 matrices_->coarsenVector(*lhs_);
609 matrices_->coarsenVector(*update_);
615 Iterator coarsest = smoothers_->
coarsest();
616 Iterator smoother = smoothers_->finest();
617 RIterator rhs = rhs_->finest();
618 DIterator lhs = lhs_->finest();
619 if(smoothers_->levels()>0) {
621 assert(lhs_->levels()==rhs_->levels());
622 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
623 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
625 if(smoother!=coarsest)
626 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
627 smoother->pre(*lhs,*rhs);
628 smoother->pre(*lhs,*rhs);
638 template<
class M,
class X,
class S,
class PI,
class A>
641 return matrices_->levels();
643 template<
class M,
class X,
class S,
class PI,
class A>
644 std::size_t AMG<M,X,S,PI,A>::maxlevels()
646 return matrices_->maxlevels();
650 template<
class M,
class X,
class S,
class PI,
class A>
653 LevelContext levelContext;
661 initIteratorsWithFineLevel(levelContext);
664 *levelContext.lhs = v;
665 *levelContext.rhs = d;
666 *levelContext.update=0;
667 levelContext.level=0;
671 if(postSteps_==0||matrices_->maxlevels()==1)
672 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
674 v=*levelContext.update;
679 template<
class M,
class X,
class S,
class PI,
class A>
682 levelContext.smoother = smoothers_->finest();
683 levelContext.matrix = matrices_->matrices().finest();
684 levelContext.pinfo = matrices_->parallelInformation().finest();
685 levelContext.redist =
686 matrices_->redistributeInformation().begin();
687 levelContext.aggregates = matrices_->aggregatesMaps().begin();
688 levelContext.lhs = lhs_->finest();
689 levelContext.update = update_->finest();
690 levelContext.rhs = rhs_->finest();
693 template<
class M,
class X,
class S,
class PI,
class A>
695 ::moveToCoarseLevel(LevelContext& levelContext)
698 bool processNextLevel=
true;
700 if(levelContext.redist->isSetup()) {
701 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
702 levelContext.rhs.getRedistributed());
703 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
704 if(processNextLevel) {
707 ++levelContext.pinfo;
708 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
709 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
710 static_cast<const Range&
>(fineRhs.getRedistributed()),
711 *levelContext.pinfo);
716 ++levelContext.pinfo;
717 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
718 ::restrictVector(*(*levelContext.aggregates),
719 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
720 *levelContext.pinfo);
723 if(processNextLevel) {
726 ++levelContext.update;
727 ++levelContext.matrix;
728 ++levelContext.level;
729 ++levelContext.redist;
731 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
733 ++levelContext.smoother;
734 ++levelContext.aggregates;
737 *levelContext.update=0;
739 return processNextLevel;
742 template<
class M,
class X,
class S,
class PI,
class A>
744 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
746 if(processNextLevel) {
747 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
749 --levelContext.smoother;
750 --levelContext.aggregates;
752 --levelContext.redist;
753 --levelContext.level;
755 --levelContext.matrix;
759 --levelContext.pinfo;
761 if(levelContext.redist->isSetup()) {
763 levelContext.lhs.getRedistributed()=0;
764 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
765 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
766 levelContext.lhs.getRedistributed(),
767 matrices_->getProlongationDampingFactor(),
768 *levelContext.pinfo, *levelContext.redist);
771 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
772 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
773 matrices_->getProlongationDampingFactor(),
774 *levelContext.pinfo);
778 if(processNextLevel) {
779 --levelContext.update;
783 *levelContext.update += *levelContext.lhs;
786 template<
class M,
class X,
class S,
class PI,
class A>
789 return IsDirectSolver< CoarseSolver>::value;
792 template<
class M,
class X,
class S,
class PI,
class A>
794 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
798 if(levelContext.redist->isSetup()) {
799 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
800 if(levelContext.rhs.getRedistributed().size()>0) {
802 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
803 levelContext.rhs.getRedistributed());
804 solver_->apply(levelContext.update.getRedistributed(),
805 levelContext.rhs.getRedistributed(), res);
807 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
808 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
810 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
811 solver_->apply(*levelContext.update, *levelContext.rhs, res);
815 coarsesolverconverged =
false;
820#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
821 bool processNextLevel = moveToCoarseLevel(levelContext);
823 if(processNextLevel) {
825 for(std::size_t i=0; i<gamma_; i++)
829 moveToFineLevel(levelContext, processNextLevel);
834 if(levelContext.matrix == matrices_->matrices().finest()) {
835 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
836 if(!coarsesolverconverged)
837 DUNE_THROW(MathError,
"Coarse solver did not converge");
845 template<
class M,
class X,
class S,
class PI,
class A>
846 void AMG<M,X,S,PI,A>::additiveMgc(){
849 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
852 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
856 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
857 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
863 lhs = lhs_->finest();
866 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
869 smoother->apply(*lhs, *rhs);
873#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
874 InverseOperatorResult res;
875 pinfo->copyOwnerToAll(*rhs, *rhs);
876 solver_->apply(*lhs, *rhs, res);
879 DUNE_THROW(MathError,
"Coarse solver did not converge");
888 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
889 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
895 template<
class M,
class X,
class S,
class PI,
class A>
902 Iterator coarsest = smoothers_->
coarsest();
903 Iterator smoother = smoothers_->finest();
904 DIterator lhs = lhs_->finest();
905 if(smoothers_->levels()>0) {
906 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
907 smoother->post(*lhs);
908 if(smoother!=coarsest)
909 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
910 smoother->post(*lhs);
911 smoother->post(*lhs);
924 template<
class M,
class X,
class S,
class PI,
class A>
928 matrices_->getCoarsestAggregatesOnFinest(cont);
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:59
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:250
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:253
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:142
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:31
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:310
All parameters for AMG.
Definition: parameters.hh:391
ConstIterator class for sequential access.
Definition: matrix.hh:398
A generic dynamic dense matrix.
Definition: matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
T block_type
Export the type representing the components.
Definition: matrix.hh:562
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:95
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
A few common exception classes.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr GeometryType none(unsigned int dim)
Returns a GeometryType representing a singular of dimension dim.
Definition: type.hh:714
AMG(const AMG &amg)
Copy constructor.
Definition: amg.hh:297
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:318
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:559
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:222
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:226
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:787
X Domain
The domain type.
Definition: amg.hh:81
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:206
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:214
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:94
S Smoother
The type of the smoother.
Definition: amg.hh:91
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:198
M Operator
The matrix operator type.
Definition: amg.hh:67
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:202
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:210
X Range
The range type.
Definition: amg.hh:83
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:454
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:926
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:218
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:165
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1369
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:78
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:85
void post(Domain &x)
Clean up.
Definition: amg.hh:896
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:476
static T * construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
std::size_t level
The level index.
Definition: amg.hh:230
AMG(const Operator &fineOperator, const C &criterion, const SmootherArgs &smootherArgs=SmootherArgs(), const ParallelInformation &pinfo=ParallelInformation())
Construct an AMG with an inexact coarse solver based on the smoother.
Definition: amg.hh:340
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:651
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:76
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amg.hh:138
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:74
Provides a classes representing the hierarchies in AMG.
Dune namespace.
Definition: alignedallocator.hh:10
shared_ptr< T > stackobject_to_shared_ptr(T &t)
Create a shared_ptr for a stack-allocated object.
Definition: shared_ptr.hh:75
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:41
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
Categories for the solvers.
Definition: solvercategory.hh:20
Category
Definition: solvercategory.hh:21
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:32
Classes for using SuperLU with ISTL matrices.
Prolongation and restriction for amg.
Traits for type conversions and type information.
Classes for using UMFPack with ISTL matrices.