3 #ifndef DUNE_AMG_AMG_HH
4 #define DUNE_AMG_AMG_HH
8 #include <dune/common/exceptions.hh>
17 #include <dune/common/typetraits.hh>
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/scalarvectorview.hh>
20 #include <dune/common/scalarmatrixview.hh>
21 #include <dune/common/parametertree.hh>
44 template<
class M,
class X,
class S,
class P,
class K,
class A>
60 template<
class M,
class X,
class S,
class PI=SequentialInformation,
61 class A=std::allocator<X> >
64 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
221 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
239 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
240 const PI& pinfo,
const Norm&,
241 const ParameterTree& configuration,
242 std::true_type compiles = std::true_type());
244 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
245 const PI& pinfo,
const Norm&,
246 const ParameterTree& configuration,
253 void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
254 const PI& pinfo,
const ParameterTree& configuration);
262 void createHierarchies(C& criterion,
263 const std::shared_ptr<const Operator>& matrixptr,
289 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
293 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
317 void mgc(LevelContext& levelContext);
327 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel);
333 bool moveToCoarseLevel(LevelContext& levelContext);
339 void initIteratorsWithFineLevel(LevelContext& levelContext);
342 std::shared_ptr<OperatorHierarchy> matrices_;
346 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
348 std::shared_ptr<CoarseSolver> solver_;
350 std::shared_ptr<Hierarchy<Range,A>> rhs_;
352 std::shared_ptr<Hierarchy<Domain,A>> lhs_;
354 std::shared_ptr<Hierarchy<Domain,A>> update_;
358 std::shared_ptr<ScalarProduct> scalarProduct_;
362 std::size_t preSteps_;
364 std::size_t postSteps_;
365 bool buildHierarchy_;
367 bool coarsesolverconverged;
368 std::shared_ptr<Smoother> coarseSmoother_;
372 std::size_t verbosity_;
378 std::stringstream retval;
379 std::ostream_iterator<char> out(retval);
380 std::transform(str.begin(), str.end(), out,
382 return std::tolower(c, std::locale::classic());
389 template<
class M,
class X,
class S,
class PI,
class A>
391 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
392 smoothers_(amg.smoothers_), solver_(amg.solver_),
393 rhs_(), lhs_(), update_(),
394 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
395 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
396 buildHierarchy_(amg.buildHierarchy_),
397 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
398 coarseSmoother_(amg.coarseSmoother_),
399 category_(amg.category_),
400 verbosity_(amg.verbosity_)
403 template<
class M,
class X,
class S,
class PI,
class A>
407 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
409 rhs_(), lhs_(), update_(), scalarProduct_(0),
410 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
411 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
412 additive(parms.getAdditive()), coarsesolverconverged(true),
416 verbosity_(parms.debugLevel())
418 assert(matrices_->isBuilt());
421 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
424 template<
class M,
class X,
class S,
class PI,
class A>
430 : smootherArgs_(smootherArgs),
432 rhs_(), lhs_(), update_(), scalarProduct_(),
433 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
434 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
435 additive(criterion.getAdditive()), coarsesolverconverged(true),
438 verbosity_(criterion.debugLevel())
445 auto matrixptr = stackobject_to_shared_ptr(matrix);
446 createHierarchies(criterion, matrixptr, pinfo);
449 template<
class M,
class X,
class S,
class PI,
class A>
451 const ParameterTree& configuration,
454 solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
455 coarsesolverconverged(true), coarseSmoother_(),
459 if (configuration.hasKey (
"smootherIterations"))
460 smootherArgs_.iterations = configuration.get<
int>(
"smootherIterations");
462 if (configuration.hasKey (
"smootherRelaxation"))
463 smootherArgs_.relaxationFactor = configuration.get<
typename SmootherArgs::RelaxationFactor>(
"smootherRelaxation");
465 auto normName = ToLower()(configuration.get(
"strengthMeasure",
"diagonal"));
466 auto index = configuration.get<
int>(
"diagonalRowIndex", 0);
468 if ( normName ==
"diagonal")
471 using real_type =
typename FieldTraits<field_type>::real_type;
472 std::is_convertible<field_type, real_type> compiles;
477 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<0>(), configuration, compiles);
480 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<1>(), configuration, compiles);
483 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<2>(), configuration, compiles);
486 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<3>(), configuration, compiles);
489 createCriterionAndHierarchies(matrixptr, pinfo,
Diagonal<4>(), configuration, compiles);
492 DUNE_THROW(InvalidStateException,
"Currently strengthIndex>4 is not supported.");
495 else if (normName ==
"rowsum")
496 createCriterionAndHierarchies(matrixptr, pinfo,
RowSum(), configuration);
497 else if (normName ==
"frobenius")
498 createCriterionAndHierarchies(matrixptr, pinfo,
FrobeniusNorm(), configuration);
499 else if (normName ==
"one")
500 createCriterionAndHierarchies(matrixptr, pinfo,
AlwaysOneNorm(), configuration);
502 DUNE_THROW(Dune::NotImplemented,
"Wrong config file: strengthMeasure "<<normName<<
" is not supported by AMG");
505 template<
class M,
class X,
class S,
class PI,
class A>
509 DUNE_THROW(InvalidStateException,
"Strength of connection measure does not support this type ("
510 << className<typename M::field_type>() <<
") as it is lacking a conversion to"
511 << className<
typename FieldTraits<typename M::field_type>::real_type>() <<
".");
514 template<
class M,
class X,
class S,
class PI,
class A>
516 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)
518 if (configuration.get<
bool>(
"criterionSymmetric",
true))
523 createHierarchies(criterion, matrixptr, pinfo, configuration);
530 createHierarchies(criterion, matrixptr, pinfo, configuration);
534 template<
class M,
class X,
class S,
class PI,
class A>
536 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
const PI& pinfo,
const ParameterTree& configuration)
538 if (configuration.hasKey (
"maxLevel"))
539 criterion.setMaxLevel(configuration.get<
int>(
"maxLevel"));
541 if (configuration.hasKey (
"minCoarseningRate"))
542 criterion.setMinCoarsenRate(configuration.get<
int>(
"minCoarseningRate"));
544 if (configuration.hasKey (
"coarsenTarget"))
545 criterion.setCoarsenTarget (configuration.get<
int>(
"coarsenTarget"));
547 if (configuration.hasKey (
"accumulationMode"))
549 std::string mode = ToLower()(configuration.get<std::string>(
"accumulationMode"));
552 else if ( mode ==
"atonce" )
554 else if ( mode ==
"successive")
557 DUNE_THROW(InvalidSolverFactoryConfiguration,
"Parameter accumulationMode does not allow value "
561 if (configuration.hasKey (
"prolongationDampingFactor"))
562 criterion.setProlongationDampingFactor (configuration.get<
double>(
"prolongationDampingFactor"));
564 if (configuration.hasKey(
"defaultAggregationSizeMode"))
566 auto mode = ToLower()(configuration.get<std::string>(
"defaultAggregationSizeMode"));
567 auto dim = configuration.get<std::size_t>(
"defaultAggregationDimension");
568 std::size_t maxDistance = 2;
569 if (configuration.hasKey(
"MaxAggregateDistance"))
570 maxDistance = configuration.get<std::size_t>(
"maxAggregateDistance");
571 if (mode ==
"isotropic")
572 criterion.setDefaultValuesIsotropic(dim, maxDistance);
573 else if(mode ==
"anisotropic")
574 criterion.setDefaultValuesAnisotropic(dim, maxDistance);
576 DUNE_THROW(InvalidSolverFactoryConfiguration,
"Parameter accumulationMode does not allow value "
580 if (configuration.hasKey(
"maxAggregateDistance"))
581 criterion.setMaxDistance(configuration.get<std::size_t>(
"maxAggregateDistance"));
583 if (configuration.hasKey(
"minAggregateSize"))
584 criterion.setMinAggregateSize(configuration.get<std::size_t>(
"minAggregateSize"));
586 if (configuration.hasKey(
"maxAggregateSize"))
587 criterion.setMaxAggregateSize(configuration.get<std::size_t>(
"maxAggregateSize"));
589 if (configuration.hasKey(
"maxAggregateConnectivity"))
590 criterion.setMaxConnectivity(configuration.get<std::size_t>(
"maxAggregateConnectivity"));
592 if (configuration.hasKey (
"alpha"))
593 criterion.setAlpha (configuration.get<
double> (
"alpha"));
595 if (configuration.hasKey (
"beta"))
596 criterion.setBeta (configuration.get<
double> (
"beta"));
598 if (configuration.hasKey (
"gamma"))
599 criterion.setGamma (configuration.get<std::size_t> (
"gamma"));
600 gamma_ = criterion.getGamma();
602 if (configuration.hasKey (
"additive"))
603 criterion.setAdditive (configuration.get<
bool>(
"additive"));
604 additive = criterion.getAdditive();
606 if (configuration.hasKey (
"preSteps"))
607 criterion.setNoPreSmoothSteps (configuration.get<std::size_t> (
"preSteps"));
608 preSteps_ = criterion.getNoPreSmoothSteps ();
610 if (configuration.hasKey (
"postSteps"))
611 criterion.setNoPostSmoothSteps (configuration.get<std::size_t> (
"preSteps"));
612 postSteps_ = criterion.getNoPostSmoothSteps ();
614 verbosity_ = configuration.get(
"verbosity", 0);
615 criterion.setDebugLevel (verbosity_);
617 createHierarchies(criterion, matrixptr, pinfo);
620 template <
class Matrix,
628 #if DISABLE_AMG_DIRECTSOLVER
630 #elif HAVE_SUITESPARSE_UMFPACK
638 template <
class M, SolverType>
644 DUNE_THROW(NotImplemented,
"DirectSolver not selected");
647 static std::string
name () {
return "None"; }
649 #if HAVE_SUITESPARSE_UMFPACK
654 static type*
create(
const M&
mat,
bool verbose,
bool reusevector )
656 return new type(
mat, verbose, reusevector );
658 static std::string
name () {
return "UMFPack"; }
668 return new type(
mat, verbose, reusevector );
670 static std::string
name () {
return "SuperLU"; }
685 template<
class M,
class X,
class S,
class PI,
class A>
687 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
688 const std::shared_ptr<const Operator>& matrixptr,
692 matrices_ = std::make_shared<OperatorHierarchy>(
693 std::const_pointer_cast<Operator>(matrixptr),
694 stackobject_to_shared_ptr(
const_cast<PI&
>(pinfo)));
696 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
699 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
705 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
706 && ( ! matrices_->redistributeInformation().back().isSetup() ||
707 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
710 SmootherArgs sargs(smootherArgs_);
711 sargs.iterations = 1;
714 cargs.setArgs(sargs);
715 if(matrices_->redistributeInformation().back().isSetup()) {
717 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
718 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
720 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
721 cargs.setComm(*matrices_->parallelInformation().coarsest());
725 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
727 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
730 if( SolverSelector::isDirectSolver &&
731 (std::is_same<ParallelInformation,SequentialInformation>::value
732 || matrices_->parallelInformation().coarsest()->communicator().size()==1
733 || (matrices_->parallelInformation().coarsest().isRedistributed()
734 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
735 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
738 if(matrices_->parallelInformation().coarsest().isRedistributed())
740 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
743 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
750 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(),
false,
false));
752 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
753 std::cout<<
"Using a direct coarse solver (" << SolverSelector::name() <<
")" << std::endl;
757 if(matrices_->parallelInformation().coarsest().isRedistributed())
759 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
764 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
766 *coarseSmoother_, 1E-2, 1000, 0));
771 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
773 *coarseSmoother_, 1E-2, 1000, 0));
792 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
793 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
794 <<
"(including coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
798 template<
class M,
class X,
class S,
class PI,
class A>
805 typedef typename M::matrix_type
Matrix;
812 const Matrix&
mat=matrices_->matrices().finest()->getmat();
813 for(RowIter row=
mat.begin(); row!=
mat.end(); ++row) {
814 bool isDirichlet =
true;
815 bool hasDiagonal =
false;
817 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
818 if(row.index()==
col.index()) {
826 if(isDirichlet && hasDiagonal)
828 auto&& xEntry = Impl::asVector(x[row.index()]);
829 auto&& bEntry = Impl::asVector(b[row.index()]);
830 Impl::asMatrix(diagonal).solve(xEntry, bEntry);
834 if(smoothers_->levels()>0)
835 smoothers_->finest()->pre(x,b);
838 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
839 rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
840 lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
841 update_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
842 matrices_->coarsenVector(*rhs_);
843 matrices_->coarsenVector(*lhs_);
844 matrices_->coarsenVector(*update_);
850 Iterator coarsest = smoothers_->
coarsest();
851 Iterator smoother = smoothers_->finest();
852 RIterator rhs = rhs_->finest();
853 DIterator lhs = lhs_->finest();
854 if(smoothers_->levels()>1) {
856 assert(lhs_->levels()==rhs_->levels());
857 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
858 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
860 if(smoother!=coarsest)
861 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
862 smoother->pre(*lhs,*rhs);
863 smoother->pre(*lhs,*rhs);
873 template<
class M,
class X,
class S,
class PI,
class A>
876 return matrices_->levels();
878 template<
class M,
class X,
class S,
class PI,
class A>
881 return matrices_->maxlevels();
885 template<
class M,
class X,
class S,
class PI,
class A>
888 LevelContext levelContext;
896 initIteratorsWithFineLevel(levelContext);
899 *levelContext.lhs = v;
900 *levelContext.rhs = d;
901 *levelContext.update=0;
902 levelContext.level=0;
906 if(postSteps_==0||matrices_->maxlevels()==1)
907 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
909 v=*levelContext.update;
914 template<
class M,
class X,
class S,
class PI,
class A>
917 levelContext.smoother = smoothers_->finest();
918 levelContext.matrix = matrices_->matrices().finest();
919 levelContext.pinfo = matrices_->parallelInformation().finest();
920 levelContext.redist =
921 matrices_->redistributeInformation().begin();
922 levelContext.aggregates = matrices_->aggregatesMaps().begin();
923 levelContext.lhs = lhs_->finest();
924 levelContext.update = update_->finest();
925 levelContext.rhs = rhs_->finest();
928 template<
class M,
class X,
class S,
class PI,
class A>
930 ::moveToCoarseLevel(LevelContext& levelContext)
933 bool processNextLevel=
true;
935 if(levelContext.redist->isSetup()) {
936 levelContext.redist->redistribute(
static_cast<const Range&
>(*levelContext.rhs),
937 levelContext.rhs.getRedistributed());
938 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
939 if(processNextLevel) {
942 ++levelContext.pinfo;
945 static_cast<const Range&
>(fineRhs.getRedistributed()),
946 *levelContext.pinfo);
951 ++levelContext.pinfo;
954 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
955 *levelContext.pinfo);
958 if(processNextLevel) {
961 ++levelContext.update;
962 ++levelContext.matrix;
963 ++levelContext.level;
964 ++levelContext.redist;
966 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
968 ++levelContext.smoother;
969 ++levelContext.aggregates;
972 *levelContext.update=0;
974 return processNextLevel;
977 template<
class M,
class X,
class S,
class PI,
class A>
979 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
981 if(processNextLevel) {
982 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
984 --levelContext.smoother;
985 --levelContext.aggregates;
987 --levelContext.redist;
988 --levelContext.level;
990 --levelContext.matrix;
994 --levelContext.pinfo;
996 if(levelContext.redist->isSetup()) {
998 levelContext.lhs.getRedistributed()=0;
1000 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1001 levelContext.lhs.getRedistributed(),
1002 matrices_->getProlongationDampingFactor(),
1003 *levelContext.pinfo, *levelContext.redist);
1005 *levelContext.lhs=0;
1007 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1008 matrices_->getProlongationDampingFactor(),
1009 *levelContext.pinfo);
1013 if(processNextLevel) {
1014 --levelContext.update;
1018 *levelContext.update += *levelContext.lhs;
1021 template<
class M,
class X,
class S,
class PI,
class A>
1027 template<
class M,
class X,
class S,
class PI,
class A>
1029 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
1033 if(levelContext.redist->isSetup()) {
1034 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
1035 if(levelContext.rhs.getRedistributed().size()>0) {
1037 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
1038 levelContext.rhs.getRedistributed());
1039 solver_->apply(levelContext.update.getRedistributed(),
1040 levelContext.rhs.getRedistributed(), res);
1042 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
1043 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
1045 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
1046 solver_->apply(*levelContext.update, *levelContext.rhs, res);
1050 coarsesolverconverged =
false;
1055 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1056 bool processNextLevel = moveToCoarseLevel(levelContext);
1058 if(processNextLevel) {
1060 for(std::size_t i=0; i<gamma_; i++)
1064 moveToFineLevel(levelContext, processNextLevel);
1069 if(levelContext.matrix == matrices_->matrices().finest()) {
1070 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
1071 if(!coarsesolverconverged)
1072 DUNE_THROW(MathError,
"Coarse solver did not converge");
1080 template<
class M,
class X,
class S,
class PI,
class A>
1081 void AMG<M,X,S,PI,A>::additiveMgc(){
1084 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
1087 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
1092 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
1098 lhs = lhs_->finest();
1101 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
1104 smoother->apply(*lhs, *rhs);
1108 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1109 InverseOperatorResult res;
1110 pinfo->copyOwnerToAll(*rhs, *rhs);
1111 solver_->apply(*lhs, *rhs, res);
1114 DUNE_THROW(MathError,
"Coarse solver did not converge");
1130 template<
class M,
class X,
class S,
class PI,
class A>
1136 Iterator coarsest = smoothers_->
coarsest();
1137 Iterator smoother = smoothers_->finest();
1138 DIterator lhs = lhs_->finest();
1139 if(smoothers_->levels()>0) {
1140 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
1141 smoother->post(*lhs);
1142 if(smoother!=coarsest)
1143 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
1144 smoother->post(*lhs);
1145 smoother->post(*lhs);
1152 template<
class M,
class X,
class S,
class PI,
class A>
1156 matrices_->getCoarsestAggregatesOnFinest(cont);
1166 std::shared_ptr<Dune::Preconditioner<typename OP::element_type::domain_type, typename OP::element_type::range_type> >
1167 makeAMG(
const OP& op,
const std::string& smoother,
const Dune::ParameterTree& config)
const
1169 DUNE_THROW(Dune::Exception,
"Operator type not supported by AMG");
1172 template<
class M,
class X,
class Y>
1173 std::shared_ptr<Dune::Preconditioner<X,Y> >
1175 const Dune::ParameterTree& config)
const
1179 if(smoother ==
"ssor")
1180 return std::make_shared<Amg::AMG<OP, X, SeqSSOR<M,X,Y>>>(op, config);
1181 if(smoother ==
"sor")
1182 return std::make_shared<Amg::AMG<OP, X, SeqSOR<M,X,Y>>>(op, config);
1183 if(smoother ==
"jac")
1184 return std::make_shared<Amg::AMG<OP, X, SeqJac<M,X,Y>>>(op, config);
1185 if(smoother ==
"gs")
1186 return std::make_shared<Amg::AMG<OP, X, SeqGS<M,X,Y>>>(op, config);
1187 if(smoother ==
"ilu")
1188 return std::make_shared<Amg::AMG<OP, X, SeqILU<M,X,Y>>>(op, config);
1190 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1193 template<
class M,
class X,
class Y,
class C>
1194 std::shared_ptr<Dune::Preconditioner<X,Y> >
1196 const Dune::ParameterTree& config)
const
1200 auto cop = std::static_pointer_cast<const OP>(op);
1202 if(smoother ==
"ssor")
1203 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1204 if(smoother ==
"sor")
1205 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1206 if(smoother ==
"jac")
1207 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqJac<M,X,Y>>,C>>(cop, config, op->getCommunication());
1208 if(smoother ==
"gs")
1209 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqGS<M,X,Y>>,C>>(cop, config, op->getCommunication());
1210 if(smoother ==
"ilu")
1211 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqILU<M,X,Y>>,C>>(cop, config, op->getCommunication());
1213 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1216 template<
class M,
class X,
class Y,
class C>
1217 std::shared_ptr<Dune::Preconditioner<X,Y> >
1219 const Dune::ParameterTree& config)
const
1223 if(smoother ==
"ssor")
1224 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1225 if(smoother ==
"sor")
1226 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1227 if(smoother ==
"jac")
1228 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqJac<M,X,Y>>,C>>(op, config, op->getCommunication());
1229 if(smoother ==
"gs")
1230 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqGS<M,X,Y>>,C>>(op, config, op->getCommunication());
1231 if(smoother ==
"ilu")
1232 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqILU<M,X,Y>>,C>>(op, config, op->getCommunication());
1234 DUNE_THROW(Dune::Exception,
"Unknown smoother for AMG");
1237 template<
typename TL,
typename OP>
1238 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1239 typename Dune::TypeListElement<2, TL>::type>>
1240 operator() (TL tl,
const std::shared_ptr<OP>& op,
const Dune::ParameterTree& config,
1243 using field_type =
typename OP::matrix_type::field_type;
1244 using real_type =
typename FieldTraits<field_type>::real_type;
1245 if (!std::is_convertible<field_type, real_type>())
1247 className<field_type>() <<
1248 ") to be convertible to its real_type (" <<
1249 className<real_type>() <<
1251 using D =
typename Dune::TypeListElement<1, decltype(tl)>::type;
1252 using R =
typename Dune::TypeListElement<2, decltype(tl)>::type;
1253 std::shared_ptr<Preconditioner<D,R>> amg;
1254 std::string smoother = config.get(
"smoother",
"ssor");
1255 return makeAMG(op, smoother, config);
1258 template<
typename TL,
typename OP>
1259 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1260 typename Dune::TypeListElement<2, TL>::type>>
1261 operator() (TL ,
const std::shared_ptr<OP>& ,
const Dune::ParameterTree& ,
1264 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:349
Matrix & mat
Definition: matrixmatrix.hh:345
AMG(const AMG &amg)
Copy constructor.
Definition: amg.hh:390
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:799
static std::string name()
Definition: amg.hh:678
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:301
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:305
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:1167
static std::string name()
Definition: amg.hh:670
Matrix ::field_type field_type
Definition: amg.hh:624
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:1022
X Domain
The domain type.
Definition: amg.hh:85
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:679
AMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:404
AMG(std::shared_ptr< const Operator > fineOperator, const ParameterTree &configuration, const ParallelInformation &pinfo=ParallelInformation())
Constructor an AMG via ParameterTree.
Definition: amg.hh:450
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:285
SolverType
Definition: amg.hh:625
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:293
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:98
Solver< Matrix, solver > SelectedSolver
Definition: amg.hh:675
std::string operator()(const std::string &str)
Definition: amg.hh:376
S Smoother
The type of the smoother.
Definition: amg.hh:95
static std::string name()
Definition: amg.hh:647
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:277
M Operator
The matrix operator type.
Definition: amg.hh:71
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:281
SelectedSolver ::type DirectSolver
Definition: amg.hh:676
SuperLU< M > type
Definition: amg.hh:665
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:289
X Range
The range type.
Definition: amg.hh:87
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:1154
std::size_t levels()
Definition: amg.hh:874
InverseOperator< Vector, Vector > type
Definition: amg.hh:641
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:1174
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:297
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:1218
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:42
static constexpr SolverType solver
Definition: amg.hh:627
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:50
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:1195
static constexpr bool isDirectSolver
Definition: amg.hh:677
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:219
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:381
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:82
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:89
void post(Domain &x)
Clean up.
Definition: amg.hh:1131
std::size_t maxlevels()
Definition: amg.hh:879
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:1240
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:642
std::size_t level
The level index.
Definition: amg.hh:309
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:426
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:886
Smoother SmootherType
Definition: amg.hh:273
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:80
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amg.hh:192
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:666
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:78
@ none
Definition: amg.hh:625
@ umfpack
Definition: amg.hh:625
@ superlu
Definition: amg.hh:625
@ atOnceAccu
Accumulate data to one process at once.
Definition: parameters.hh:242
@ noAccu
No data accumulution.
Definition: parameters.hh:236
@ successiveAccu
Successively accumulate to fewer processes.
Definition: parameters.hh:246
Definition: allocator.hh:9
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
Definition: matrixutils.hh:25
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:62
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:461
Definition: aggregates.hh:478
Definition: aggregates.hh:494
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:517
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:537
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:138
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:31
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:63
An overlapping Schwarz operator.
Definition: schwarz.hh:76
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:214
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:217
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:59
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:281
All parameters for AMG.
Definition: parameters.hh:391
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
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:30
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioner.hh:37
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:50
Statistics about the application of an inverse operator.
Definition: solver.hh:46
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Categories for the solvers.
Definition: solvercategory.hh:20
Category
Definition: solvercategory.hh:21
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:32
Definition: solvercategory.hh:52
Definition: solverregistry.hh:75
Definition: solvertype.hh:14
SuperLu Solver.
Definition: superlu.hh:269
Definition: umfpack.hh:47
The UMFPack direct sparse solver.
Definition: umfpack.hh:213