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/istl/paamg/smoother.hh>
13 #include <dune/istl/paamg/transfer.hh>
14 #include <dune/istl/paamg/hierarchy.hh>
15 #include <dune/istl/solvers.hh>
16 #include <dune/istl/scalarproducts.hh>
17 #include <dune/istl/superlu.hh>
18 #include <dune/istl/umfpack.hh>
19 #include <dune/istl/solvertype.hh>
20 #include <dune/common/typetraits.hh>
21 #include <dune/common/exceptions.hh>
32 template<
class M,
class T>
33 void redistributeMatrixAmg(M&, M&, SequentialInformation&, SequentialInformation&, T&)
38 template<
class M,
class PI>
39 typename std::enable_if<!std::is_same<PI,SequentialInformation>::value,
void>::type
40 redistributeMatrixAmg(M& mat, M& matRedist, PI& info, PI& infoRedist,
41 Dune::RedistributeInformation<PI>& redistInfo)
43 info.buildGlobalLookup(mat.N());
44 redistributeMatrixEntries(mat, matRedist, info, infoRedist, redistInfo);
45 info.freeGlobalLookup();
66 template<
class M,
class X,
class S,
class P,
class K,
class A>
82 template<
class M,
class X,
class S,
class PI=SequentialInformation,
83 class A=std::allocator<X> >
86 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
132 const SmootherArgs& smootherArgs,
const Parameters& parms);
179 std::size_t levels();
181 std::size_t maxlevels();
193 auto copyFlags = NegateSet<typename PI::OwnerSet>();
194 const auto& matrices = matrices_->matrices();
195 const auto& aggregatesMapHierarchy = matrices_->aggregatesMaps();
196 const auto& infoHierarchy = matrices_->parallelInformation();
197 const auto& redistInfoHierarchy = matrices_->redistributeInformation();
198 BaseGalerkinProduct productBuilder;
199 auto aggregatesMap = aggregatesMapHierarchy.begin();
200 auto info = infoHierarchy.finest();
201 auto redistInfo = redistInfoHierarchy.begin();
202 auto matrix = matrices.finest();
203 auto coarsestMatrix = matrices.coarsest();
205 using Matrix =
typename M::matrix_type;
208 if(matrix.isRedistributed()) {
209 redistributeMatrixAmg(
const_cast<Matrix&
>(matrix->getmat()),
210 const_cast<Matrix&
>(matrix.getRedistributed().getmat()),
211 const_cast<PI&
>(*info),
const_cast<PI&
>(info.getRedistributed()),
212 const_cast<Dune::RedistributeInformation<PI>&
>(*redistInfo));
216 for(; matrix!=coarsestMatrix; ++aggregatesMap) {
217 const Matrix& fine = (matrix.isRedistributed() ? matrix.getRedistributed() : *matrix).getmat();
221 productBuilder.calculate(fine, *(*aggregatesMap),
const_cast<Matrix&
>(matrix->getmat()), *info, copyFlags);
223 if(matrix.isRedistributed()) {
224 redistributeMatrixAmg(
const_cast<Matrix&
>(matrix->getmat()),
225 const_cast<Matrix&
>(matrix.getRedistributed().getmat()),
226 const_cast<PI&
>(*info),
const_cast<PI&
>(info.getRedistributed()),
227 const_cast<Dune::RedistributeInformation<PI>&
>(*redistInfo));
257 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7 )
259 void createHierarchies(C& criterion,
Operator& matrix,
263 std::shared_ptr< Operator > op( &matrix, [](
Operator* ) {});
264 std::shared_ptr< PI > pifo(
const_cast< PI*
> (&pinfo), []( PI * ) {});
265 createHierarchies( criterion, op, pifo);
269 void createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
270 std::shared_ptr< PI > pinfo );
274 void createHierarchies(C& criterion,
Operator& matrix,
278 void setupCoarseSolver();
296 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
300 typename ParallelInformationHierarchy::Iterator
pinfo;
304 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
308 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
312 typename Hierarchy<Domain,A>::Iterator
lhs;
316 typename Hierarchy<Domain,A>::Iterator
update;
320 typename Hierarchy<Range,A>::Iterator
rhs;
332 void mgc(LevelContext& levelContext);
342 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
348 bool moveToCoarseLevel(LevelContext& levelContext);
354 void initIteratorsWithFineLevel(LevelContext& levelContext);
357 std::shared_ptr<OperatorHierarchy> matrices_;
361 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
363 std::shared_ptr<CoarseSolver> solver_;
365 std::shared_ptr< Hierarchy<Range,A> > rhs_;
367 std::shared_ptr< Hierarchy<Domain,A> > lhs_;
369 std::shared_ptr< Hierarchy<Domain,A> > update_;
371 using ScalarProduct = Dune::ScalarProduct<X>;
373 std::shared_ptr<ScalarProduct> scalarProduct_;
377 std::size_t preSteps_;
379 std::size_t postSteps_;
380 bool buildHierarchy_;
382 bool coarsesolverconverged;
383 std::shared_ptr<Smoother> coarseSmoother_;
385 SolverCategory::Category category_;
387 std::size_t verbosity_;
390 template<
class M,
class X,
class S,
class PI,
class A>
392 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
393 smoothers_(amg.smoothers_), solver_(amg.solver_),
394 rhs_(), lhs_(), update_(),
395 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
396 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
397 buildHierarchy_(amg.buildHierarchy_),
398 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
399 coarseSmoother_(amg.coarseSmoother_),
400 category_(amg.category_),
401 verbosity_(amg.verbosity_)
404 rhs_.reset(
new Hierarchy<Range,A>(*amg.rhs_) );
406 lhs_.reset(
new Hierarchy<Domain,A>(*amg.lhs_) );
408 update_.reset(
new Hierarchy<Domain,A>(*amg.update_) );
411 template<
class M,
class X,
class S,
class PI,
class A>
414 const Parameters& parms)
415 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
416 smoothers_(new Hierarchy<
Smoother,A>), solver_(&coarseSolver),
417 rhs_(), lhs_(), update_(), scalarProduct_(0),
418 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
419 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
420 additive(parms.getAdditive()), coarsesolverconverged(true),
423 category_(SolverCategory::category(*smoothers_->coarsest())),
424 verbosity_(parms.debugLevel())
426 assert(matrices_->isBuilt());
429 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
432 template<
class M,
class X,
class S,
class PI,
class A>
438 : smootherArgs_(smootherArgs),
439 smoothers_(new Hierarchy<
Smoother,A>), solver_(),
440 rhs_(), lhs_(), update_(), scalarProduct_(),
441 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
442 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
443 additive(criterion.getAdditive()), coarsesolverconverged(true),
445 category_(SolverCategory::category(pinfo)),
446 verbosity_(criterion.debugLevel())
448 if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
449 DUNE_THROW(InvalidSolverCategory,
"Matrix and Communication must have the same SolverCategory!");
450 createHierarchies(criterion,
const_cast<Operator&
>(matrix), pinfo);
453 template<
class M,
class X,
class S,
class PI,
class A>
456 if(buildHierarchy_) {
460 coarseSmoother_.reset();
464 template<
class M,
class X,
class S,
class PI,
class A>
471 template<
class M,
class X,
class S,
class PI,
class A>
475 smoothers_.reset(
new Hierarchy<Smoother,A>);
477 coarseSmoother_.reset();
478 scalarProduct_.reset();
479 buildHierarchy_=
true;
480 coarsesolverconverged =
true;
481 smoothers_.reset(
new Hierarchy<Smoother,A>);
482 recalculateHierarchy();
483 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
485 if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
486 std::cout <<
"Recalculating galerkin and coarse smoothers "<< matrices_->maxlevels() <<
" levels "
487 << watch.elapsed() <<
" seconds." << std::endl;
491 template<
class M,
class X,
class S,
class PI,
class A>
493 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
495 std::shared_ptr< PI > pinfo )
501 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
503 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
506 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
508 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
509 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
510 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
513 template<
class M,
class X,
class S,
class PI,
class A>
514 void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
520 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
521 && ( ! matrices_->redistributeInformation().back().isSetup() ||
522 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
525 SmootherArgs sargs(smootherArgs_);
526 sargs.iterations = 1;
528 typename ConstructionTraits<Smoother>::Arguments cargs;
529 cargs.setArgs(sargs);
530 if(matrices_->redistributeInformation().back().isSetup()) {
532 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
533 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
535 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
536 cargs.setComm(*matrices_->parallelInformation().coarsest());
539 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
540 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
542 coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
545 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
548 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
551 if( SolverSelector::isDirectSolver &&
552 (std::is_same<ParallelInformation,SequentialInformation>::value
553 || matrices_->parallelInformation().coarsest()->communicator().size()==1
554 || (matrices_->parallelInformation().coarsest().isRedistributed()
555 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
556 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
559 if(matrices_->parallelInformation().coarsest().isRedistributed())
561 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
564 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
571 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
573 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
574 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
578 if(matrices_->parallelInformation().coarsest().isRedistributed())
580 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
585 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
587 reinterpret_cast<typename
588 std::conditional<std::is_same<PI,SequentialInformation>::value,
589 Dune::SeqScalarProduct<X>,
590 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
591 *coarseSmoother_, 1E-2, 1000, 0));
596 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
598 reinterpret_cast<typename
599 std::conditional<std::is_same<PI,SequentialInformation>::value,
600 Dune::SeqScalarProduct<X>,
601 Dune::OverlappingSchwarzScalarProduct<X,PI>
>::type&>(*scalarProduct_),
602 *coarseSmoother_, 1E-2, 1000, 0));
622 template<
class M,
class X,
class S,
class PI,
class A>
629 typedef typename M::matrix_type Matrix;
630 typedef typename Matrix::ConstRowIterator RowIter;
631 typedef typename Matrix::ConstColIterator ColIter;
632 typedef typename Matrix::block_type Block;
634 zero=
typename Matrix::field_type();
636 const Matrix& mat=matrices_->matrices().finest()->getmat();
637 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
638 bool isDirichlet =
true;
639 bool hasDiagonal =
false;
641 for(ColIter col=row->begin(); col!=row->end(); ++col) {
642 if(row.index()==col.index()) {
650 if(isDirichlet && hasDiagonal)
651 diagonal.solve(x[row.index()], b[row.index()]);
654 if(smoothers_->levels()>0)
655 smoothers_->finest()->pre(x,b);
658 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
661 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
662 typedef std::shared_ptr< Range > RangePtr ;
663 typedef std::shared_ptr< Domain > DomainPtr;
665 typedef Range* RangePtr;
666 typedef Domain* DomainPtr;
670 RangePtr copy(
new Range(b) );
671 rhs_.reset(
new Hierarchy<Range,A>(copy) );
672 DomainPtr dcopy(
new Domain(x) );
673 lhs_.reset(
new Hierarchy<Domain,A>(dcopy) );
674 DomainPtr dcopy2(
new Domain(x) );
675 update_.reset(
new Hierarchy<Domain,A>(dcopy2) );
677 matrices_->coarsenVector(*rhs_);
678 matrices_->coarsenVector(*lhs_);
679 matrices_->coarsenVector(*update_);
682 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
683 typedef typename Hierarchy<Range,A>::Iterator RIterator;
684 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
685 Iterator coarsest = smoothers_->coarsest();
686 Iterator smoother = smoothers_->finest();
687 RIterator rhs = rhs_->finest();
688 DIterator lhs = lhs_->finest();
689 if(smoothers_->levels()>0) {
691 assert(lhs_->levels()==rhs_->levels());
692 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
693 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
695 if(smoother!=coarsest)
696 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
697 smoother->pre(*lhs,*rhs);
698 smoother->pre(*lhs,*rhs);
708 template<
class M,
class X,
class S,
class PI,
class A>
711 return matrices_->levels();
713 template<
class M,
class X,
class S,
class PI,
class A>
714 std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
716 return matrices_->maxlevels();
720 template<
class M,
class X,
class S,
class PI,
class A>
723 LevelContext levelContext;
731 initIteratorsWithFineLevel(levelContext);
734 *levelContext.lhs = v;
735 *levelContext.rhs = d;
736 *levelContext.update=0;
737 levelContext.level=0;
741 if(postSteps_==0||matrices_->maxlevels()==1)
742 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
744 v=*levelContext.update;
749 template<
class M,
class X,
class S,
class PI,
class A>
752 levelContext.smoother = smoothers_->finest();
753 levelContext.matrix = matrices_->matrices().finest();
754 levelContext.pinfo = matrices_->parallelInformation().finest();
755 levelContext.redist =
756 matrices_->redistributeInformation().begin();
757 levelContext.aggregates = matrices_->aggregatesMaps().begin();
758 levelContext.lhs = lhs_->finest();
759 levelContext.
update = update_->finest();
760 levelContext.rhs = rhs_->finest();
763 template<
class M,
class X,
class S,
class PI,
class A>
764 bool AMGCPR<M,X,S,PI,A>
765 ::moveToCoarseLevel(LevelContext& levelContext)
768 bool processNextLevel=
true;
770 if(levelContext.redist->isSetup()) {
771 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
772 levelContext.rhs.getRedistributed());
773 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
774 if(processNextLevel) {
776 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
777 ++levelContext.pinfo;
778 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
779 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
780 static_cast<const Range&
>(fineRhs.getRedistributed()),
781 *levelContext.pinfo);
785 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
786 ++levelContext.pinfo;
787 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
788 ::restrictVector(*(*levelContext.aggregates),
789 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
790 *levelContext.pinfo);
793 if(processNextLevel) {
796 ++levelContext.update;
797 ++levelContext.matrix;
798 ++levelContext.level;
799 ++levelContext.redist;
801 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
803 ++levelContext.smoother;
804 ++levelContext.aggregates;
807 *levelContext.update=0;
809 return processNextLevel;
812 template<
class M,
class X,
class S,
class PI,
class A>
813 void AMGCPR<M,X,S,PI,A>
814 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
816 if(processNextLevel) {
817 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
819 --levelContext.smoother;
820 --levelContext.aggregates;
822 --levelContext.redist;
823 --levelContext.level;
825 --levelContext.matrix;
829 --levelContext.pinfo;
831 if(levelContext.redist->isSetup()) {
833 levelContext.lhs.getRedistributed()=0;
834 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
835 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
836 levelContext.lhs.getRedistributed(),
837 matrices_->getProlongationDampingFactor(),
838 *levelContext.pinfo, *levelContext.redist);
841 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
842 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
843 matrices_->getProlongationDampingFactor(),
844 *levelContext.pinfo);
848 if(processNextLevel) {
849 --levelContext.update;
853 *levelContext.update += *levelContext.lhs;
856 template<
class M,
class X,
class S,
class PI,
class A>
859 return IsDirectSolver< CoarseSolver>::value;
862 template<
class M,
class X,
class S,
class PI,
class A>
864 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
866 InverseOperatorResult res;
868 if(levelContext.redist->isSetup()) {
869 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
870 if(levelContext.rhs.getRedistributed().size()>0) {
872 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
873 levelContext.rhs.getRedistributed());
874 solver_->
apply(levelContext.update.getRedistributed(),
875 levelContext.rhs.getRedistributed(), res);
877 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
878 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
880 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
881 solver_->apply(*levelContext.update, *levelContext.rhs, res);
885 coarsesolverconverged =
false;
888 presmooth(levelContext, preSteps_);
890 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
891 bool processNextLevel = moveToCoarseLevel(levelContext);
893 if(processNextLevel) {
895 for(std::size_t i=0; i<gamma_; i++)
899 moveToFineLevel(levelContext, processNextLevel);
904 if(levelContext.matrix == matrices_->matrices().finest()) {
905 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
906 if(!coarsesolverconverged){
911 postsmooth(levelContext, postSteps_);
916 template<
class M,
class X,
class S,
class PI,
class A>
917 void AMGCPR<M,X,S,PI,A>::additiveMgc(){
920 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
921 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
922 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
923 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
925 for(
typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
927 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
928 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
934 lhs = lhs_->finest();
935 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
937 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
940 smoother->apply(*lhs, *rhs);
944 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
945 InverseOperatorResult res;
946 pinfo->copyOwnerToAll(*rhs, *rhs);
947 solver_->apply(*lhs, *rhs, res);
950 DUNE_THROW(MathError,
"Coarse solver did not converge");
958 for(
typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
959 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
960 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
966 template<
class M,
class X,
class S,
class PI,
class A>
969 DUNE_UNUSED_PARAMETER(x);
971 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
972 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
973 Iterator coarsest = smoothers_->coarsest();
974 Iterator smoother = smoothers_->finest();
975 DIterator lhs = lhs_->finest();
976 if(smoothers_->levels()>0) {
977 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
978 smoother->post(*lhs);
979 if(smoother!=coarsest)
980 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
981 smoother->post(*lhs);
982 smoother->post(*lhs);
992 template<
class M,
class X,
class S,
class PI,
class A>
996 matrices_->getCoarsestAggregatesOnFinest(cont);
Parallel algebraic multigrid based on agglomeration.
Definition: amgcpr.hh:85
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
M Operator
The matrix operator type.
Definition: amgcpr.hh:93
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amgcpr.hh:292
AMGCPR(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amgcpr.hh:412
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amgcpr.hh:316
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amgcpr.hh:320
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amgcpr.hh:102
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amgcpr.hh:300
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amgcpr.hh:164
void apply(Domain &v, const Range &d)
Definition: amgcpr.hh:721
void pre(Domain &x, Range &b)
Definition: amgcpr.hh:623
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amgcpr.hh:296
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amgcpr.hh:994
void post(Domain &x)
Definition: amgcpr.hh:967
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amgcpr.hh:104
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amgcpr.hh:191
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amgcpr.hh:120
X Range
The range type.
Definition: amgcpr.hh:109
X Domain
The domain type.
Definition: amgcpr.hh:107
S Smoother
The type of the smoother.
Definition: amgcpr.hh:117
void updateSolver(C &criterion, const Operator &, const PI &pinfo)
Update the coarse solver and the hierarchies.
Definition: amgcpr.hh:466
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amgcpr.hh:308
virtual void update()
Update the coarse solver and the hierarchies.
Definition: amgcpr.hh:472
PI ParallelInformation
The type of the parallel information.
Definition: amgcpr.hh:100
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amgcpr.hh:111
std::size_t level
The level index.
Definition: amgcpr.hh:324
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amgcpr.hh:312
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amgcpr.hh:857
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amgcpr.hh:304