dune-istl  2.8.0
amg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_AMG_AMG_HH
4 #define DUNE_AMG_AMG_HH
5 
6 #include <memory>
7 #include <sstream>
8 #include <dune/common/exceptions.hh>
12 #include <dune/istl/solvers.hh>
14 #include <dune/istl/superlu.hh>
15 #include <dune/istl/umfpack.hh>
16 #include <dune/istl/solvertype.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>
22 
23 namespace Dune
24 {
25  namespace Amg
26  {
44  template<class M, class X, class S, class P, class K, class A>
45  class KAMG;
46 
47  template<class T>
48  class KAmgTwoGrid;
49 
60  template<class M, class X, class S, class PI=SequentialInformation,
61  class A=std::allocator<X> >
62  class AMG : public Preconditioner<X,X>
63  {
64  template<class M1, class X1, class S1, class P1, class K1, class A1>
65  friend class KAMG;
66 
67  friend class KAmgTwoGrid<AMG>;
68 
69  public:
71  typedef M Operator;
78  typedef PI ParallelInformation;
83 
85  typedef X Domain;
87  typedef X Range;
95  typedef S Smoother;
96 
99 
109  AMG(OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
110  const SmootherArgs& smootherArgs, const Parameters& parms);
111 
123  template<class C>
124  AMG(const Operator& fineOperator, const C& criterion,
125  const SmootherArgs& smootherArgs=SmootherArgs(),
127 
178  AMG(std::shared_ptr<const Operator> fineOperator, const ParameterTree& configuration, const ParallelInformation& pinfo=ParallelInformation());
179 
183  AMG(const AMG& amg);
184 
186  void pre(Domain& x, Range& b);
187 
189  void apply(Domain& v, const Range& d);
190 
193  {
194  return category_;
195  }
196 
198  void post(Domain& x);
199 
204  template<class A1>
205  void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
206 
207  std::size_t levels();
208 
209  std::size_t maxlevels();
210 
220  {
221  matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
222  }
223 
229 
230  private:
231  /*
232  * @brief Helper function to create hierarchies with parameter tree.
233  *
234  * Will create the coarsen criterion with the norm and create the
235  * Hierarchies
236  * \tparam Norm Type of the norm to use.
237  */
238  template<class Norm>
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());
243  template<class Norm>
244  void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
245  const PI& pinfo, const Norm&,
246  const ParameterTree& configuration,
247  std::false_type);
252  template<class C>
253  void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
254  const PI& pinfo, const ParameterTree& configuration);
261  template<class C>
262  void createHierarchies(C& criterion,
263  const std::shared_ptr<const Operator>& matrixptr,
264  const PI& pinfo);
271  struct LevelContext
272  {
289  typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
293  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
309  std::size_t level;
310  };
311 
312 
317  void mgc(LevelContext& levelContext);
318 
319  void additiveMgc();
320 
327  void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
328 
333  bool moveToCoarseLevel(LevelContext& levelContext);
334 
339  void initIteratorsWithFineLevel(LevelContext& levelContext);
340 
342  std::shared_ptr<OperatorHierarchy> matrices_;
344  SmootherArgs smootherArgs_;
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_;
360  std::size_t gamma_;
362  std::size_t preSteps_;
364  std::size_t postSteps_;
365  bool buildHierarchy_;
366  bool additive;
367  bool coarsesolverconverged;
368  std::shared_ptr<Smoother> coarseSmoother_;
370  SolverCategory::Category category_;
372  std::size_t verbosity_;
373 
374  struct ToLower
375  {
376  std::string operator()(const std::string& str)
377  {
378  std::stringstream retval;
379  std::ostream_iterator<char> out(retval);
380  std::transform(str.begin(), str.end(), out,
381  [](char c){
382  return std::tolower(c, std::locale::classic());
383  });
384  return retval.str();
385  }
386  };
387  };
388 
389  template<class M, class X, class S, class PI, class A>
390  inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
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_)
401  {}
402 
403  template<class M, class X, class S, class PI, class A>
405  const SmootherArgs& smootherArgs,
406  const Parameters& parms)
407  : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
408  smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
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),
413  coarseSmoother_(),
414 // #warning should category be retrieved from matrices?
415  category_(SolverCategory::category(*smoothers_->coarsest())),
416  verbosity_(parms.debugLevel())
417  {
418  assert(matrices_->isBuilt());
419 
420  // build the necessary smoother hierarchies
421  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
422  }
423 
424  template<class M, class X, class S, class PI, class A>
425  template<class C>
427  const C& criterion,
428  const SmootherArgs& smootherArgs,
429  const PI& pinfo)
430  : smootherArgs_(smootherArgs),
431  smoothers_(new Hierarchy<Smoother,A>), solver_(),
432  rhs_(), lhs_(), update_(), scalarProduct_(),
433  gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
434  postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
435  additive(criterion.getAdditive()), coarsesolverconverged(true),
436  coarseSmoother_(),
437  category_(SolverCategory::category(pinfo)),
438  verbosity_(criterion.debugLevel())
439  {
441  DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
442  // TODO: reestablish compile time checks.
443  //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
444  // "Matrix and Solver must match in terms of category!");
445  auto matrixptr = stackobject_to_shared_ptr(matrix);
446  createHierarchies(criterion, matrixptr, pinfo);
447  }
448 
449  template<class M, class X, class S, class PI, class A>
450  AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
451  const ParameterTree& configuration,
452  const ParallelInformation& pinfo) :
453  smoothers_(new Hierarchy<Smoother,A>),
454  solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
455  coarsesolverconverged(true), coarseSmoother_(),
456  category_(SolverCategory::category(pinfo))
457  {
458 
459  if (configuration.hasKey ("smootherIterations"))
460  smootherArgs_.iterations = configuration.get<int>("smootherIterations");
461 
462  if (configuration.hasKey ("smootherRelaxation"))
463  smootherArgs_.relaxationFactor = configuration.get<typename SmootherArgs::RelaxationFactor>("smootherRelaxation");
464 
465  auto normName = ToLower()(configuration.get("strengthMeasure", "diagonal"));
466  auto index = configuration.get<int>("diagonalRowIndex", 0);
467 
468  if ( normName == "diagonal")
469  {
470  using field_type = typename M::field_type;
471  using real_type = typename FieldTraits<field_type>::real_type;
472  std::is_convertible<field_type, real_type> compiles;
473 
474  switch (index)
475  {
476  case 0:
477  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<0>(), configuration, compiles);
478  break;
479  case 1:
480  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<1>(), configuration, compiles);
481  break;
482  case 2:
483  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<2>(), configuration, compiles);
484  break;
485  case 3:
486  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<3>(), configuration, compiles);
487  break;
488  case 4:
489  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<4>(), configuration, compiles);
490  break;
491  default:
492  DUNE_THROW(InvalidStateException, "Currently strengthIndex>4 is not supported.");
493  }
494  }
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);
501  else
502  DUNE_THROW(Dune::NotImplemented, "Wrong config file: strengthMeasure "<<normName<<" is not supported by AMG");
503  }
504 
505  template<class M, class X, class S, class PI, class A>
506  template<class Norm>
507  void AMG<M,X,S,PI,A>::createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const Norm&, const ParameterTree& configuration, std::false_type)
508  {
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>() << ".");
512  }
513 
514  template<class M, class X, class S, class PI, class A>
515  template<class Norm>
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)
517  {
518  if (configuration.get<bool>("criterionSymmetric", true))
519  {
520  using Criterion = Dune::Amg::CoarsenCriterion<
522  Criterion criterion;
523  createHierarchies(criterion, matrixptr, pinfo, configuration);
524  }
525  else
526  {
527  using Criterion = Dune::Amg::CoarsenCriterion<
529  Criterion criterion;
530  createHierarchies(criterion, matrixptr, pinfo, configuration);
531  }
532  }
533 
534  template<class M, class X, class S, class PI, class A>
535  template<class C>
536  void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const ParameterTree& configuration)
537  {
538  if (configuration.hasKey ("maxLevel"))
539  criterion.setMaxLevel(configuration.get<int>("maxLevel"));
540 
541  if (configuration.hasKey ("minCoarseningRate"))
542  criterion.setMinCoarsenRate(configuration.get<int>("minCoarseningRate"));
543 
544  if (configuration.hasKey ("coarsenTarget"))
545  criterion.setCoarsenTarget (configuration.get<int>("coarsenTarget"));
546 
547  if (configuration.hasKey ("accumulationMode"))
548  {
549  std::string mode = ToLower()(configuration.get<std::string>("accumulationMode"));
550  if ( mode == "none")
551  criterion.setAccumulate(AccumulationMode::noAccu);
552  else if ( mode == "atonce" )
553  criterion.setAccumulate(AccumulationMode::atOnceAccu);
554  else if ( mode == "successive")
555  criterion.setCoarsenTarget (AccumulationMode::successiveAccu);
556  else
557  DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
558  << mode <<".");
559  }
560 
561  if (configuration.hasKey ("prolongationDampingFactor"))
562  criterion.setProlongationDampingFactor (configuration.get<double>("prolongationDampingFactor"));
563 
564  if (configuration.hasKey("defaultAggregationSizeMode"))
565  {
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);
575  else
576  DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
577  << mode <<".");
578  }
579 
580  if (configuration.hasKey("maxAggregateDistance"))
581  criterion.setMaxDistance(configuration.get<std::size_t>("maxAggregateDistance"));
582 
583  if (configuration.hasKey("minAggregateSize"))
584  criterion.setMinAggregateSize(configuration.get<std::size_t>("minAggregateSize"));
585 
586  if (configuration.hasKey("maxAggregateSize"))
587  criterion.setMaxAggregateSize(configuration.get<std::size_t>("maxAggregateSize"));
588 
589  if (configuration.hasKey("maxAggregateConnectivity"))
590  criterion.setMaxConnectivity(configuration.get<std::size_t>("maxAggregateConnectivity"));
591 
592  if (configuration.hasKey ("alpha"))
593  criterion.setAlpha (configuration.get<double> ("alpha"));
594 
595  if (configuration.hasKey ("beta"))
596  criterion.setBeta (configuration.get<double> ("beta"));
597 
598  if (configuration.hasKey ("gamma"))
599  criterion.setGamma (configuration.get<std::size_t> ("gamma"));
600  gamma_ = criterion.getGamma();
601 
602  if (configuration.hasKey ("additive"))
603  criterion.setAdditive (configuration.get<bool>("additive"));
604  additive = criterion.getAdditive();
605 
606  if (configuration.hasKey ("preSteps"))
607  criterion.setNoPreSmoothSteps (configuration.get<std::size_t> ("preSteps"));
608  preSteps_ = criterion.getNoPreSmoothSteps ();
609 
610  if (configuration.hasKey ("postSteps"))
611  criterion.setNoPostSmoothSteps (configuration.get<std::size_t> ("preSteps"));
612  postSteps_ = criterion.getNoPostSmoothSteps ();
613 
614  verbosity_ = configuration.get("verbosity", 0);
615  criterion.setDebugLevel (verbosity_);
616 
617  createHierarchies(criterion, matrixptr, pinfo);
618  }
619 
620  template <class Matrix,
621  class Vector>
623  {
626 
627  static constexpr SolverType solver =
628 #if DISABLE_AMG_DIRECTSOLVER
629  none;
630 #elif HAVE_SUITESPARSE_UMFPACK
632 #elif HAVE_SUPERLU
633  superlu ;
634 #else
635  none;
636 #endif
637 
638  template <class M, SolverType>
639  struct Solver
640  {
642  static type* create(const M& mat, bool verbose, bool reusevector )
643  {
644  DUNE_THROW(NotImplemented,"DirectSolver not selected");
645  return nullptr;
646  }
647  static std::string name () { return "None"; }
648  };
649 #if HAVE_SUITESPARSE_UMFPACK
650  template <class M>
651  struct Solver< M, umfpack >
652  {
653  typedef UMFPack< M > type;
654  static type* create(const M& mat, bool verbose, bool reusevector )
655  {
656  return new type(mat, verbose, reusevector );
657  }
658  static std::string name () { return "UMFPack"; }
659  };
660 #endif
661 #if HAVE_SUPERLU
662  template <class M>
663  struct Solver< M, superlu >
664  {
666  static type* create(const M& mat, bool verbose, bool reusevector )
667  {
668  return new type(mat, verbose, reusevector );
669  }
670  static std::string name () { return "SuperLU"; }
671  };
672 #endif
673 
674  // define direct solver type to be used
677  static constexpr bool isDirectSolver = solver != none;
678  static std::string name() { return SelectedSolver :: name (); }
679  static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
680  {
681  return SelectedSolver :: create( mat, verbose, reusevector );
682  }
683  };
684 
685  template<class M, class X, class S, class PI, class A>
686  template<class C>
687  void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
688  const std::shared_ptr<const Operator>& matrixptr,
689  const PI& pinfo)
690  {
691  Timer watch;
692  matrices_ = std::make_shared<OperatorHierarchy>(
693  std::const_pointer_cast<Operator>(matrixptr),
694  stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
695 
696  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
697 
698  // build the necessary smoother hierarchies
699  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
700 
701  // test whether we should solve on the coarse level. That is the case if we
702  // have that level and if there was a redistribution on this level then our
703  // communicator has to be valid (size()>0) as the smoother might try to communicate
704  // in the constructor.
705  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
706  && ( ! matrices_->redistributeInformation().back().isSetup() ||
707  matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
708  {
709  // We have the carsest level. Create the coarse Solver
710  SmootherArgs sargs(smootherArgs_);
711  sargs.iterations = 1;
712 
714  cargs.setArgs(sargs);
715  if(matrices_->redistributeInformation().back().isSetup()) {
716  // Solve on the redistributed partitioning
717  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
718  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
719  }else{
720  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
721  cargs.setComm(*matrices_->parallelInformation().coarsest());
722  }
723 
724  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
725  scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
726 
727  typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
728 
729  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
730  if( SolverSelector::isDirectSolver &&
731  (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
732  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
733  || (matrices_->parallelInformation().coarsest().isRedistributed()
734  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
735  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
736  )
737  { // redistribute and 1 proc
738  if(matrices_->parallelInformation().coarsest().isRedistributed())
739  {
740  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
741  {
742  // We are still participating on this level
743  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
744  }
745  else
746  solver_.reset();
747  }
748  else
749  {
750  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
751  }
752  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
753  std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
754  }
755  else
756  {
757  if(matrices_->parallelInformation().coarsest().isRedistributed())
758  {
759  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
760  // We are still participating on this level
761 
762  // we have to allocate these types using the rebound allocator
763  // in order to ensure that we fulfill the alignment requirements
764  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
765  *scalarProduct_,
766  *coarseSmoother_, 1E-2, 1000, 0));
767  else
768  solver_.reset();
769  }else
770  {
771  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
772  *scalarProduct_,
773  *coarseSmoother_, 1E-2, 1000, 0));
774  // // we have to allocate these types using the rebound allocator
775  // // in order to ensure that we fulfill the alignment requirements
776  // using Alloc = typename std::allocator_traits<A>::template rebind_alloc<BiCGSTABSolver<X>>;
777  // Alloc alloc;
778  // auto p = alloc.allocate(1);
779  // std::allocator_traits<Alloc>::construct(alloc, p,
780  // const_cast<M&>(*matrices_->matrices().coarsest()),
781  // *scalarProduct_,
782  // *coarseSmoother_, 1E-2, 1000, 0);
783  // solver_.reset(p,[](BiCGSTABSolver<X>* p){
784  // Alloc alloc;
785  // std::allocator_traits<Alloc>::destroy(alloc, p);
786  // alloc.deallocate(p,1);
787  // });
788  }
789  }
790  }
791 
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;
795  }
796 
797 
798  template<class M, class X, class S, class PI, class A>
800  {
801  // Detect Matrix rows where all offdiagonal entries are
802  // zero and set x such that A_dd*x_d=b_d
803  // Thus users can be more careless when setting up their linear
804  // systems.
805  typedef typename M::matrix_type Matrix;
806  typedef typename Matrix::ConstRowIterator RowIter;
807  typedef typename Matrix::ConstColIterator ColIter;
808  typedef typename Matrix::block_type Block;
809  Block zero;
810  zero=typename Matrix::field_type();
811 
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;
816  Block diagonal{};
817  for(ColIter col=row->begin(); col!=row->end(); ++col) {
818  if(row.index()==col.index()) {
819  diagonal = *col;
820  hasDiagonal = true;
821  }else{
822  if(*col!=zero)
823  isDirichlet = false;
824  }
825  }
826  if(isDirichlet && hasDiagonal)
827  {
828  auto&& xEntry = Impl::asVector(x[row.index()]);
829  auto&& bEntry = Impl::asVector(b[row.index()]);
830  Impl::asMatrix(diagonal).solve(xEntry, bEntry);
831  }
832  }
833 
834  if(smoothers_->levels()>0)
835  smoothers_->finest()->pre(x,b);
836  else
837  // No smoother to make x consistent! Do it by hand
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_);
845 
846  // Preprocess all smoothers
847  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
848  typedef typename Hierarchy<Range,A>::Iterator RIterator;
849  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
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) {
855 
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());
859 
860  if(smoother!=coarsest)
861  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
862  smoother->pre(*lhs,*rhs);
863  smoother->pre(*lhs,*rhs);
864  }
865 
866 
867  // The preconditioner might change x and b. So we have to
868  // copy the changes to the original vectors.
869  x = *lhs_->finest();
870  b = *rhs_->finest();
871 
872  }
873  template<class M, class X, class S, class PI, class A>
875  {
876  return matrices_->levels();
877  }
878  template<class M, class X, class S, class PI, class A>
880  {
881  return matrices_->maxlevels();
882  }
883 
885  template<class M, class X, class S, class PI, class A>
887  {
888  LevelContext levelContext;
889 
890  if(additive) {
891  *(rhs_->finest())=d;
892  additiveMgc();
893  v=*lhs_->finest();
894  }else{
895  // Init all iterators for the current level
896  initIteratorsWithFineLevel(levelContext);
897 
898 
899  *levelContext.lhs = v;
900  *levelContext.rhs = d;
901  *levelContext.update=0;
902  levelContext.level=0;
903 
904  mgc(levelContext);
905 
906  if(postSteps_==0||matrices_->maxlevels()==1)
907  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
908 
909  v=*levelContext.update;
910  }
911 
912  }
913 
914  template<class M, class X, class S, class PI, class A>
915  void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
916  {
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();
926  }
927 
928  template<class M, class X, class S, class PI, class A>
929  bool AMG<M,X,S,PI,A>
930  ::moveToCoarseLevel(LevelContext& levelContext)
931  {
932 
933  bool processNextLevel=true;
934 
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) {
940  //restrict defect to coarse level right hand side.
941  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
942  ++levelContext.pinfo;
944  ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
945  static_cast<const Range&>(fineRhs.getRedistributed()),
946  *levelContext.pinfo);
947  }
948  }else{
949  //restrict defect to coarse level right hand side.
950  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
951  ++levelContext.pinfo;
953  ::restrictVector(*(*levelContext.aggregates),
954  *levelContext.rhs, static_cast<const Range&>(*fineRhs),
955  *levelContext.pinfo);
956  }
957 
958  if(processNextLevel) {
959  // prepare coarse system
960  ++levelContext.lhs;
961  ++levelContext.update;
962  ++levelContext.matrix;
963  ++levelContext.level;
964  ++levelContext.redist;
965 
966  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
967  // next level is not the globally coarsest one
968  ++levelContext.smoother;
969  ++levelContext.aggregates;
970  }
971  // prepare the update on the next level
972  *levelContext.update=0;
973  }
974  return processNextLevel;
975  }
976 
977  template<class M, class X, class S, class PI, class A>
978  void AMG<M,X,S,PI,A>
979  ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
980  {
981  if(processNextLevel) {
982  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
983  // previous level is not the globally coarsest one
984  --levelContext.smoother;
985  --levelContext.aggregates;
986  }
987  --levelContext.redist;
988  --levelContext.level;
989  //prolongate and add the correction (update is in coarse left hand side)
990  --levelContext.matrix;
991 
992  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
993  --levelContext.lhs;
994  --levelContext.pinfo;
995  }
996  if(levelContext.redist->isSetup()) {
997  // Need to redistribute during prolongateVector
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);
1004  }else{
1005  *levelContext.lhs=0;
1007  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1008  matrices_->getProlongationDampingFactor(),
1009  *levelContext.pinfo);
1010  }
1011 
1012 
1013  if(processNextLevel) {
1014  --levelContext.update;
1015  --levelContext.rhs;
1016  }
1017 
1018  *levelContext.update += *levelContext.lhs;
1019  }
1020 
1021  template<class M, class X, class S, class PI, class A>
1023  {
1025  }
1026 
1027  template<class M, class X, class S, class PI, class A>
1028  void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
1029  if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
1030  // Solve directly
1032  res.converged=true; // If we do not compute this flag will not get updated
1033  if(levelContext.redist->isSetup()) {
1034  levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
1035  if(levelContext.rhs.getRedistributed().size()>0) {
1036  // We are still participating in the computation
1037  levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
1038  levelContext.rhs.getRedistributed());
1039  solver_->apply(levelContext.update.getRedistributed(),
1040  levelContext.rhs.getRedistributed(), res);
1041  }
1042  levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
1043  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
1044  }else{
1045  levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
1046  solver_->apply(*levelContext.update, *levelContext.rhs, res);
1047  }
1048 
1049  if (!res.converged)
1050  coarsesolverconverged = false;
1051  }else{
1052  // presmoothing
1053  presmooth(levelContext, preSteps_);
1054 
1055 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1056  bool processNextLevel = moveToCoarseLevel(levelContext);
1057 
1058  if(processNextLevel) {
1059  // next level
1060  for(std::size_t i=0; i<gamma_; i++)
1061  mgc(levelContext);
1062  }
1063 
1064  moveToFineLevel(levelContext, processNextLevel);
1065 #else
1066  *lhs=0;
1067 #endif
1068 
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");
1073  }
1074  // postsmoothing
1075  postsmooth(levelContext, postSteps_);
1076 
1077  }
1078  }
1079 
1080  template<class M, class X, class S, class PI, class A>
1081  void AMG<M,X,S,PI,A>::additiveMgc(){
1082 
1083  // restrict residual to all levels
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();
1088 
1089  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
1090  ++pinfo;
1092  ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
1093  }
1094 
1095  // pinfo is invalid, set to coarsest level
1096  //pinfo = matrices_->parallelInformation().coarsest
1097  // calculate correction for all levels
1098  lhs = lhs_->finest();
1099  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
1100 
1101  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
1102  // presmoothing
1103  *lhs=0;
1104  smoother->apply(*lhs, *rhs);
1105  }
1106 
1107  // Coarse level solve
1108 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1109  InverseOperatorResult res;
1110  pinfo->copyOwnerToAll(*rhs, *rhs);
1111  solver_->apply(*lhs, *rhs, res);
1112 
1113  if(!res.converged)
1114  DUNE_THROW(MathError, "Coarse solver did not converge");
1115 #else
1116  *lhs=0;
1117 #endif
1118  // Prologate and add up corrections from all levels
1119  --pinfo;
1120  --aggregates;
1121 
1122  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
1124  ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
1125  }
1126  }
1127 
1128 
1130  template<class M, class X, class S, class PI, class A>
1131  void AMG<M,X,S,PI,A>::post([[maybe_unused]] Domain& x)
1132  {
1133  // Postprocess all smoothers
1134  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
1135  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
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);
1146  }
1147  lhs_ = nullptr;
1148  update_ = nullptr;
1149  rhs_ = nullptr;
1150  }
1151 
1152  template<class M, class X, class S, class PI, class A>
1153  template<class A1>
1154  void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
1155  {
1156  matrices_->getCoarsestAggregatesOnFinest(cont);
1157  }
1158 
1159  } // end namespace Amg
1160 
1161  struct AMGCreator{
1162  template<class> struct isValidBlockType : std::false_type{};
1163  template<class T, int n, int m> struct isValidBlockType<FieldMatrix<T,n,m>> : std::true_type{};
1164 
1165  template<class OP>
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
1168  {
1169  DUNE_THROW(Dune::Exception, "Operator type not supported by AMG");
1170  }
1171 
1172  template<class M, class X, class Y>
1173  std::shared_ptr<Dune::Preconditioner<X,Y> >
1174  makeAMG(const std::shared_ptr<MatrixAdapter<M,X,Y>>& op, const std::string& smoother,
1175  const Dune::ParameterTree& config) const
1176  {
1177  using OP = MatrixAdapter<M,X,Y>;
1178 
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);
1189  else
1190  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1191  }
1192 
1193  template<class M, class X, class Y, class C>
1194  std::shared_ptr<Dune::Preconditioner<X,Y> >
1195  makeAMG(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1196  const Dune::ParameterTree& config) const
1197  {
1199 
1200  auto cop = std::static_pointer_cast<const OP>(op);
1201 
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());
1212  else
1213  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1214  }
1215 
1216  template<class M, class X, class Y, class C>
1217  std::shared_ptr<Dune::Preconditioner<X,Y> >
1218  makeAMG(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1219  const Dune::ParameterTree& config) const
1220  {
1222 
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());
1233  else
1234  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1235  }
1236 
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,
1241  std::enable_if_t<isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
1242  {
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>())
1246  DUNE_THROW(UnsupportedType, "AMG needs field_type(" <<
1247  className<field_type>() <<
1248  ") to be convertible to its real_type (" <<
1249  className<real_type>() <<
1250  ").");
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);
1256  }
1257 
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 /*tl*/, const std::shared_ptr<OP>& /*mat*/, const Dune::ParameterTree& /*config*/,
1262  std::enable_if_t<!isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
1263  {
1264  DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
1265  }
1266  };
1267 
1269 } // end namespace Dune
1270 
1271 #endif
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
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
Definition: amg.hh:623
Definition: amg.hh:1161
Definition: amg.hh:1162
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
Definition: pinfo.hh:26
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