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,
281 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
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"; }
678 static std::string
name() {
return SelectedSolver :: name (); }
681 return SelectedSolver :: create(
mat, verbose, reusevector );
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) {
941 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
942 ++levelContext.pinfo;
943 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
944 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
945 static_cast<const Range&
>(fineRhs.getRedistributed()),
946 *levelContext.pinfo);
950 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
951 ++levelContext.pinfo;
952 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
953 ::restrictVector(*(*levelContext.aggregates),
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;
999 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
1000 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1001 levelContext.lhs.getRedistributed(),
1002 matrices_->getProlongationDampingFactor(),
1003 *levelContext.pinfo, *levelContext.redist);
1005 *levelContext.lhs=0;
1006 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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();
1085 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
1086 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
1087 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
1089 for(
typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
1091 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
1092 ::restrictVector(*(*aggregates), *rhs,
static_cast<const Range&
>(*fineRhs), *pinfo);
1098 lhs = lhs_->finest();
1099 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->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");
1122 for(
typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
1123 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
1124 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
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.
#define DUNE_REGISTER_PRECONDITIONER(name,...)
Definition: solverregistry.hh:14
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 DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:679
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
static std::string name()
Definition: amg.hh:670
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 type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:642
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::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
std::string operator()(const std::string &str)
Definition: amg.hh:376
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
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
SuperLU< M > type
Definition: amg.hh:665
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:666
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< OverlappingSchwarzOperator< M, X, Y, C > > &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1195
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:297
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
static constexpr bool isDirectSolver
Definition: amg.hh:677
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:219
Matrix::field_type field_type
Definition: amg.hh:624
SelectedSolver::type DirectSolver
Definition: amg.hh:676
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
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< 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
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
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
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
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
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition: aggregates.hh:377
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
Iterator over the levels in the hierarchy.
Definition: hierarchy.hh:118
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
The default class for the smoother arguments.
Definition: smoother.hh:36
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
Abstract base class for all solvers.
Definition: solver.hh:97
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