40 template<
class M,
class X,
class S,
class P,
class K,
class A>
53 template<
class M,
class X,
class S,
class PI=SequentialInformation,
54 class A=std::allocator<X> >
57 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
116 std::size_t preSmoothingSteps,
117 std::size_t postSmoothingSteps,
154 std::
size_t preSmoothingSteps,
155 std::
size_t postSmoothingSteps,
198 std::
size_t levels();
200 std::
size_t maxlevels();
229 void createHierarchies(C& criterion,
Operator& matrix,
255 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
259 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
283 void mgc(LevelContext& levelContext);
293 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
299 bool moveToCoarseLevel(LevelContext& levelContext);
305 void initIteratorsWithFineLevel(LevelContext& levelContext);
308 std::shared_ptr<OperatorHierarchy> matrices_;
312 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
314 std::shared_ptr<CoarseSolver> solver_;
324 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
325 typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
327 ScalarProductPointer scalarProduct_;
331 std::size_t preSteps_;
333 std::size_t postSteps_;
334 bool buildHierarchy_;
336 bool coarsesolverconverged;
337 std::shared_ptr<Smoother> coarseSmoother_;
339 std::size_t verbosity_;
342 template<
class M,
class X,
class S,
class PI,
class A>
344 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
345 smoothers_(amg.smoothers_), solver_(amg.solver_),
346 rhs_(), lhs_(), update_(),
347 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
348 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
349 buildHierarchy_(amg.buildHierarchy_),
350 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
351 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
361 template<
class M,
class X,
class S,
class PI,
class A>
364 std::size_t gamma, std::size_t preSmoothingSteps,
365 std::size_t postSmoothingSteps,
bool additive_)
366 : matrices_(&matrices), smootherArgs_(smootherArgs),
368 rhs_(), lhs_(), update_(), scalarProduct_(0),
369 gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
370 additive(additive_), coarsesolverconverged(true),
371 coarseSmoother_(), verbosity_(2)
373 assert(matrices_->isBuilt());
376 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
379 template<
class M,
class X,
class S,
class PI,
class A>
383 : matrices_(&matrices), smootherArgs_(smootherArgs),
385 rhs_(), lhs_(), update_(), scalarProduct_(0),
386 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
387 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
388 additive(parms.getAdditive()), coarsesolverconverged(true),
389 coarseSmoother_(), verbosity_(parms.debugLevel())
391 assert(matrices_->isBuilt());
394 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
397 template<
class M,
class X,
class S,
class PI,
class A>
402 std::size_t gamma, std::size_t preSmoothingSteps,
403 std::size_t postSmoothingSteps,
406 : smootherArgs_(smootherArgs),
408 rhs_(), lhs_(), update_(), scalarProduct_(), gamma_(gamma),
409 preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
410 additive(additive_), coarsesolverconverged(true),
411 coarseSmoother_(), verbosity_(criterion.debugLevel())
413 static_assert(
static_cast<int>(M::category)==
static_cast<int>(S::category),
414 "Matrix and Solver must match in terms of category!");
418 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
421 template<
class M,
class X,
class S,
class PI,
class A>
427 : smootherArgs_(smootherArgs),
429 rhs_(), lhs_(), update_(), scalarProduct_(),
430 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
431 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
432 additive(criterion.getAdditive()), coarsesolverconverged(true),
433 coarseSmoother_(), verbosity_(criterion.debugLevel())
435 static_assert(
static_cast<int>(M::category)==
static_cast<int>(S::category),
436 "Matrix and Solver must match in terms of category!");
440 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
444 template<
class M,
class X,
class S,
class PI,
class A>
447 if(buildHierarchy_) {
451 coarseSmoother_.reset();
464 template<
class M,
class X,
class S,
class PI,
class A>
466 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
470 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
472 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
475 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
477 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
479 SmootherArgs sargs(smootherArgs_);
480 sargs.iterations = 1;
483 cargs.setArgs(sargs);
484 if(matrices_->redistributeInformation().back().isSetup()) {
486 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
487 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
489 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
490 cargs.setComm(*matrices_->parallelInformation().coarsest());
494 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
496#if HAVE_SUPERLU || HAVE_UMFPACK
498#define DIRECTSOLVER UMFPack
500#define DIRECTSOLVER SuperLU
503 if(is_same<ParallelInformation,SequentialInformation>::value
504 || matrices_->parallelInformation().coarsest()->communicator().size()==1
505 || (matrices_->parallelInformation().coarsest().isRedistributed()
506 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
507 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
508 if(matrices_->parallelInformation().coarsest().isRedistributed())
510 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
512 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
516 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
517 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
518 std::cout<<
"Using a direct coarse solver (" <<
static_cast< DIRECTSOLVER<typename M::matrix_type>*
>(solver_.get())->name() <<
")" << std::endl;
523 if(matrices_->parallelInformation().coarsest().isRedistributed())
525 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
527 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
529 *coarseSmoother_, 1E-2, 1000, 0));
533 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
535 *coarseSmoother_, 1E-2, 1000, 0));
539 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
540 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
541 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
545 template<
class M,
class X,
class S,
class PI,
class A>
552 typedef typename M::matrix_type
Matrix;
559 const Matrix& mat=matrices_->matrices().finest()->getmat();
560 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
561 bool isDirichlet =
true;
562 bool hasDiagonal =
false;
564 for(ColIter col=row->begin(); col!=row->end(); ++col) {
565 if(row.index()==col.index()) {
573 if(isDirichlet && hasDiagonal)
574 diagonal.solve(x[row.index()], b[row.index()]);
577 if(smoothers_->levels()>0)
578 smoothers_->finest()->pre(x,b);
581 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
594 matrices_->coarsenVector(*rhs_);
595 matrices_->coarsenVector(*lhs_);
596 matrices_->coarsenVector(*update_);
602 Iterator coarsest = smoothers_->
coarsest();
603 Iterator smoother = smoothers_->finest();
604 RIterator rhs = rhs_->finest();
605 DIterator lhs = lhs_->finest();
606 if(smoothers_->levels()>0) {
608 assert(lhs_->levels()==rhs_->levels());
609 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
610 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
612 if(smoother!=coarsest)
613 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
614 smoother->pre(*lhs,*rhs);
615 smoother->pre(*lhs,*rhs);
625 template<
class M,
class X,
class S,
class PI,
class A>
628 return matrices_->levels();
630 template<
class M,
class X,
class S,
class PI,
class A>
631 std::size_t AMG<M,X,S,PI,A>::maxlevels()
633 return matrices_->maxlevels();
637 template<
class M,
class X,
class S,
class PI,
class A>
640 LevelContext levelContext;
648 initIteratorsWithFineLevel(levelContext);
651 *levelContext.lhs = v;
652 *levelContext.rhs = d;
653 *levelContext.update=0;
654 levelContext.level=0;
658 if(postSteps_==0||matrices_->maxlevels()==1)
659 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
661 v=*levelContext.update;
666 template<
class M,
class X,
class S,
class PI,
class A>
669 levelContext.smoother = smoothers_->finest();
670 levelContext.matrix = matrices_->matrices().finest();
671 levelContext.pinfo = matrices_->parallelInformation().finest();
672 levelContext.redist =
673 matrices_->redistributeInformation().begin();
674 levelContext.aggregates = matrices_->aggregatesMaps().begin();
675 levelContext.lhs = lhs_->finest();
676 levelContext.update = update_->finest();
677 levelContext.rhs = rhs_->finest();
680 template<
class M,
class X,
class S,
class PI,
class A>
682 ::moveToCoarseLevel(LevelContext& levelContext)
685 bool processNextLevel=
true;
687 if(levelContext.redist->isSetup()) {
688 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
689 levelContext.rhs.getRedistributed());
690 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
691 if(processNextLevel) {
694 ++levelContext.pinfo;
695 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
696 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
697 static_cast<const Range&
>(fineRhs.getRedistributed()),
698 *levelContext.pinfo);
703 ++levelContext.pinfo;
704 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
705 ::restrictVector(*(*levelContext.aggregates),
706 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
707 *levelContext.pinfo);
710 if(processNextLevel) {
713 ++levelContext.update;
714 ++levelContext.matrix;
715 ++levelContext.level;
716 ++levelContext.redist;
718 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
720 ++levelContext.smoother;
721 ++levelContext.aggregates;
724 *levelContext.update=0;
726 return processNextLevel;
729 template<
class M,
class X,
class S,
class PI,
class A>
731 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
733 if(processNextLevel) {
734 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
736 --levelContext.smoother;
737 --levelContext.aggregates;
739 --levelContext.redist;
740 --levelContext.level;
742 --levelContext.matrix;
746 --levelContext.pinfo;
748 if(levelContext.redist->isSetup()) {
750 levelContext.lhs.getRedistributed()=0;
751 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
752 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
753 levelContext.lhs.getRedistributed(),
754 matrices_->getProlongationDampingFactor(),
755 *levelContext.pinfo, *levelContext.redist);
758 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
759 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
760 matrices_->getProlongationDampingFactor(),
761 *levelContext.pinfo);
765 if(processNextLevel) {
766 --levelContext.update;
770 *levelContext.update += *levelContext.lhs;
773 template<
class M,
class X,
class S,
class PI,
class A>
776 return IsDirectSolver< CoarseSolver>::value;
779 template<
class M,
class X,
class S,
class PI,
class A>
781 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
785 if(levelContext.redist->isSetup()) {
786 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
787 if(levelContext.rhs.getRedistributed().size()>0) {
789 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
790 levelContext.rhs.getRedistributed());
791 solver_->apply(levelContext.update.getRedistributed(),
792 levelContext.rhs.getRedistributed(), res);
794 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
795 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
797 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
798 solver_->apply(*levelContext.update, *levelContext.rhs, res);
802 coarsesolverconverged =
false;
807#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
808 bool processNextLevel = moveToCoarseLevel(levelContext);
810 if(processNextLevel) {
812 for(std::size_t i=0; i<gamma_; i++)
816 moveToFineLevel(levelContext, processNextLevel);
821 if(levelContext.matrix == matrices_->matrices().finest()) {
822 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
823 if(!coarsesolverconverged)
824 DUNE_THROW(MathError,
"Coarse solver did not converge");
832 template<
class M,
class X,
class S,
class PI,
class A>
833 void AMG<M,X,S,PI,A>::additiveMgc(){
836 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
839 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
843 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
844 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
850 lhs = lhs_->finest();
853 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
856 smoother->apply(*lhs, *rhs);
860#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
861 InverseOperatorResult res;
862 pinfo->copyOwnerToAll(*rhs, *rhs);
863 solver_->apply(*lhs, *rhs, res);
866 DUNE_THROW(MathError,
"Coarse solver did not converge");
875 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
876 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
882 template<
class M,
class X,
class S,
class PI,
class A>
889 Iterator coarsest = smoothers_->
coarsest();
890 Iterator smoother = smoothers_->finest();
891 DIterator lhs = lhs_->finest();
892 if(smoothers_->levels()>0) {
893 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
894 smoother->post(*lhs);
895 if(smoother!=coarsest)
896 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
897 smoother->post(*lhs);
898 smoother->post(*lhs);
911 template<
class M,
class X,
class S,
class PI,
class A>
915 matrices_->getCoarsestAggregatesOnFinest(cont);
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:56
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:257
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:260
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:141
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:31
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:317
All parameters for AMG.
Definition: parameters.hh:391
A generic dynamic dense matrix.
Definition: matrix.hh:25
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
T block_type
Export the type representing the components.
Definition: matrix.hh:32
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:26
ConstIterator class for sequential access.
Definition: vbvector.hh:647
A few common exception classes.
#define DUNE_DEPRECATED
Mark some entity as deprecated.
Definition: deprecated.hh:84
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:546
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:267
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:271
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:774
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, bool additive=false) DUNE_DEPRECATED
Construct a new amg with a specific coarse solver.
Definition: amg.hh:362
X Domain
The domain type.
Definition: amg.hh:78
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:251
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:259
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:91
S Smoother
The type of the smoother.
Definition: amg.hh:88
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:243
M Operator
The matrix operator type.
Definition: amg.hh:64
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:247
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:255
X Range
The range type.
Definition: amg.hh:80
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:408
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:913
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:263
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:210
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1385
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:75
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:82
void post(Domain &x)
Clean up.
Definition: amg.hh:883
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:430
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:275
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:638
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:73
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:71
@ category
The solver category.
Definition: amg.hh:95
Provides a classes representing the hierarchies in AMG.
Dune namespace.
Definition: alignment.hh:10
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:32
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Choose the approriate scalar product for a solver category.
Definition: scalarproducts.hh:77
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.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18