3#ifndef DUNE_AMG_AMG_CPR_HH
4#define DUNE_AMG_AMG_CPR_HH
9#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
11#include <dune/common/exceptions.hh>
12#include <dune/common/version.hh>
13#include <dune/istl/paamg/amg.hh>
14#include <dune/istl/paamg/smoother.hh>
15#include <dune/istl/paamg/transfer.hh>
16#include <dune/istl/paamg/hierarchy.hh>
17#include <dune/istl/solvers.hh>
18#include <dune/istl/scalarproducts.hh>
19#include <dune/istl/superlu.hh>
20#include <dune/istl/umfpack.hh>
21#include <dune/istl/solvertype.hh>
22#include <dune/common/typetraits.hh>
23#include <dune/common/exceptions.hh>
34 template<
class M,
class T>
35 void redistributeMatrixAmg(M&, M&, SequentialInformation&, SequentialInformation&, T&)
40 template<
class M,
class PI>
41 typename std::enable_if<!std::is_same<PI,SequentialInformation>::value,
void>::type
42 redistributeMatrixAmg(M& mat, M& matRedist, PI& info, PI& infoRedist,
43 Dune::RedistributeInformation<PI>& redistInfo)
45 info.buildGlobalLookup(mat.N());
46 redistributeMatrixEntries(mat, matRedist, info, infoRedist, redistInfo);
47 info.freeGlobalLookup();
68 template<
class M,
class X,
class S,
class P,
class K,
class A>
84 template<
class M,
class X,
class S,
class PI=SequentialInformation,
85 class A=std::allocator<X> >
88 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
134 const SmootherArgs& smootherArgs,
const Parameters& parms);
181 std::size_t levels();
183 std::size_t maxlevels();
195 auto copyFlags = NegateSet<typename PI::OwnerSet>();
196 const auto& matrices = matrices_->matrices();
197 const auto& aggregatesMapHierarchy = matrices_->aggregatesMaps();
198 const auto& infoHierarchy = matrices_->parallelInformation();
199 const auto& redistInfoHierarchy = matrices_->redistributeInformation();
200 BaseGalerkinProduct productBuilder;
201 auto aggregatesMap = aggregatesMapHierarchy.begin();
202 auto info = infoHierarchy.finest();
203 auto redistInfo = redistInfoHierarchy.begin();
204 auto matrix = matrices.finest();
205 auto coarsestMatrix = matrices.coarsest();
207 using Matrix =
typename M::matrix_type;
210 if(matrix.isRedistributed()) {
211 redistributeMatrixAmg(
const_cast<Matrix&
>(matrix->getmat()),
212 const_cast<Matrix&
>(matrix.getRedistributed().getmat()),
213 const_cast<PI&
>(*info),
const_cast<PI&
>(info.getRedistributed()),
214 const_cast<Dune::RedistributeInformation<PI>&
>(*redistInfo));
218 for(; matrix!=coarsestMatrix; ++aggregatesMap) {
219 const Matrix& fine = (matrix.isRedistributed() ? matrix.getRedistributed() : *matrix).getmat();
223 productBuilder.calculate(fine, *(*aggregatesMap),
const_cast<Matrix&
>(matrix->getmat()), *info, copyFlags);
225 if(matrix.isRedistributed()) {
226 redistributeMatrixAmg(
const_cast<Matrix&
>(matrix->getmat()),
227 const_cast<Matrix&
>(matrix.getRedistributed().getmat()),
228 const_cast<PI&
>(*info),
const_cast<PI&
>(info.getRedistributed()),
229 const_cast<Dune::RedistributeInformation<PI>&
>(*redistInfo));
260 void createHierarchies(C& criterion,
Operator& matrix,
264 std::shared_ptr< Operator > op( &matrix, [](
Operator* ) {});
265 std::shared_ptr< PI > pifo(
const_cast< PI*
> (&pinfo), []( PI * ) {});
266 createHierarchies( criterion, op, pifo);
270 void createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
271 std::shared_ptr< PI > pinfo );
273 void setupCoarseSolver();
291 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
295 typename ParallelInformationHierarchy::Iterator
pinfo;
299 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
303 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
307 typename Hierarchy<Domain,A>::Iterator
lhs;
311 typename Hierarchy<Domain,A>::Iterator
update;
315 typename Hierarchy<Range,A>::Iterator
rhs;
327 void mgc(LevelContext& levelContext);
337 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
343 bool moveToCoarseLevel(LevelContext& levelContext);
349 void initIteratorsWithFineLevel(LevelContext& levelContext);
352 std::shared_ptr<OperatorHierarchy> matrices_;
356 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
358 std::shared_ptr<CoarseSolver> solver_;
360 std::shared_ptr< Hierarchy<Range,A> > rhs_;
362 std::shared_ptr< Hierarchy<Domain,A> > lhs_;
364 std::shared_ptr< Hierarchy<Domain,A> > update_;
366 using ScalarProduct = Dune::ScalarProduct<X>;
368 std::shared_ptr<ScalarProduct> scalarProduct_;
372 std::size_t preSteps_;
374 std::size_t postSteps_;
375 bool buildHierarchy_;
377 bool coarsesolverconverged;
378 std::shared_ptr<Smoother> coarseSmoother_;
380 SolverCategory::Category category_;
382 std::size_t verbosity_;
385 template<
class M,
class X,
class S,
class PI,
class A>
387 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
388 smoothers_(amg.smoothers_), solver_(amg.solver_),
389 rhs_(), lhs_(), update_(),
390 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
391 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
392 buildHierarchy_(amg.buildHierarchy_),
393 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
394 coarseSmoother_(amg.coarseSmoother_),
395 category_(amg.category_),
396 verbosity_(amg.verbosity_)
399 rhs_.reset(
new Hierarchy<Range,A>(*amg.rhs_) );
401 lhs_.reset(
new Hierarchy<Domain,A>(*amg.lhs_) );
403 update_.reset(
new Hierarchy<Domain,A>(*amg.update_) );
406 template<
class M,
class X,
class S,
class PI,
class A>
409 const Parameters& parms)
410 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
411 smoothers_(new Hierarchy<
Smoother,A>), solver_(&coarseSolver),
412 rhs_(), lhs_(), update_(), scalarProduct_(0),
413 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
414 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
415 additive(parms.getAdditive()), coarsesolverconverged(true),
418 category_(SolverCategory::category(*smoothers_->coarsest())),
419 verbosity_(parms.debugLevel())
421 assert(matrices_->isBuilt());
424 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
427 template<
class M,
class X,
class S,
class PI,
class A>
433 : smootherArgs_(smootherArgs),
434 smoothers_(new Hierarchy<
Smoother,A>), solver_(),
435 rhs_(), lhs_(), update_(), scalarProduct_(),
436 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
437 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
438 additive(criterion.getAdditive()), coarsesolverconverged(true),
440 category_(SolverCategory::category(pinfo)),
441 verbosity_(criterion.debugLevel())
443 if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
444 DUNE_THROW(InvalidSolverCategory,
"Matrix and Communication must have the same SolverCategory!");
445 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
448 template<
class M,
class X,
class S,
class PI,
class A>
451 if(buildHierarchy_) {
455 coarseSmoother_.reset();
459 template<
class M,
class X,
class S,
class PI,
class A>
466 template<
class M,
class X,
class S,
class PI,
class A>
470 smoothers_.reset(
new Hierarchy<Smoother,A>);
472 coarseSmoother_.reset();
473 scalarProduct_.reset();
474 buildHierarchy_=
true;
475 coarsesolverconverged =
true;
476 smoothers_.reset(
new Hierarchy<Smoother,A>);
477 recalculateHierarchy();
478 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
480 if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
481 std::cout <<
"Recalculating galerkin and coarse smoothers "<< matrices_->maxlevels() <<
" levels "
482 << watch.elapsed() <<
" seconds." << std::endl;
486 template<
class M,
class X,
class S,
class PI,
class A>
489 std::shared_ptr< PI > pinfo )
492 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
494 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
497 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
499 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
500 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
501 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
504 template<
class M,
class X,
class S,
class PI,
class A>
505 void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
511 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
512 && ( ! matrices_->redistributeInformation().back().isSetup() ||
513 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
516 SmootherArgs sargs(smootherArgs_);
517 sargs.iterations = 1;
519 typename ConstructionTraits<Smoother>::Arguments cargs;
520 cargs.setArgs(sargs);
521 if(matrices_->redistributeInformation().back().isSetup()) {
523 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
524 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
526 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
527 cargs.setComm(*matrices_->parallelInformation().coarsest());
530 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
531 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
533 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
536 if( SolverSelector::isDirectSolver &&
537 (std::is_same<ParallelInformation,SequentialInformation>::value
538 || matrices_->parallelInformation().coarsest()->communicator().size()==1
539 || (matrices_->parallelInformation().coarsest().isRedistributed()
540 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
541 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
544 if(matrices_->parallelInformation().coarsest().isRedistributed())
546 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
549 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
556 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
558 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
559 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
563 if(matrices_->parallelInformation().coarsest().isRedistributed())
565 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
570 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
572 reinterpret_cast<typename
573 std::conditional<std::is_same<PI,SequentialInformation>::value,
574 Dune::SeqScalarProduct<X>,
575 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
576 *coarseSmoother_, 1E-2, 1000, 0));
581 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
583 reinterpret_cast<typename
584 std::conditional<std::is_same<PI,SequentialInformation>::value,
585 Dune::SeqScalarProduct<X>,
586 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
587 *coarseSmoother_, 1E-2, 1000, 0));
607 template<
class M,
class X,
class S,
class PI,
class A>
614 typedef typename M::matrix_type Matrix;
615 typedef typename Matrix::ConstRowIterator RowIter;
616 typedef typename Matrix::ConstColIterator ColIter;
617 typedef typename Matrix::block_type Block;
619 zero=
typename Matrix::field_type();
621 const Matrix& mat=matrices_->matrices().finest()->getmat();
622 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
623 bool isDirichlet =
true;
624 bool hasDiagonal =
false;
626 for(ColIter col=row->begin(); col!=row->end(); ++col) {
627 if(row.index()==col.index()) {
635 if(isDirichlet && hasDiagonal)
636 diagonal.solve(x[row.index()], b[row.index()]);
639 if(smoothers_->levels()>0)
640 smoothers_->finest()->pre(x,b);
643 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
646 typedef std::shared_ptr< Range > RangePtr ;
647 typedef std::shared_ptr< Domain > DomainPtr;
650 RangePtr copy(
new Range(b) );
651 rhs_.reset(
new Hierarchy<Range,A>(copy) );
652 DomainPtr dcopy(
new Domain(x) );
653 lhs_.reset(
new Hierarchy<Domain,A>(dcopy) );
654 DomainPtr dcopy2(
new Domain(x) );
655 update_.reset(
new Hierarchy<Domain,A>(dcopy2) );
657 matrices_->coarsenVector(*rhs_);
658 matrices_->coarsenVector(*lhs_);
659 matrices_->coarsenVector(*update_);
662 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
663 typedef typename Hierarchy<Range,A>::Iterator RIterator;
664 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
665 Iterator coarsest = smoothers_->coarsest();
666 Iterator smoother = smoothers_->finest();
667 RIterator rhs = rhs_->finest();
668 DIterator lhs = lhs_->finest();
669 if(smoothers_->levels()>0) {
671 assert(lhs_->levels()==rhs_->levels());
672 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
673 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
675 if(smoother!=coarsest)
676 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
677 smoother->pre(*lhs,*rhs);
678 smoother->pre(*lhs,*rhs);
688 template<
class M,
class X,
class S,
class PI,
class A>
691 return matrices_->levels();
693 template<
class M,
class X,
class S,
class PI,
class A>
694 std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
696 return matrices_->maxlevels();
700 template<
class M,
class X,
class S,
class PI,
class A>
703 LevelContext levelContext;
711 initIteratorsWithFineLevel(levelContext);
714 *levelContext.lhs = v;
715 *levelContext.rhs = d;
716 *levelContext.update=0;
717 levelContext.level=0;
721 if(postSteps_==0||matrices_->maxlevels()==1)
722 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
724 v=*levelContext.update;
729 template<
class M,
class X,
class S,
class PI,
class A>
732 levelContext.smoother = smoothers_->finest();
733 levelContext.matrix = matrices_->matrices().finest();
734 levelContext.pinfo = matrices_->parallelInformation().finest();
735 levelContext.redist =
736 matrices_->redistributeInformation().begin();
737 levelContext.aggregates = matrices_->aggregatesMaps().begin();
738 levelContext.lhs = lhs_->finest();
739 levelContext.update = update_->finest();
740 levelContext.rhs = rhs_->finest();
743 template<
class M,
class X,
class S,
class PI,
class A>
744 bool AMGCPR<M,X,S,PI,A>
745 ::moveToCoarseLevel(LevelContext& levelContext)
748 bool processNextLevel=
true;
750 if(levelContext.redist->isSetup()) {
751 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
752 levelContext.rhs.getRedistributed());
753 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
754 if(processNextLevel) {
756 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
757 ++levelContext.pinfo;
758 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
759 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
760 static_cast<const Range&
>(fineRhs.getRedistributed()),
761 *levelContext.pinfo);
765 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
766 ++levelContext.pinfo;
767 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
768 ::restrictVector(*(*levelContext.aggregates),
769 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
770 *levelContext.pinfo);
773 if(processNextLevel) {
776 ++levelContext.update;
777 ++levelContext.matrix;
778 ++levelContext.level;
779 ++levelContext.redist;
781 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
783 ++levelContext.smoother;
784 ++levelContext.aggregates;
787 *levelContext.update=0;
789 return processNextLevel;
792 template<
class M,
class X,
class S,
class PI,
class A>
793 void AMGCPR<M,X,S,PI,A>
794 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
796 if(processNextLevel) {
797 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
799 --levelContext.smoother;
800 --levelContext.aggregates;
802 --levelContext.redist;
803 --levelContext.level;
805 --levelContext.matrix;
809 --levelContext.pinfo;
811 if(levelContext.redist->isSetup()) {
813 levelContext.lhs.getRedistributed()=0;
814 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
815 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
816 levelContext.lhs.getRedistributed(),
817 matrices_->getProlongationDampingFactor(),
818 *levelContext.pinfo, *levelContext.redist);
821 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
822 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
823 matrices_->getProlongationDampingFactor(),
824 *levelContext.pinfo);
828 if(processNextLevel) {
829 --levelContext.update;
833 *levelContext.update += *levelContext.lhs;
836 template<
class M,
class X,
class S,
class PI,
class A>
839 return IsDirectSolver< CoarseSolver>::value;
842 template<
class M,
class X,
class S,
class PI,
class A>
844 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
846 InverseOperatorResult res;
848 if(levelContext.redist->isSetup()) {
849 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
850 if(levelContext.rhs.getRedistributed().size()>0) {
852 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
853 levelContext.rhs.getRedistributed());
854 solver_->apply(levelContext.update.getRedistributed(),
855 levelContext.rhs.getRedistributed(), res);
857 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
858 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
860 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
861 solver_->apply(*levelContext.update, *levelContext.rhs, res);
865 coarsesolverconverged =
false;
868 presmooth(levelContext, preSteps_);
870#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
871 bool processNextLevel = moveToCoarseLevel(levelContext);
873 if(processNextLevel) {
875 for(std::size_t i=0; i<gamma_; i++)
879 moveToFineLevel(levelContext, processNextLevel);
884 if(levelContext.matrix == matrices_->matrices().finest()) {
885 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
886 if(!coarsesolverconverged){
891 postsmooth(levelContext, postSteps_);
896 template<
class M,
class X,
class S,
class PI,
class A>
897 void AMGCPR<M,X,S,PI,A>::additiveMgc(){
900 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
901 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
902 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
903 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
905 for(
typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
907 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
908 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
914 lhs = lhs_->finest();
915 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
917 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
920 smoother->apply(*lhs, *rhs);
924#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
925 InverseOperatorResult res;
926 pinfo->copyOwnerToAll(*rhs, *rhs);
927 solver_->apply(*lhs, *rhs, res);
930 DUNE_THROW(MathError,
"Coarse solver did not converge");
938 for(
typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
939 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
940 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
946 template<
class M,
class X,
class S,
class PI,
class A>
950 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
951 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
952 Iterator coarsest = smoothers_->coarsest();
953 Iterator smoother = smoothers_->finest();
954 DIterator lhs = lhs_->finest();
955 if(smoothers_->levels()>0) {
956 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
957 smoother->post(*lhs);
958 if(smoother!=coarsest)
959 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
960 smoother->post(*lhs);
961 smoother->post(*lhs);
971 template<
class M,
class X,
class S,
class PI,
class A>
975 matrices_->getCoarsestAggregatesOnFinest(cont);
Parallel algebraic multigrid based on agglomeration.
Definition amgcpr.hh:87
Interface class adding the update() method to the preconditioner interface.
Definition PreconditionerWithUpdate.hpp:32
M Operator
The matrix operator type.
Definition amgcpr.hh:95
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition amgcpr.hh:287
AMGCPR(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition amgcpr.hh:407
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition amgcpr.hh:311
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition amgcpr.hh:315
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition amgcpr.hh:104
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition amgcpr.hh:295
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition amgcpr.hh:166
void apply(Domain &v, const Range &d)
Definition amgcpr.hh:701
void pre(Domain &x, Range &b)
Definition amgcpr.hh:608
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition amgcpr.hh:291
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition amgcpr.hh:973
void post(Domain &x)
Definition amgcpr.hh:947
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition amgcpr.hh:106
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition amgcpr.hh:193
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition amgcpr.hh:122
X Range
The range type.
Definition amgcpr.hh:111
X Domain
The domain type.
Definition amgcpr.hh:109
S Smoother
The type of the smoother.
Definition amgcpr.hh:119
void updateSolver(C &criterion, const Operator &, const PI &pinfo)
Update the coarse solver and the hierarchies.
Definition amgcpr.hh:461
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition amgcpr.hh:303
virtual void update()
Update the coarse solver and the hierarchies.
Definition amgcpr.hh:467
PI ParallelInformation
The type of the parallel information.
Definition amgcpr.hh:102
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition amgcpr.hh:113
std::size_t level
The level index.
Definition amgcpr.hh:319
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition amgcpr.hh:307
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition amgcpr.hh:837
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition amgcpr.hh:299