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));
259 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7 )
261 void createHierarchies(C& criterion,
Operator& matrix,
265 std::shared_ptr< Operator > op( &matrix, [](
Operator* ) {});
266 std::shared_ptr< PI > pifo(
const_cast< PI*
> (&pinfo), []( PI * ) {});
267 createHierarchies( criterion, op, pifo);
271 void createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
272 std::shared_ptr< PI > pinfo );
276 void createHierarchies(C& criterion,
Operator& matrix,
280 void setupCoarseSolver();
298 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
302 typename ParallelInformationHierarchy::Iterator
pinfo;
306 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
310 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
314 typename Hierarchy<Domain,A>::Iterator
lhs;
318 typename Hierarchy<Domain,A>::Iterator
update;
322 typename Hierarchy<Range,A>::Iterator
rhs;
334 void mgc(LevelContext& levelContext);
344 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
350 bool moveToCoarseLevel(LevelContext& levelContext);
356 void initIteratorsWithFineLevel(LevelContext& levelContext);
359 std::shared_ptr<OperatorHierarchy> matrices_;
363 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
365 std::shared_ptr<CoarseSolver> solver_;
367 std::shared_ptr< Hierarchy<Range,A> > rhs_;
369 std::shared_ptr< Hierarchy<Domain,A> > lhs_;
371 std::shared_ptr< Hierarchy<Domain,A> > update_;
373 using ScalarProduct = Dune::ScalarProduct<X>;
375 std::shared_ptr<ScalarProduct> scalarProduct_;
379 std::size_t preSteps_;
381 std::size_t postSteps_;
382 bool buildHierarchy_;
384 bool coarsesolverconverged;
385 std::shared_ptr<Smoother> coarseSmoother_;
387 SolverCategory::Category category_;
389 std::size_t verbosity_;
392 template<
class M,
class X,
class S,
class PI,
class A>
394 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
395 smoothers_(amg.smoothers_), solver_(amg.solver_),
396 rhs_(), lhs_(), update_(),
397 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
398 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
399 buildHierarchy_(amg.buildHierarchy_),
400 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
401 coarseSmoother_(amg.coarseSmoother_),
402 category_(amg.category_),
403 verbosity_(amg.verbosity_)
406 rhs_.reset(
new Hierarchy<Range,A>(*amg.rhs_) );
408 lhs_.reset(
new Hierarchy<Domain,A>(*amg.lhs_) );
410 update_.reset(
new Hierarchy<Domain,A>(*amg.update_) );
413 template<
class M,
class X,
class S,
class PI,
class A>
416 const Parameters& parms)
417 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
418 smoothers_(new Hierarchy<
Smoother,A>), solver_(&coarseSolver),
419 rhs_(), lhs_(), update_(), scalarProduct_(0),
420 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
421 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
422 additive(parms.getAdditive()), coarsesolverconverged(true),
425 category_(SolverCategory::category(*smoothers_->coarsest())),
426 verbosity_(parms.debugLevel())
428 assert(matrices_->isBuilt());
431 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
434 template<
class M,
class X,
class S,
class PI,
class A>
440 : smootherArgs_(smootherArgs),
441 smoothers_(new Hierarchy<
Smoother,A>), solver_(),
442 rhs_(), lhs_(), update_(), scalarProduct_(),
443 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
444 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
445 additive(criterion.getAdditive()), coarsesolverconverged(true),
447 category_(SolverCategory::category(pinfo)),
448 verbosity_(criterion.debugLevel())
450 if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
451 DUNE_THROW(InvalidSolverCategory,
"Matrix and Communication must have the same SolverCategory!");
452 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
455 template<
class M,
class X,
class S,
class PI,
class A>
458 if(buildHierarchy_) {
462 coarseSmoother_.reset();
466 template<
class M,
class X,
class S,
class PI,
class A>
473 template<
class M,
class X,
class S,
class PI,
class A>
477 smoothers_.reset(
new Hierarchy<Smoother,A>);
479 coarseSmoother_.reset();
480 scalarProduct_.reset();
481 buildHierarchy_=
true;
482 coarsesolverconverged =
true;
483 smoothers_.reset(
new Hierarchy<Smoother,A>);
484 recalculateHierarchy();
485 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
487 if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
488 std::cout <<
"Recalculating galerkin and coarse smoothers "<< matrices_->maxlevels() <<
" levels "
489 << watch.elapsed() <<
" seconds." << std::endl;
493 template<
class M,
class X,
class S,
class PI,
class A>
495 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
497 std::shared_ptr< PI > pinfo )
503 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
505 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
508 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
510 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
511 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
512 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
515 template<
class M,
class X,
class S,
class PI,
class A>
516 void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
522 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
523 && ( ! matrices_->redistributeInformation().back().isSetup() ||
524 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
527 SmootherArgs sargs(smootherArgs_);
528 sargs.iterations = 1;
530 typename ConstructionTraits<Smoother>::Arguments cargs;
531 cargs.setArgs(sargs);
532 if(matrices_->redistributeInformation().back().isSetup()) {
534 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
535 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
537 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
538 cargs.setComm(*matrices_->parallelInformation().coarsest());
541 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
542 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
544 coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
547 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
550 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
553 if( SolverSelector::isDirectSolver &&
554 (std::is_same<ParallelInformation,SequentialInformation>::value
555 || matrices_->parallelInformation().coarsest()->communicator().size()==1
556 || (matrices_->parallelInformation().coarsest().isRedistributed()
557 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
558 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
561 if(matrices_->parallelInformation().coarsest().isRedistributed())
563 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
566 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
573 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
575 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
576 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
580 if(matrices_->parallelInformation().coarsest().isRedistributed())
582 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
587 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
589 reinterpret_cast<typename
590 std::conditional<std::is_same<PI,SequentialInformation>::value,
591 Dune::SeqScalarProduct<X>,
592 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
593 *coarseSmoother_, 1E-2, 1000, 0));
598 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
600 reinterpret_cast<typename
601 std::conditional<std::is_same<PI,SequentialInformation>::value,
602 Dune::SeqScalarProduct<X>,
603 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
604 *coarseSmoother_, 1E-2, 1000, 0));
624 template<
class M,
class X,
class S,
class PI,
class A>
631 typedef typename M::matrix_type Matrix;
632 typedef typename Matrix::ConstRowIterator RowIter;
633 typedef typename Matrix::ConstColIterator ColIter;
634 typedef typename Matrix::block_type Block;
636 zero=
typename Matrix::field_type();
638 const Matrix& mat=matrices_->matrices().finest()->getmat();
639 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
640 bool isDirichlet =
true;
641 bool hasDiagonal =
false;
643 for(ColIter col=row->begin(); col!=row->end(); ++col) {
644 if(row.index()==col.index()) {
652 if(isDirichlet && hasDiagonal)
653 diagonal.solve(x[row.index()], b[row.index()]);
656 if(smoothers_->levels()>0)
657 smoothers_->finest()->pre(x,b);
660 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
663 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
664 typedef std::shared_ptr< Range > RangePtr ;
665 typedef std::shared_ptr< Domain > DomainPtr;
667 typedef Range* RangePtr;
668 typedef Domain* DomainPtr;
672 RangePtr copy(
new Range(b) );
673 rhs_.reset(
new Hierarchy<Range,A>(copy) );
674 DomainPtr dcopy(
new Domain(x) );
675 lhs_.reset(
new Hierarchy<Domain,A>(dcopy) );
676 DomainPtr dcopy2(
new Domain(x) );
677 update_.reset(
new Hierarchy<Domain,A>(dcopy2) );
679 matrices_->coarsenVector(*rhs_);
680 matrices_->coarsenVector(*lhs_);
681 matrices_->coarsenVector(*update_);
684 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
685 typedef typename Hierarchy<Range,A>::Iterator RIterator;
686 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
687 Iterator coarsest = smoothers_->coarsest();
688 Iterator smoother = smoothers_->finest();
689 RIterator rhs = rhs_->finest();
690 DIterator lhs = lhs_->finest();
691 if(smoothers_->levels()>0) {
693 assert(lhs_->levels()==rhs_->levels());
694 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
695 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
697 if(smoother!=coarsest)
698 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
699 smoother->pre(*lhs,*rhs);
700 smoother->pre(*lhs,*rhs);
710 template<
class M,
class X,
class S,
class PI,
class A>
713 return matrices_->levels();
715 template<
class M,
class X,
class S,
class PI,
class A>
716 std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
718 return matrices_->maxlevels();
722 template<
class M,
class X,
class S,
class PI,
class A>
725 LevelContext levelContext;
733 initIteratorsWithFineLevel(levelContext);
736 *levelContext.lhs = v;
737 *levelContext.rhs = d;
738 *levelContext.update=0;
739 levelContext.level=0;
743 if(postSteps_==0||matrices_->maxlevels()==1)
744 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
746 v=*levelContext.update;
751 template<
class M,
class X,
class S,
class PI,
class A>
754 levelContext.smoother = smoothers_->finest();
755 levelContext.matrix = matrices_->matrices().finest();
756 levelContext.pinfo = matrices_->parallelInformation().finest();
757 levelContext.redist =
758 matrices_->redistributeInformation().begin();
759 levelContext.aggregates = matrices_->aggregatesMaps().begin();
760 levelContext.lhs = lhs_->finest();
761 levelContext.
update = update_->finest();
762 levelContext.rhs = rhs_->finest();
765 template<
class M,
class X,
class S,
class PI,
class A>
766 bool AMGCPR<M,X,S,PI,A>
767 ::moveToCoarseLevel(LevelContext& levelContext)
770 bool processNextLevel=
true;
772 if(levelContext.redist->isSetup()) {
773 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
774 levelContext.rhs.getRedistributed());
775 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
776 if(processNextLevel) {
778 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
779 ++levelContext.pinfo;
780 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
781 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
782 static_cast<const Range&
>(fineRhs.getRedistributed()),
783 *levelContext.pinfo);
787 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
788 ++levelContext.pinfo;
789 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
790 ::restrictVector(*(*levelContext.aggregates),
791 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
792 *levelContext.pinfo);
795 if(processNextLevel) {
798 ++levelContext.update;
799 ++levelContext.matrix;
800 ++levelContext.level;
801 ++levelContext.redist;
803 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
805 ++levelContext.smoother;
806 ++levelContext.aggregates;
809 *levelContext.update=0;
811 return processNextLevel;
814 template<
class M,
class X,
class S,
class PI,
class A>
815 void AMGCPR<M,X,S,PI,A>
816 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
818 if(processNextLevel) {
819 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
821 --levelContext.smoother;
822 --levelContext.aggregates;
824 --levelContext.redist;
825 --levelContext.level;
827 --levelContext.matrix;
831 --levelContext.pinfo;
833 if(levelContext.redist->isSetup()) {
835 levelContext.lhs.getRedistributed()=0;
836 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
837 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
838 levelContext.lhs.getRedistributed(),
839 matrices_->getProlongationDampingFactor(),
840 *levelContext.pinfo, *levelContext.redist);
843 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
844 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
845 matrices_->getProlongationDampingFactor(),
846 *levelContext.pinfo);
850 if(processNextLevel) {
851 --levelContext.update;
855 *levelContext.update += *levelContext.lhs;
858 template<
class M,
class X,
class S,
class PI,
class A>
861 return IsDirectSolver< CoarseSolver>::value;
864 template<
class M,
class X,
class S,
class PI,
class A>
866 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
868 InverseOperatorResult res;
870 if(levelContext.redist->isSetup()) {
871 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
872 if(levelContext.rhs.getRedistributed().size()>0) {
874 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
875 levelContext.rhs.getRedistributed());
876 solver_->
apply(levelContext.update.getRedistributed(),
877 levelContext.rhs.getRedistributed(), res);
879 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
880 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
882 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
883 solver_->apply(*levelContext.update, *levelContext.rhs, res);
887 coarsesolverconverged =
false;
890 presmooth(levelContext, preSteps_);
892 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
893 bool processNextLevel = moveToCoarseLevel(levelContext);
895 if(processNextLevel) {
897 for(std::size_t i=0; i<gamma_; i++)
901 moveToFineLevel(levelContext, processNextLevel);
906 if(levelContext.matrix == matrices_->matrices().finest()) {
907 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
908 if(!coarsesolverconverged){
913 postsmooth(levelContext, postSteps_);
918 template<
class M,
class X,
class S,
class PI,
class A>
919 void AMGCPR<M,X,S,PI,A>::additiveMgc(){
922 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
923 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
924 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
925 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
927 for(
typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
929 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
930 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
936 lhs = lhs_->finest();
937 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
939 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
942 smoother->apply(*lhs, *rhs);
946 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
947 InverseOperatorResult res;
948 pinfo->copyOwnerToAll(*rhs, *rhs);
949 solver_->apply(*lhs, *rhs, res);
952 DUNE_THROW(MathError,
"Coarse solver did not converge");
960 for(
typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
961 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
962 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
968 template<
class M,
class X,
class S,
class PI,
class A>
972 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
973 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
974 Iterator coarsest = smoothers_->coarsest();
975 Iterator smoother = smoothers_->finest();
976 DIterator lhs = lhs_->finest();
977 if(smoothers_->levels()>0) {
978 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
979 smoother->post(*lhs);
980 if(smoother!=coarsest)
981 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
982 smoother->post(*lhs);
983 smoother->post(*lhs);
993 template<
class M,
class X,
class S,
class PI,
class A>
997 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:294
AMGCPR(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amgcpr.hh:414
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amgcpr.hh:318
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amgcpr.hh:322
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:302
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:723
void pre(Domain &x, Range &b)
Definition: amgcpr.hh:625
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amgcpr.hh:298
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amgcpr.hh:995
void post(Domain &x)
Definition: amgcpr.hh:969
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:468
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amgcpr.hh:310
virtual void update()
Update the coarse solver and the hierarchies.
Definition: amgcpr.hh:474
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:326
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amgcpr.hh:314
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amgcpr.hh:859
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amgcpr.hh:306