5 #ifndef DUNE_AMG_AMG_HH
6 #define DUNE_AMG_AMG_HH
10 #include <dune/common/exceptions.hh>
19 #include <dune/common/typetraits.hh>
20 #include <dune/common/exceptions.hh>
21 #include <dune/common/scalarvectorview.hh>
22 #include <dune/common/scalarmatrixview.hh>
23 #include <dune/common/parametertree.hh>
46 template<
class M,
class X,
class S,
class P,
class K,
class A>
62 template<
class M,
class X,
class S,
class PI=SequentialInformation,
63 class A=std::allocator<X> >
66 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
223 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
241 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
242 const PI& pinfo,
const Norm&,
243 const ParameterTree& configuration,
244 std::true_type compiles = std::true_type());
246 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
247 const PI& pinfo,
const Norm&,
248 const ParameterTree& configuration,
255 void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
256 const PI& pinfo,
const ParameterTree& configuration);
264 void createHierarchies(C& criterion,
265 const std::shared_ptr<const Operator>& matrixptr,
291 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
295 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
319 void mgc(LevelContext& levelContext);
329 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
335 bool moveToCoarseLevel(LevelContext& levelContext);
341 void initIteratorsWithFineLevel(LevelContext& levelContext);
344 std::shared_ptr<OperatorHierarchy> matrices_;
348 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
350 std::shared_ptr<CoarseSolver> solver_;
352 std::shared_ptr<Hierarchy<Range,A>> rhs_;
354 std::shared_ptr<Hierarchy<Domain,A>> lhs_;
356 std::shared_ptr<Hierarchy<Domain,A>> update_;
360 std::shared_ptr<ScalarProduct> scalarProduct_;
364 std::size_t preSteps_;
366 std::size_t postSteps_;
367 bool buildHierarchy_;
369 bool coarsesolverconverged;
370 std::shared_ptr<Smoother> coarseSmoother_;
374 std::size_t verbosity_;
380 std::stringstream retval;
381 std::ostream_iterator<char> out(retval);
382 std::transform(str.begin(), str.end(), out,
384 return std::tolower(c, std::locale::classic());
391 template<
class M,
class X,
class S,
class PI,
class A>
393 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
394 smoothers_(amg.smoothers_), solver_(amg.solver_),
395 rhs_(), lhs_(), update_(),
396 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
397 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
398 buildHierarchy_(amg.buildHierarchy_),
399 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
400 coarseSmoother_(amg.coarseSmoother_),
401 category_(amg.category_),
402 verbosity_(amg.verbosity_)
405 template<
class M,
class X,
class S,
class PI,
class A>
409 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
411 rhs_(), lhs_(), update_(), scalarProduct_(0),
412 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
413 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
414 additive(parms.getAdditive()), coarsesolverconverged(true),
418 verbosity_(parms.debugLevel())
420 assert(matrices_->isBuilt());
423 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
426 template<
class M,
class X,
class S,
class PI,
class A>
432 : smootherArgs_(smootherArgs),
434 rhs_(), lhs_(), update_(), scalarProduct_(),
435 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
436 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
437 additive(criterion.getAdditive()), coarsesolverconverged(true),
440 verbosity_(criterion.debugLevel())
447 auto matrixptr = stackobject_to_shared_ptr(matrix);
448 createHierarchies(criterion, matrixptr, pinfo);
451 template<
class M,
class X,
class S,
class PI,
class A>
453 const ParameterTree& configuration,
456 solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
457 coarsesolverconverged(true), coarseSmoother_(),
461 if (configuration.hasKey (
"smootherIterations"))
462 smootherArgs_.iterations = configuration.get<
int>(
"smootherIterations");
464 if (configuration.hasKey (
"smootherRelaxation"))
465 smootherArgs_.relaxationFactor = configuration.get<
typename SmootherArgs::RelaxationFactor>(
"smootherRelaxation");
467 auto normName = ToLower()(configuration.get(
"strengthMeasure",
"diagonal"));
468 auto index = configuration.get<
int>(
"diagonalRowIndex", 0);
470 if ( normName ==
"diagonal")
473 using real_type =
typename FieldTraits<field_type>::real_type;
474 std::is_convertible<field_type, real_type> compiles;
479 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<0>(), configuration, compiles);
482 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<1>(), configuration, compiles);
485 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<2>(), configuration, compiles);
488 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<3>(), configuration, compiles);
491 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<4>(), configuration, compiles);
494 DUNE_THROW(InvalidStateException,
"Currently strengthIndex>4 is not supported.");
497 else if (normName ==
"rowsum")
498 createCriterionAndHierarchies(matrixptr, pinfo,
RowSum(), configuration);
499 else if (normName ==
"frobenius")
500 createCriterionAndHierarchies(matrixptr, pinfo,
FrobeniusNorm(), configuration);
501 else if (normName ==
"one")
502 createCriterionAndHierarchies(matrixptr, pinfo,
AlwaysOneNorm(), configuration);
504 DUNE_THROW(Dune::NotImplemented,
"Wrong config file: strengthMeasure "<<normName<<
" is not supported by AMG");
507 template<
class M,
class X,
class S,
class PI,
class A>
511 DUNE_THROW(InvalidStateException,
"Strength of connection measure does not support this type ("
512 << className<typename M::field_type>() <<
") as it is lacking a conversion to"
513 << className<
typename FieldTraits<typename M::field_type>::real_type>() <<
".");
516 template<
class M,
class X,
class S,
class PI,
class A>
518 void AMG<M,X,S,PI,A>::createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
const PI& pinfo,
const Norm&,
const ParameterTree& configuration, std::true_type)
520 if (configuration.get<
bool>(
"criterionSymmetric",
true))
525 createHierarchies(criterion, matrixptr, pinfo, configuration);
532 createHierarchies(criterion, matrixptr, pinfo, configuration);
536 template<
class M,
class X,
class S,
class PI,
class A>
538 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
const PI& pinfo,
const ParameterTree& configuration)
540 if (configuration.hasKey (
"maxLevel"))
541 criterion.setMaxLevel(configuration.get<
int>(
"maxLevel"));
543 if (configuration.hasKey (
"minCoarseningRate"))
544 criterion.setMinCoarsenRate(configuration.get<
int>(
"minCoarseningRate"));
546 if (configuration.hasKey (
"coarsenTarget"))
547 criterion.setCoarsenTarget (configuration.get<
int>(
"coarsenTarget"));
549 if (configuration.hasKey (
"accumulationMode"))
551 std::string mode = ToLower()(configuration.get<std::string>(
"accumulationMode"));
554 else if ( mode ==
"atonce" )
556 else if ( mode ==
"successive")
559 DUNE_THROW(InvalidSolverFactoryConfiguration,
"Parameter accumulationMode does not allow value "
563 if (configuration.hasKey (
"prolongationDampingFactor"))
564 criterion.setProlongationDampingFactor (configuration.get<
double>(
"prolongationDampingFactor"));
566 if (configuration.hasKey(
"defaultAggregationSizeMode"))
568 auto mode = ToLower()(configuration.get<std::string>(
"defaultAggregationSizeMode"));
569 auto dim = configuration.get<std::size_t>(
"defaultAggregationDimension");
570 std::size_t maxDistance = 2;
571 if (configuration.hasKey(
"MaxAggregateDistance"))
572 maxDistance = configuration.get<std::size_t>(
"maxAggregateDistance");
573 if (mode ==
"isotropic")
574 criterion.setDefaultValuesIsotropic(dim, maxDistance);
575 else if(mode ==
"anisotropic")
576 criterion.setDefaultValuesAnisotropic(dim, maxDistance);
578 DUNE_THROW(InvalidSolverFactoryConfiguration,
"Parameter accumulationMode does not allow value "
582 if (configuration.hasKey(
"maxAggregateDistance"))
583 criterion.setMaxDistance(configuration.get<std::size_t>(
"maxAggregateDistance"));
585 if (configuration.hasKey(
"minAggregateSize"))
586 criterion.setMinAggregateSize(configuration.get<std::size_t>(
"minAggregateSize"));
588 if (configuration.hasKey(
"maxAggregateSize"))
589 criterion.setMaxAggregateSize(configuration.get<std::size_t>(
"maxAggregateSize"));
591 if (configuration.hasKey(
"maxAggregateConnectivity"))
592 criterion.setMaxConnectivity(configuration.get<std::size_t>(
"maxAggregateConnectivity"));
594 if (configuration.hasKey (
"alpha"))
595 criterion.setAlpha (configuration.get<
double> (
"alpha"));
597 if (configuration.hasKey (
"beta"))
598 criterion.setBeta (configuration.get<
double> (
"beta"));
600 if (configuration.hasKey (
"gamma"))
601 criterion.setGamma (configuration.get<std::size_t> (
"gamma"));
602 gamma_ = criterion.getGamma();
604 if (configuration.hasKey (
"additive"))
605 criterion.setAdditive (configuration.get<
bool>(
"additive"));
606 additive = criterion.getAdditive();
608 if (configuration.hasKey (
"preSteps"))
609 criterion.setNoPreSmoothSteps (configuration.get<std::size_t> (
"preSteps"));
610 preSteps_ = criterion.getNoPreSmoothSteps ();
612 if (configuration.hasKey (
"postSteps"))
613 criterion.setNoPostSmoothSteps (configuration.get<std::size_t> (
"postSteps"));
614 postSteps_ = criterion.getNoPostSmoothSteps ();
616 verbosity_ = configuration.get(
"verbosity", 0);
617 criterion.setDebugLevel (verbosity_);
619 createHierarchies(criterion, matrixptr, pinfo);
622 template <
class Matrix,
630 #if DISABLE_AMG_DIRECTSOLVER
632 #elif HAVE_SUITESPARSE_UMFPACK
640 template <
class M, SolverType>
646 DUNE_THROW(NotImplemented,
"DirectSolver not selected");
649 static std::string
name () {
return "None"; }
651 #if HAVE_SUITESPARSE_UMFPACK
656 static type*
create(
const M&
mat,
bool verbose,
bool reusevector )
658 return new type(
mat, verbose, reusevector );
660 static std::string
name () {
return "UMFPack"; }
670 return new type(
mat, verbose, reusevector );
672 static std::string
name () {
return "SuperLU"; }
687 template<
class M,
class X,
class S,
class PI,
class A>
689 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
690 const std::shared_ptr<const Operator>& matrixptr,
694 matrices_ = std::make_shared<OperatorHierarchy>(
695 std::const_pointer_cast<Operator>(matrixptr),
696 stackobject_to_shared_ptr(
const_cast<PI&
>(pinfo)));
698 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
701 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
707 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
708 && ( ! matrices_->redistributeInformation().back().isSetup() ||
709 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
712 SmootherArgs sargs(smootherArgs_);
713 sargs.iterations = 1;
716 cargs.setArgs(sargs);
717 if(matrices_->redistributeInformation().back().isSetup()) {
719 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
720 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
722 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
723 cargs.setComm(*matrices_->parallelInformation().coarsest());
727 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
729 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
732 if( SolverSelector::isDirectSolver &&
733 (std::is_same<ParallelInformation,SequentialInformation>::value
734 || matrices_->parallelInformation().coarsest()->communicator().size()==1
735 || (matrices_->parallelInformation().coarsest().isRedistributed()
736 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
737 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
740 if(matrices_->parallelInformation().coarsest().isRedistributed())
742 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
745 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
752 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
754 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
755 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
759 if(matrices_->parallelInformation().coarsest().isRedistributed())
761 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
766 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
768 *coarseSmoother_, 1E-2, 1000, 0));
773 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
775 *coarseSmoother_, 1E-2, 1000, 0));
794 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
795 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
796 <<
"(including coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
800 template<
class M,
class X,
class S,
class PI,
class A>
807 typedef typename M::matrix_type
Matrix;
814 const Matrix&
mat=matrices_->matrices().finest()->getmat();
815 for(RowIter row=
mat.begin(); row!=
mat.end(); ++row) {
816 bool isDirichlet =
true;
817 bool hasDiagonal =
false;
819 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
820 if(row.index()==
col.index()) {
828 if(isDirichlet && hasDiagonal)
830 auto&& xEntry = Impl::asVector(x[row.index()]);
831 auto&& bEntry = Impl::asVector(b[row.index()]);
832 Impl::asMatrix(diagonal).solve(xEntry, bEntry);
836 if(smoothers_->levels()>0)
837 smoothers_->finest()->pre(x,b);
840 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
841 rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
842 lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
843 update_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
844 matrices_->coarsenVector(*rhs_);
845 matrices_->coarsenVector(*lhs_);
846 matrices_->coarsenVector(*update_);
852 Iterator coarsest = smoothers_->
coarsest();
853 Iterator smoother = smoothers_->finest();
854 RIterator rhs = rhs_->finest();
855 DIterator lhs = lhs_->finest();
856 if(smoothers_->levels()>1) {
858 assert(lhs_->levels()==rhs_->levels());
859 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
860 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
862 if(smoother!=coarsest)
863 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
864 smoother->pre(*lhs,*rhs);
865 smoother->pre(*lhs,*rhs);
875 template<
class M,
class X,
class S,
class PI,
class A>
878 return matrices_->levels();
880 template<
class M,
class X,
class S,
class PI,
class A>
883 return matrices_->maxlevels();
887 template<
class M,
class X,
class S,
class PI,
class A>
890 LevelContext levelContext;
898 initIteratorsWithFineLevel(levelContext);
901 *levelContext.lhs = v;
902 *levelContext.rhs = d;
903 *levelContext.update=0;
904 levelContext.level=0;
908 if(postSteps_==0||matrices_->maxlevels()==1)
909 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
911 v=*levelContext.update;
916 template<
class M,
class X,
class S,
class PI,
class A>
919 levelContext.smoother = smoothers_->finest();
920 levelContext.matrix = matrices_->matrices().finest();
921 levelContext.pinfo = matrices_->parallelInformation().finest();
922 levelContext.redist =
923 matrices_->redistributeInformation().begin();
924 levelContext.aggregates = matrices_->aggregatesMaps().begin();
925 levelContext.lhs = lhs_->finest();
926 levelContext.update = update_->finest();
927 levelContext.rhs = rhs_->finest();
930 template<
class M,
class X,
class S,
class PI,
class A>
932 ::moveToCoarseLevel(LevelContext& levelContext)
935 bool processNextLevel=
true;
937 if(levelContext.redist->isSetup()) {
938 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
939 levelContext.rhs.getRedistributed());
940 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
941 if(processNextLevel) {
944 ++levelContext.pinfo;
947 static_cast<const Range&
>(fineRhs.getRedistributed()),
948 *levelContext.pinfo);
953 ++levelContext.pinfo;
956 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
957 *levelContext.pinfo);
960 if(processNextLevel) {
963 ++levelContext.update;
964 ++levelContext.matrix;
965 ++levelContext.level;
966 ++levelContext.redist;
968 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
970 ++levelContext.smoother;
971 ++levelContext.aggregates;
974 *levelContext.update=0;
976 return processNextLevel;
979 template<
class M,
class X,
class S,
class PI,
class A>
981 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
983 if(processNextLevel) {
984 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
986 --levelContext.smoother;
987 --levelContext.aggregates;
989 --levelContext.redist;
990 --levelContext.level;
992 --levelContext.matrix;
996 --levelContext.pinfo;
998 if(levelContext.redist->isSetup()) {
1000 levelContext.lhs.getRedistributed()=0;
1002 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1003 levelContext.lhs.getRedistributed(),
1004 matrices_->getProlongationDampingFactor(),
1005 *levelContext.pinfo, *levelContext.redist);
1007 *levelContext.lhs=0;
1009 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1010 matrices_->getProlongationDampingFactor(),
1011 *levelContext.pinfo);
1015 if(processNextLevel) {
1016 --levelContext.update;
1020 *levelContext.update += *levelContext.lhs;
1023 template<
class M,
class X,
class S,
class PI,
class A>
1029 template<
class M,
class X,
class S,
class PI,
class A>
1031 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
1035 if(levelContext.redist->isSetup()) {
1036 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
1037 if(levelContext.rhs.getRedistributed().size()>0) {
1039 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
1040 levelContext.rhs.getRedistributed());
1041 solver_->apply(levelContext.update.getRedistributed(),
1042 levelContext.rhs.getRedistributed(), res);
1044 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
1045 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
1047 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
1048 solver_->apply(*levelContext.update, *levelContext.rhs, res);
1052 coarsesolverconverged =
false;
1057 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1058 bool processNextLevel = moveToCoarseLevel(levelContext);
1060 if(processNextLevel) {
1062 for(std::size_t i=0; i<gamma_; i++){
1064 if (levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels())
1067 levelContext.matrix->applyscaleadd(-1., *levelContext.lhs, *levelContext.rhs);
1072 moveToFineLevel(levelContext, processNextLevel);
1077 if(levelContext.matrix == matrices_->matrices().finest()) {
1078 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
1079 if(!coarsesolverconverged)
1080 DUNE_THROW(MathError,
"Coarse solver did not converge");
1088 template<
class M,
class X,
class S,
class PI,
class A>
1089 void AMG<M,X,S,PI,A>::additiveMgc(){
1092 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
1095 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
1100 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
1106 lhs = lhs_->finest();
1109 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
1112 smoother->apply(*lhs, *rhs);
1116 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1117 InverseOperatorResult res;
1118 pinfo->copyOwnerToAll(*rhs, *rhs);
1119 solver_->apply(*lhs, *rhs, res);
1122 DUNE_THROW(MathError,
"Coarse solver did not converge");
1138 template<
class M,
class X,
class S,
class PI,
class A>
1144 Iterator coarsest = smoothers_->
coarsest();
1145 Iterator smoother = smoothers_->finest();
1146 DIterator lhs = lhs_->finest();
1147 if(smoothers_->levels()>0) {
1148 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
1149 smoother->post(*lhs);
1150 if(smoother!=coarsest)
1151 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
1152 smoother->post(*lhs);
1153 smoother->post(*lhs);
1160 template<
class M,
class X,
class S,
class PI,
class A>
1164 matrices_->getCoarsestAggregatesOnFinest(cont);
1174 std::shared_ptr<Dune::Preconditioner<typename OP::element_type::domain_type, typename OP::element_type::range_type> >
1175 makeAMG(
const OP& op,
const std::string& smoother,
const Dune::ParameterTree& config)
const
1177 DUNE_THROW(Dune::Exception,
"Operator type not supported by AMG");
1180 template<
class M,
class X,
class Y>
1181 std::shared_ptr<Dune::Preconditioner<X,Y> >
1183 const Dune::ParameterTree& config)
const
1187 if(smoother ==
"ssor")
1188 return std::make_shared<Amg::AMG<OP, X, SeqSSOR<M,X,Y>>>(op, config);
1189 if(smoother ==
"sor")
1190 return std::make_shared<Amg::AMG<OP, X, SeqSOR<M,X,Y>>>(op, config);
1191 if(smoother ==
"jac")
1192 return std::make_shared<Amg::AMG<OP, X, SeqJac<M,X,Y>>>(op, config);
1193 if(smoother ==
"gs")
1194 return std::make_shared<Amg::AMG<OP, X, SeqGS<M,X,Y>>>(op, config);
1195 if(smoother ==
"ilu")
1196 return std::make_shared<Amg::AMG<OP, X, SeqILU<M,X,Y>>>(op, config);
1198 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1201 template<
class M,
class X,
class Y,
class C>
1202 std::shared_ptr<Dune::Preconditioner<X,Y> >
1204 const Dune::ParameterTree& config)
const
1208 auto cop = std::static_pointer_cast<const OP>(op);
1210 if(smoother ==
"ssor")
1211 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1212 if(smoother ==
"sor")
1213 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1214 if(smoother ==
"jac")
1215 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqJac<M,X,Y>>,C>>(cop, config, op->getCommunication());
1216 if(smoother ==
"gs")
1217 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqGS<M,X,Y>>,C>>(cop, config, op->getCommunication());
1218 if(smoother ==
"ilu")
1219 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqILU<M,X,Y>>,C>>(cop, config, op->getCommunication());
1221 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1224 template<
class M,
class X,
class Y,
class C>
1225 std::shared_ptr<Dune::Preconditioner<X,Y> >
1227 const Dune::ParameterTree& config)
const
1231 if(smoother ==
"ssor")
1232 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1233 if(smoother ==
"sor")
1234 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1235 if(smoother ==
"jac")
1236 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqJac<M,X,Y>>,C>>(op, config, op->getCommunication());
1237 if(smoother ==
"gs")
1238 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqGS<M,X,Y>>,C>>(op, config, op->getCommunication());
1239 if(smoother ==
"ilu")
1240 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqILU<M,X,Y>>,C>>(op, config, op->getCommunication());
1242 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1245 template<
typename TL,
typename OP>
1246 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1247 typename Dune::TypeListElement<2, TL>::type>>
1248 operator() (TL tl,
const std::shared_ptr<OP>& op,
const Dune::ParameterTree& config,
1251 using field_type =
typename OP::matrix_type::field_type;
1252 using real_type =
typename FieldTraits<field_type>::real_type;
1253 if (!std::is_convertible<field_type, real_type>())
1255 className<field_type>() <<
1256 ") to be convertible to its real_type (" <<
1257 className<real_type>() <<
1259 using D =
typename Dune::TypeListElement<1, decltype(tl)>::type;
1260 using R =
typename Dune::TypeListElement<2, decltype(tl)>::type;
1261 std::shared_ptr<Preconditioner<D,R>> amg;
1262 std::string smoother = config.get(
"smoother",
"ssor");
1263 return makeAMG(op, smoother, config);
1266 template<
typename TL,
typename OP>
1267 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1268 typename Dune::TypeListElement<2, TL>::type>>
1269 operator() (TL ,
const std::shared_ptr<OP>& ,
const Dune::ParameterTree& ,
1272 DUNE_THROW(
UnsupportedType,
"AMG needs a FieldMatrix as Matrix block_type");
Provides a classes representing the hierarchies in AMG.
Classes for the generic construction and application of the smoothers.
Prolongation and restriction for amg.
Define base class for scalar product and norm.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Classes for using SuperLU with ISTL matrices.
Classes for using UMFPack with ISTL matrices.
Col col
Definition: matrixmatrix.hh:351
Matrix & mat
Definition: matrixmatrix.hh:347
AMG(const AMG &amg)
Copy constructor.
Definition: amg.hh:392
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:801
static std::string name()
Definition: amg.hh:680
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:303
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:307
std::shared_ptr< Dune::Preconditioner< typename OP::element_type::domain_type, typename OP::element_type::range_type > > makeAMG(const OP &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1175
static std::string name()
Definition: amg.hh:672
Matrix ::field_type field_type
Definition: amg.hh:626
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:1024
X Domain
The domain type.
Definition: amg.hh:87
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:681
AMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:406
AMG(std::shared_ptr< const Operator > fineOperator, const ParameterTree &configuration, const ParallelInformation &pinfo=ParallelInformation())
Constructor an AMG via ParameterTree.
Definition: amg.hh:452
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:287
SolverType
Definition: amg.hh:627
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:295
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:100
Solver< Matrix, solver > SelectedSolver
Definition: amg.hh:677
std::string operator()(const std::string &str)
Definition: amg.hh:378
S Smoother
The type of the smoother.
Definition: amg.hh:97
static std::string name()
Definition: amg.hh:649
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:279
M Operator
The matrix operator type.
Definition: amg.hh:73
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:283
SelectedSolver ::type DirectSolver
Definition: amg.hh:678
SuperLU< M > type
Definition: amg.hh:667
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:291
X Range
The range type.
Definition: amg.hh:89
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:406
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:1162
std::size_t levels()
Definition: amg.hh:876
InverseOperator< Vector, Vector > type
Definition: amg.hh:643
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< MatrixAdapter< M, X, Y >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1182
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:299
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< NonoverlappingSchwarzOperator< M, X, Y, C >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1226
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
static constexpr SolverType solver
Definition: amg.hh:629
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< OverlappingSchwarzOperator< M, X, Y, C >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1203
static constexpr bool isDirectSolver
Definition: amg.hh:679
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:221
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:383
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:84
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:91
void post(Domain &x)
Clean up.
Definition: amg.hh:1139
std::size_t maxlevels()
Definition: amg.hh:881
std::shared_ptr< Dune::Preconditioner< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL tl, const std::shared_ptr< OP > &op, const Dune::ParameterTree &config, std::enable_if_t< isValidBlockType< typename OP::matrix_type::block_type >::value, int >=0) const
Definition: amg.hh:1248
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:428
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:644
std::size_t level
The level index.
Definition: amg.hh:311
AMG(const Operator &fineOperator, const C &criterion, const SmootherArgs &smootherArgs=SmootherArgs(), const ParallelInformation &pinfo=ParallelInformation())
Construct an AMG with an inexact coarse solver based on the smoother.
Definition: amg.hh:428
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:888
Smoother SmootherType
Definition: amg.hh:275
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:82
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amg.hh:194
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:668
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:80
@ none
Definition: amg.hh:627
@ umfpack
Definition: amg.hh:627
@ superlu
Definition: amg.hh:627
@ atOnceAccu
Accumulate data to one process at once.
Definition: parameters.hh:244
@ noAccu
No data accumulution.
Definition: parameters.hh:238
@ successiveAccu
Successively accumulate to fewer processes.
Definition: parameters.hh:248
Definition: allocator.hh:11
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
ConstIterator class for sequential access.
Definition: matrix.hh:404
A generic dynamic dense matrix.
Definition: matrix.hh:561
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
T block_type
Export the type representing the components.
Definition: matrix.hh:568
Definition: matrixutils.hh:27
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:61
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:137
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:463
Definition: aggregates.hh:480
Definition: aggregates.hh:496
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:539
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:140
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:33
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:65
An overlapping Schwarz operator.
Definition: schwarz.hh:75
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:216
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:219
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:61
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:283
All parameters for AMG.
Definition: parameters.hh:393
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:66
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:32
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioner.hh:39
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Statistics about the application of an inverse operator.
Definition: solver.hh:48
bool converged
True if convergence criterion has been met.
Definition: solver.hh:73
Categories for the solvers.
Definition: solvercategory.hh:22
Category
Definition: solvercategory.hh:23
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:34
Definition: solvercategory.hh:54
Definition: solverregistry.hh:77
Definition: solvertype.hh:16
SuperLu Solver.
Definition: superlu.hh:271
Definition: umfpack.hh:49
The UMFPack direct sparse solver.
Definition: umfpack.hh:215