dune-istl  2.9.0
amg.hh
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4 // vi: set et ts=4 sw=2 sts=2:
5 #ifndef DUNE_AMG_AMG_HH
6 #define DUNE_AMG_AMG_HH
7 
8 #include <memory>
9 #include <sstream>
10 #include <dune/common/exceptions.hh>
14 #include <dune/istl/solvers.hh>
16 #include <dune/istl/superlu.hh>
17 #include <dune/istl/umfpack.hh>
18 #include <dune/istl/solvertype.hh>
19 #include <dune/common/typetraits.hh>
20 #include <dune/common/exceptions.hh>
21 #include <dune/common/scalarvectorview.hh>
22 #include <dune/common/scalarmatrixview.hh>
23 #include <dune/common/parametertree.hh>
24 
25 namespace Dune
26 {
27  namespace Amg
28  {
46  template<class M, class X, class S, class P, class K, class A>
47  class KAMG;
48 
49  template<class T>
50  class KAmgTwoGrid;
51 
62  template<class M, class X, class S, class PI=SequentialInformation,
63  class A=std::allocator<X> >
64  class AMG : public Preconditioner<X,X>
65  {
66  template<class M1, class X1, class S1, class P1, class K1, class A1>
67  friend class KAMG;
68 
69  friend class KAmgTwoGrid<AMG>;
70 
71  public:
73  typedef M Operator;
80  typedef PI ParallelInformation;
85 
87  typedef X Domain;
89  typedef X Range;
97  typedef S Smoother;
98 
101 
111  AMG(OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
112  const SmootherArgs& smootherArgs, const Parameters& parms);
113 
125  template<class C>
126  AMG(const Operator& fineOperator, const C& criterion,
127  const SmootherArgs& smootherArgs=SmootherArgs(),
129 
180  AMG(std::shared_ptr<const Operator> fineOperator, const ParameterTree& configuration, const ParallelInformation& pinfo=ParallelInformation());
181 
185  AMG(const AMG& amg);
186 
188  void pre(Domain& x, Range& b);
189 
191  void apply(Domain& v, const Range& d);
192 
195  {
196  return category_;
197  }
198 
200  void post(Domain& x);
201 
206  template<class A1>
207  void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
208 
209  std::size_t levels();
210 
211  std::size_t maxlevels();
212 
222  {
223  matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
224  }
225 
231 
232  private:
233  /*
234  * @brief Helper function to create hierarchies with parameter tree.
235  *
236  * Will create the coarsen criterion with the norm and create the
237  * Hierarchies
238  * \tparam Norm Type of the norm to use.
239  */
240  template<class Norm>
241  void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
242  const PI& pinfo, const Norm&,
243  const ParameterTree& configuration,
244  std::true_type compiles = std::true_type());
245  template<class Norm>
246  void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
247  const PI& pinfo, const Norm&,
248  const ParameterTree& configuration,
249  std::false_type);
254  template<class C>
255  void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
256  const PI& pinfo, const ParameterTree& configuration);
263  template<class C>
264  void createHierarchies(C& criterion,
265  const std::shared_ptr<const Operator>& matrixptr,
266  const PI& pinfo);
273  struct LevelContext
274  {
291  typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
295  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
311  std::size_t level;
312  };
313 
314 
319  void mgc(LevelContext& levelContext);
320 
321  void additiveMgc();
322 
329  void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
330 
335  bool moveToCoarseLevel(LevelContext& levelContext);
336 
341  void initIteratorsWithFineLevel(LevelContext& levelContext);
342 
344  std::shared_ptr<OperatorHierarchy> matrices_;
346  SmootherArgs smootherArgs_;
348  std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
350  std::shared_ptr<CoarseSolver> solver_;
352  std::shared_ptr<Hierarchy<Range,A>> rhs_;
354  std::shared_ptr<Hierarchy<Domain,A>> lhs_;
356  std::shared_ptr<Hierarchy<Domain,A>> update_;
360  std::shared_ptr<ScalarProduct> scalarProduct_;
362  std::size_t gamma_;
364  std::size_t preSteps_;
366  std::size_t postSteps_;
367  bool buildHierarchy_;
368  bool additive;
369  bool coarsesolverconverged;
370  std::shared_ptr<Smoother> coarseSmoother_;
372  SolverCategory::Category category_;
374  std::size_t verbosity_;
375 
376  struct ToLower
377  {
378  std::string operator()(const std::string& str)
379  {
380  std::stringstream retval;
381  std::ostream_iterator<char> out(retval);
382  std::transform(str.begin(), str.end(), out,
383  [](char c){
384  return std::tolower(c, std::locale::classic());
385  });
386  return retval.str();
387  }
388  };
389  };
390 
391  template<class M, class X, class S, class PI, class A>
392  inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
393  : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
394  smoothers_(amg.smoothers_), solver_(amg.solver_),
395  rhs_(), lhs_(), update_(),
396  scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
397  preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
398  buildHierarchy_(amg.buildHierarchy_),
399  additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
400  coarseSmoother_(amg.coarseSmoother_),
401  category_(amg.category_),
402  verbosity_(amg.verbosity_)
403  {}
404 
405  template<class M, class X, class S, class PI, class A>
407  const SmootherArgs& smootherArgs,
408  const Parameters& parms)
409  : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
410  smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
411  rhs_(), lhs_(), update_(), scalarProduct_(0),
412  gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
413  postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
414  additive(parms.getAdditive()), coarsesolverconverged(true),
415  coarseSmoother_(),
416 // #warning should category be retrieved from matrices?
417  category_(SolverCategory::category(*smoothers_->coarsest())),
418  verbosity_(parms.debugLevel())
419  {
420  assert(matrices_->isBuilt());
421 
422  // build the necessary smoother hierarchies
423  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
424  }
425 
426  template<class M, class X, class S, class PI, class A>
427  template<class C>
429  const C& criterion,
430  const SmootherArgs& smootherArgs,
431  const PI& pinfo)
432  : smootherArgs_(smootherArgs),
433  smoothers_(new Hierarchy<Smoother,A>), solver_(),
434  rhs_(), lhs_(), update_(), scalarProduct_(),
435  gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
436  postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
437  additive(criterion.getAdditive()), coarsesolverconverged(true),
438  coarseSmoother_(),
439  category_(SolverCategory::category(pinfo)),
440  verbosity_(criterion.debugLevel())
441  {
443  DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
444  // TODO: reestablish compile time checks.
445  //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
446  // "Matrix and Solver must match in terms of category!");
447  auto matrixptr = stackobject_to_shared_ptr(matrix);
448  createHierarchies(criterion, matrixptr, pinfo);
449  }
450 
451  template<class M, class X, class S, class PI, class A>
452  AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
453  const ParameterTree& configuration,
454  const ParallelInformation& pinfo) :
455  smoothers_(new Hierarchy<Smoother,A>),
456  solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
457  coarsesolverconverged(true), coarseSmoother_(),
458  category_(SolverCategory::category(pinfo))
459  {
460 
461  if (configuration.hasKey ("smootherIterations"))
462  smootherArgs_.iterations = configuration.get<int>("smootherIterations");
463 
464  if (configuration.hasKey ("smootherRelaxation"))
465  smootherArgs_.relaxationFactor = configuration.get<typename SmootherArgs::RelaxationFactor>("smootherRelaxation");
466 
467  auto normName = ToLower()(configuration.get("strengthMeasure", "diagonal"));
468  auto index = configuration.get<int>("diagonalRowIndex", 0);
469 
470  if ( normName == "diagonal")
471  {
472  using field_type = typename M::field_type;
473  using real_type = typename FieldTraits<field_type>::real_type;
474  std::is_convertible<field_type, real_type> compiles;
475 
476  switch (index)
477  {
478  case 0:
479  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<0>(), configuration, compiles);
480  break;
481  case 1:
482  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<1>(), configuration, compiles);
483  break;
484  case 2:
485  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<2>(), configuration, compiles);
486  break;
487  case 3:
488  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<3>(), configuration, compiles);
489  break;
490  case 4:
491  createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<4>(), configuration, compiles);
492  break;
493  default:
494  DUNE_THROW(InvalidStateException, "Currently strengthIndex>4 is not supported.");
495  }
496  }
497  else if (normName == "rowsum")
498  createCriterionAndHierarchies(matrixptr, pinfo, RowSum(), configuration);
499  else if (normName == "frobenius")
500  createCriterionAndHierarchies(matrixptr, pinfo, FrobeniusNorm(), configuration);
501  else if (normName == "one")
502  createCriterionAndHierarchies(matrixptr, pinfo, AlwaysOneNorm(), configuration);
503  else
504  DUNE_THROW(Dune::NotImplemented, "Wrong config file: strengthMeasure "<<normName<<" is not supported by AMG");
505  }
506 
507  template<class M, class X, class S, class PI, class A>
508  template<class Norm>
509  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)
510  {
511  DUNE_THROW(InvalidStateException, "Strength of connection measure does not support this type ("
512  << className<typename M::field_type>() << ") as it is lacking a conversion to"
513  << className<typename FieldTraits<typename M::field_type>::real_type>() << ".");
514  }
515 
516  template<class M, class X, class S, class PI, class A>
517  template<class Norm>
518  void AMG<M,X,S,PI,A>::createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const Norm&, const ParameterTree& configuration, std::true_type)
519  {
520  if (configuration.get<bool>("criterionSymmetric", true))
521  {
522  using Criterion = Dune::Amg::CoarsenCriterion<
524  Criterion criterion;
525  createHierarchies(criterion, matrixptr, pinfo, configuration);
526  }
527  else
528  {
529  using Criterion = Dune::Amg::CoarsenCriterion<
531  Criterion criterion;
532  createHierarchies(criterion, matrixptr, pinfo, configuration);
533  }
534  }
535 
536  template<class M, class X, class S, class PI, class A>
537  template<class C>
538  void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const ParameterTree& configuration)
539  {
540  if (configuration.hasKey ("maxLevel"))
541  criterion.setMaxLevel(configuration.get<int>("maxLevel"));
542 
543  if (configuration.hasKey ("minCoarseningRate"))
544  criterion.setMinCoarsenRate(configuration.get<int>("minCoarseningRate"));
545 
546  if (configuration.hasKey ("coarsenTarget"))
547  criterion.setCoarsenTarget (configuration.get<int>("coarsenTarget"));
548 
549  if (configuration.hasKey ("accumulationMode"))
550  {
551  std::string mode = ToLower()(configuration.get<std::string>("accumulationMode"));
552  if ( mode == "none")
553  criterion.setAccumulate(AccumulationMode::noAccu);
554  else if ( mode == "atonce" )
555  criterion.setAccumulate(AccumulationMode::atOnceAccu);
556  else if ( mode == "successive")
557  criterion.setCoarsenTarget (AccumulationMode::successiveAccu);
558  else
559  DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
560  << mode <<".");
561  }
562 
563  if (configuration.hasKey ("prolongationDampingFactor"))
564  criterion.setProlongationDampingFactor (configuration.get<double>("prolongationDampingFactor"));
565 
566  if (configuration.hasKey("defaultAggregationSizeMode"))
567  {
568  auto mode = ToLower()(configuration.get<std::string>("defaultAggregationSizeMode"));
569  auto dim = configuration.get<std::size_t>("defaultAggregationDimension");
570  std::size_t maxDistance = 2;
571  if (configuration.hasKey("MaxAggregateDistance"))
572  maxDistance = configuration.get<std::size_t>("maxAggregateDistance");
573  if (mode == "isotropic")
574  criterion.setDefaultValuesIsotropic(dim, maxDistance);
575  else if(mode == "anisotropic")
576  criterion.setDefaultValuesAnisotropic(dim, maxDistance);
577  else
578  DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
579  << mode <<".");
580  }
581 
582  if (configuration.hasKey("maxAggregateDistance"))
583  criterion.setMaxDistance(configuration.get<std::size_t>("maxAggregateDistance"));
584 
585  if (configuration.hasKey("minAggregateSize"))
586  criterion.setMinAggregateSize(configuration.get<std::size_t>("minAggregateSize"));
587 
588  if (configuration.hasKey("maxAggregateSize"))
589  criterion.setMaxAggregateSize(configuration.get<std::size_t>("maxAggregateSize"));
590 
591  if (configuration.hasKey("maxAggregateConnectivity"))
592  criterion.setMaxConnectivity(configuration.get<std::size_t>("maxAggregateConnectivity"));
593 
594  if (configuration.hasKey ("alpha"))
595  criterion.setAlpha (configuration.get<double> ("alpha"));
596 
597  if (configuration.hasKey ("beta"))
598  criterion.setBeta (configuration.get<double> ("beta"));
599 
600  if (configuration.hasKey ("gamma"))
601  criterion.setGamma (configuration.get<std::size_t> ("gamma"));
602  gamma_ = criterion.getGamma();
603 
604  if (configuration.hasKey ("additive"))
605  criterion.setAdditive (configuration.get<bool>("additive"));
606  additive = criterion.getAdditive();
607 
608  if (configuration.hasKey ("preSteps"))
609  criterion.setNoPreSmoothSteps (configuration.get<std::size_t> ("preSteps"));
610  preSteps_ = criterion.getNoPreSmoothSteps ();
611 
612  if (configuration.hasKey ("postSteps"))
613  criterion.setNoPostSmoothSteps (configuration.get<std::size_t> ("postSteps"));
614  postSteps_ = criterion.getNoPostSmoothSteps ();
615 
616  verbosity_ = configuration.get("verbosity", 0);
617  criterion.setDebugLevel (verbosity_);
618 
619  createHierarchies(criterion, matrixptr, pinfo);
620  }
621 
622  template <class Matrix,
623  class Vector>
625  {
628 
629  static constexpr SolverType solver =
630 #if DISABLE_AMG_DIRECTSOLVER
631  none;
632 #elif HAVE_SUITESPARSE_UMFPACK
634 #elif HAVE_SUPERLU
635  superlu ;
636 #else
637  none;
638 #endif
639 
640  template <class M, SolverType>
641  struct Solver
642  {
644  static type* create(const M& mat, bool verbose, bool reusevector )
645  {
646  DUNE_THROW(NotImplemented,"DirectSolver not selected");
647  return nullptr;
648  }
649  static std::string name () { return "None"; }
650  };
651 #if HAVE_SUITESPARSE_UMFPACK
652  template <class M>
653  struct Solver< M, umfpack >
654  {
655  typedef UMFPack< M > type;
656  static type* create(const M& mat, bool verbose, bool reusevector )
657  {
658  return new type(mat, verbose, reusevector );
659  }
660  static std::string name () { return "UMFPack"; }
661  };
662 #endif
663 #if HAVE_SUPERLU
664  template <class M>
665  struct Solver< M, superlu >
666  {
668  static type* create(const M& mat, bool verbose, bool reusevector )
669  {
670  return new type(mat, verbose, reusevector );
671  }
672  static std::string name () { return "SuperLU"; }
673  };
674 #endif
675 
676  // define direct solver type to be used
679  static constexpr bool isDirectSolver = solver != none;
680  static std::string name() { return SelectedSolver :: name (); }
681  static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
682  {
683  return SelectedSolver :: create( mat, verbose, reusevector );
684  }
685  };
686 
687  template<class M, class X, class S, class PI, class A>
688  template<class C>
689  void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
690  const std::shared_ptr<const Operator>& matrixptr,
691  const PI& pinfo)
692  {
693  Timer watch;
694  matrices_ = std::make_shared<OperatorHierarchy>(
695  std::const_pointer_cast<Operator>(matrixptr),
696  stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
697 
698  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
699 
700  // build the necessary smoother hierarchies
701  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
702 
703  // test whether we should solve on the coarse level. That is the case if we
704  // have that level and if there was a redistribution on this level then our
705  // communicator has to be valid (size()>0) as the smoother might try to communicate
706  // in the constructor.
707  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
708  && ( ! matrices_->redistributeInformation().back().isSetup() ||
709  matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
710  {
711  // We have the carsest level. Create the coarse Solver
712  SmootherArgs sargs(smootherArgs_);
713  sargs.iterations = 1;
714 
716  cargs.setArgs(sargs);
717  if(matrices_->redistributeInformation().back().isSetup()) {
718  // Solve on the redistributed partitioning
719  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
720  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
721  }else{
722  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
723  cargs.setComm(*matrices_->parallelInformation().coarsest());
724  }
725 
726  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
727  scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
728 
729  typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
730 
731  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
732  if( SolverSelector::isDirectSolver &&
733  (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
734  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
735  || (matrices_->parallelInformation().coarsest().isRedistributed()
736  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
737  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
738  )
739  { // redistribute and 1 proc
740  if(matrices_->parallelInformation().coarsest().isRedistributed())
741  {
742  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
743  {
744  // We are still participating on this level
745  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
746  }
747  else
748  solver_.reset();
749  }
750  else
751  {
752  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
753  }
754  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
755  std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
756  }
757  else
758  {
759  if(matrices_->parallelInformation().coarsest().isRedistributed())
760  {
761  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
762  // We are still participating on this level
763 
764  // we have to allocate these types using the rebound allocator
765  // in order to ensure that we fulfill the alignment requirements
766  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
767  *scalarProduct_,
768  *coarseSmoother_, 1E-2, 1000, 0));
769  else
770  solver_.reset();
771  }else
772  {
773  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
774  *scalarProduct_,
775  *coarseSmoother_, 1E-2, 1000, 0));
776  // // we have to allocate these types using the rebound allocator
777  // // in order to ensure that we fulfill the alignment requirements
778  // using Alloc = typename std::allocator_traits<A>::template rebind_alloc<BiCGSTABSolver<X>>;
779  // Alloc alloc;
780  // auto p = alloc.allocate(1);
781  // std::allocator_traits<Alloc>::construct(alloc, p,
782  // const_cast<M&>(*matrices_->matrices().coarsest()),
783  // *scalarProduct_,
784  // *coarseSmoother_, 1E-2, 1000, 0);
785  // solver_.reset(p,[](BiCGSTABSolver<X>* p){
786  // Alloc alloc;
787  // std::allocator_traits<Alloc>::destroy(alloc, p);
788  // alloc.deallocate(p,1);
789  // });
790  }
791  }
792  }
793 
794  if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
795  std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
796  <<"(including coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
797  }
798 
799 
800  template<class M, class X, class S, class PI, class A>
802  {
803  // Detect Matrix rows where all offdiagonal entries are
804  // zero and set x such that A_dd*x_d=b_d
805  // Thus users can be more careless when setting up their linear
806  // systems.
807  typedef typename M::matrix_type Matrix;
808  typedef typename Matrix::ConstRowIterator RowIter;
809  typedef typename Matrix::ConstColIterator ColIter;
810  typedef typename Matrix::block_type Block;
811  Block zero;
812  zero=typename Matrix::field_type();
813 
814  const Matrix& mat=matrices_->matrices().finest()->getmat();
815  for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
816  bool isDirichlet = true;
817  bool hasDiagonal = false;
818  Block diagonal{};
819  for(ColIter col=row->begin(); col!=row->end(); ++col) {
820  if(row.index()==col.index()) {
821  diagonal = *col;
822  hasDiagonal = true;
823  }else{
824  if(*col!=zero)
825  isDirichlet = false;
826  }
827  }
828  if(isDirichlet && hasDiagonal)
829  {
830  auto&& xEntry = Impl::asVector(x[row.index()]);
831  auto&& bEntry = Impl::asVector(b[row.index()]);
832  Impl::asMatrix(diagonal).solve(xEntry, bEntry);
833  }
834  }
835 
836  if(smoothers_->levels()>0)
837  smoothers_->finest()->pre(x,b);
838  else
839  // No smoother to make x consistent! Do it by hand
840  matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
841  rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
842  lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
843  update_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
844  matrices_->coarsenVector(*rhs_);
845  matrices_->coarsenVector(*lhs_);
846  matrices_->coarsenVector(*update_);
847 
848  // Preprocess all smoothers
849  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
850  typedef typename Hierarchy<Range,A>::Iterator RIterator;
851  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
852  Iterator coarsest = smoothers_->coarsest();
853  Iterator smoother = smoothers_->finest();
854  RIterator rhs = rhs_->finest();
855  DIterator lhs = lhs_->finest();
856  if(smoothers_->levels()>1) {
857 
858  assert(lhs_->levels()==rhs_->levels());
859  assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
860  assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
861 
862  if(smoother!=coarsest)
863  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
864  smoother->pre(*lhs,*rhs);
865  smoother->pre(*lhs,*rhs);
866  }
867 
868 
869  // The preconditioner might change x and b. So we have to
870  // copy the changes to the original vectors.
871  x = *lhs_->finest();
872  b = *rhs_->finest();
873 
874  }
875  template<class M, class X, class S, class PI, class A>
877  {
878  return matrices_->levels();
879  }
880  template<class M, class X, class S, class PI, class A>
882  {
883  return matrices_->maxlevels();
884  }
885 
887  template<class M, class X, class S, class PI, class A>
889  {
890  LevelContext levelContext;
891 
892  if(additive) {
893  *(rhs_->finest())=d;
894  additiveMgc();
895  v=*lhs_->finest();
896  }else{
897  // Init all iterators for the current level
898  initIteratorsWithFineLevel(levelContext);
899 
900 
901  *levelContext.lhs = v;
902  *levelContext.rhs = d;
903  *levelContext.update=0;
904  levelContext.level=0;
905 
906  mgc(levelContext);
907 
908  if(postSteps_==0||matrices_->maxlevels()==1)
909  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
910 
911  v=*levelContext.update;
912  }
913 
914  }
915 
916  template<class M, class X, class S, class PI, class A>
917  void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
918  {
919  levelContext.smoother = smoothers_->finest();
920  levelContext.matrix = matrices_->matrices().finest();
921  levelContext.pinfo = matrices_->parallelInformation().finest();
922  levelContext.redist =
923  matrices_->redistributeInformation().begin();
924  levelContext.aggregates = matrices_->aggregatesMaps().begin();
925  levelContext.lhs = lhs_->finest();
926  levelContext.update = update_->finest();
927  levelContext.rhs = rhs_->finest();
928  }
929 
930  template<class M, class X, class S, class PI, class A>
931  bool AMG<M,X,S,PI,A>
932  ::moveToCoarseLevel(LevelContext& levelContext)
933  {
934 
935  bool processNextLevel=true;
936 
937  if(levelContext.redist->isSetup()) {
938  levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
939  levelContext.rhs.getRedistributed());
940  processNextLevel = levelContext.rhs.getRedistributed().size()>0;
941  if(processNextLevel) {
942  //restrict defect to coarse level right hand side.
943  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
944  ++levelContext.pinfo;
946  ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
947  static_cast<const Range&>(fineRhs.getRedistributed()),
948  *levelContext.pinfo);
949  }
950  }else{
951  //restrict defect to coarse level right hand side.
952  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
953  ++levelContext.pinfo;
955  ::restrictVector(*(*levelContext.aggregates),
956  *levelContext.rhs, static_cast<const Range&>(*fineRhs),
957  *levelContext.pinfo);
958  }
959 
960  if(processNextLevel) {
961  // prepare coarse system
962  ++levelContext.lhs;
963  ++levelContext.update;
964  ++levelContext.matrix;
965  ++levelContext.level;
966  ++levelContext.redist;
967 
968  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
969  // next level is not the globally coarsest one
970  ++levelContext.smoother;
971  ++levelContext.aggregates;
972  }
973  // prepare the update on the next level
974  *levelContext.update=0;
975  }
976  return processNextLevel;
977  }
978 
979  template<class M, class X, class S, class PI, class A>
980  void AMG<M,X,S,PI,A>
981  ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
982  {
983  if(processNextLevel) {
984  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
985  // previous level is not the globally coarsest one
986  --levelContext.smoother;
987  --levelContext.aggregates;
988  }
989  --levelContext.redist;
990  --levelContext.level;
991  //prolongate and add the correction (update is in coarse left hand side)
992  --levelContext.matrix;
993 
994  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
995  --levelContext.lhs;
996  --levelContext.pinfo;
997  }
998  if(levelContext.redist->isSetup()) {
999  // Need to redistribute during prolongateVector
1000  levelContext.lhs.getRedistributed()=0;
1002  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1003  levelContext.lhs.getRedistributed(),
1004  matrices_->getProlongationDampingFactor(),
1005  *levelContext.pinfo, *levelContext.redist);
1006  }else{
1007  *levelContext.lhs=0;
1009  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1010  matrices_->getProlongationDampingFactor(),
1011  *levelContext.pinfo);
1012  }
1013 
1014 
1015  if(processNextLevel) {
1016  --levelContext.update;
1017  --levelContext.rhs;
1018  }
1019 
1020  *levelContext.update += *levelContext.lhs;
1021  }
1022 
1023  template<class M, class X, class S, class PI, class A>
1025  {
1027  }
1028 
1029  template<class M, class X, class S, class PI, class A>
1030  void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
1031  if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
1032  // Solve directly
1034  res.converged=true; // If we do not compute this flag will not get updated
1035  if(levelContext.redist->isSetup()) {
1036  levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
1037  if(levelContext.rhs.getRedistributed().size()>0) {
1038  // We are still participating in the computation
1039  levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
1040  levelContext.rhs.getRedistributed());
1041  solver_->apply(levelContext.update.getRedistributed(),
1042  levelContext.rhs.getRedistributed(), res);
1043  }
1044  levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
1045  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
1046  }else{
1047  levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
1048  solver_->apply(*levelContext.update, *levelContext.rhs, res);
1049  }
1050 
1051  if (!res.converged)
1052  coarsesolverconverged = false;
1053  }else{
1054  // presmoothing
1055  presmooth(levelContext, preSteps_);
1056 
1057 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1058  bool processNextLevel = moveToCoarseLevel(levelContext);
1059 
1060  if(processNextLevel) {
1061  // next level
1062  for(std::size_t i=0; i<gamma_; i++){
1063  mgc(levelContext);
1064  if (levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels())
1065  break;
1066  if(i+1 < gamma_){
1067  levelContext.matrix->applyscaleadd(-1., *levelContext.lhs, *levelContext.rhs);
1068  }
1069  }
1070  }
1071 
1072  moveToFineLevel(levelContext, processNextLevel);
1073 #else
1074  *lhs=0;
1075 #endif
1076 
1077  if(levelContext.matrix == matrices_->matrices().finest()) {
1078  coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
1079  if(!coarsesolverconverged)
1080  DUNE_THROW(MathError, "Coarse solver did not converge");
1081  }
1082  // postsmoothing
1083  postsmooth(levelContext, postSteps_);
1084 
1085  }
1086  }
1087 
1088  template<class M, class X, class S, class PI, class A>
1089  void AMG<M,X,S,PI,A>::additiveMgc(){
1090 
1091  // restrict residual to all levels
1092  typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
1093  typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
1094  typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
1095  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
1096 
1097  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
1098  ++pinfo;
1100  ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
1101  }
1102 
1103  // pinfo is invalid, set to coarsest level
1104  //pinfo = matrices_->parallelInformation().coarsest
1105  // calculate correction for all levels
1106  lhs = lhs_->finest();
1107  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
1108 
1109  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
1110  // presmoothing
1111  *lhs=0;
1112  smoother->apply(*lhs, *rhs);
1113  }
1114 
1115  // Coarse level solve
1116 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1117  InverseOperatorResult res;
1118  pinfo->copyOwnerToAll(*rhs, *rhs);
1119  solver_->apply(*lhs, *rhs, res);
1120 
1121  if(!res.converged)
1122  DUNE_THROW(MathError, "Coarse solver did not converge");
1123 #else
1124  *lhs=0;
1125 #endif
1126  // Prologate and add up corrections from all levels
1127  --pinfo;
1128  --aggregates;
1129 
1130  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
1132  ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
1133  }
1134  }
1135 
1136 
1138  template<class M, class X, class S, class PI, class A>
1139  void AMG<M,X,S,PI,A>::post([[maybe_unused]] Domain& x)
1140  {
1141  // Postprocess all smoothers
1142  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
1143  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
1144  Iterator coarsest = smoothers_->coarsest();
1145  Iterator smoother = smoothers_->finest();
1146  DIterator lhs = lhs_->finest();
1147  if(smoothers_->levels()>0) {
1148  if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
1149  smoother->post(*lhs);
1150  if(smoother!=coarsest)
1151  for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
1152  smoother->post(*lhs);
1153  smoother->post(*lhs);
1154  }
1155  lhs_ = nullptr;
1156  update_ = nullptr;
1157  rhs_ = nullptr;
1158  }
1159 
1160  template<class M, class X, class S, class PI, class A>
1161  template<class A1>
1162  void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
1163  {
1164  matrices_->getCoarsestAggregatesOnFinest(cont);
1165  }
1166 
1167  } // end namespace Amg
1168 
1169  struct AMGCreator{
1170  template<class> struct isValidBlockType : std::false_type{};
1171  template<class T, int n, int m> struct isValidBlockType<FieldMatrix<T,n,m>> : std::true_type{};
1172 
1173  template<class OP>
1174  std::shared_ptr<Dune::Preconditioner<typename OP::element_type::domain_type, typename OP::element_type::range_type> >
1175  makeAMG(const OP& op, const std::string& smoother, const Dune::ParameterTree& config) const
1176  {
1177  DUNE_THROW(Dune::Exception, "Operator type not supported by AMG");
1178  }
1179 
1180  template<class M, class X, class Y>
1181  std::shared_ptr<Dune::Preconditioner<X,Y> >
1182  makeAMG(const std::shared_ptr<MatrixAdapter<M,X,Y>>& op, const std::string& smoother,
1183  const Dune::ParameterTree& config) const
1184  {
1185  using OP = MatrixAdapter<M,X,Y>;
1186 
1187  if(smoother == "ssor")
1188  return std::make_shared<Amg::AMG<OP, X, SeqSSOR<M,X,Y>>>(op, config);
1189  if(smoother == "sor")
1190  return std::make_shared<Amg::AMG<OP, X, SeqSOR<M,X,Y>>>(op, config);
1191  if(smoother == "jac")
1192  return std::make_shared<Amg::AMG<OP, X, SeqJac<M,X,Y>>>(op, config);
1193  if(smoother == "gs")
1194  return std::make_shared<Amg::AMG<OP, X, SeqGS<M,X,Y>>>(op, config);
1195  if(smoother == "ilu")
1196  return std::make_shared<Amg::AMG<OP, X, SeqILU<M,X,Y>>>(op, config);
1197  else
1198  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1199  }
1200 
1201  template<class M, class X, class Y, class C>
1202  std::shared_ptr<Dune::Preconditioner<X,Y> >
1203  makeAMG(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1204  const Dune::ParameterTree& config) const
1205  {
1207 
1208  auto cop = std::static_pointer_cast<const OP>(op);
1209 
1210  if(smoother == "ssor")
1211  return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1212  if(smoother == "sor")
1213  return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1214  if(smoother == "jac")
1215  return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqJac<M,X,Y>>,C>>(cop, config, op->getCommunication());
1216  if(smoother == "gs")
1217  return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqGS<M,X,Y>>,C>>(cop, config, op->getCommunication());
1218  if(smoother == "ilu")
1219  return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqILU<M,X,Y>>,C>>(cop, config, op->getCommunication());
1220  else
1221  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1222  }
1223 
1224  template<class M, class X, class Y, class C>
1225  std::shared_ptr<Dune::Preconditioner<X,Y> >
1226  makeAMG(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1227  const Dune::ParameterTree& config) const
1228  {
1230 
1231  if(smoother == "ssor")
1232  return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1233  if(smoother == "sor")
1234  return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1235  if(smoother == "jac")
1236  return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqJac<M,X,Y>>,C>>(op, config, op->getCommunication());
1237  if(smoother == "gs")
1238  return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqGS<M,X,Y>>,C>>(op, config, op->getCommunication());
1239  if(smoother == "ilu")
1240  return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqILU<M,X,Y>>,C>>(op, config, op->getCommunication());
1241  else
1242  DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1243  }
1244 
1245  template<typename TL, typename OP>
1246  std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1247  typename Dune::TypeListElement<2, TL>::type>>
1248  operator() (TL tl, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
1249  std::enable_if_t<isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
1250  {
1251  using field_type = typename OP::matrix_type::field_type;
1252  using real_type = typename FieldTraits<field_type>::real_type;
1253  if (!std::is_convertible<field_type, real_type>())
1254  DUNE_THROW(UnsupportedType, "AMG needs field_type(" <<
1255  className<field_type>() <<
1256  ") to be convertible to its real_type (" <<
1257  className<real_type>() <<
1258  ").");
1259  using D = typename Dune::TypeListElement<1, decltype(tl)>::type;
1260  using R = typename Dune::TypeListElement<2, decltype(tl)>::type;
1261  std::shared_ptr<Preconditioner<D,R>> amg;
1262  std::string smoother = config.get("smoother", "ssor");
1263  return makeAMG(op, smoother, config);
1264  }
1265 
1266  template<typename TL, typename OP>
1267  std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1268  typename Dune::TypeListElement<2, TL>::type>>
1269  operator() (TL /*tl*/, const std::shared_ptr<OP>& /*mat*/, const Dune::ParameterTree& /*config*/,
1270  std::enable_if_t<!isValidBlockType<typename OP::matrix_type::block_type>::value,int> = 0) const
1271  {
1272  DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
1273  }
1274  };
1275 
1277 } // end namespace Dune
1278 
1279 #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:351
Matrix & mat
Definition: matrixmatrix.hh:347
AMG(const AMG &amg)
Copy constructor.
Definition: amg.hh:392
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:801
static std::string name()
Definition: amg.hh:680
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:303
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:307
std::shared_ptr< Dune::Preconditioner< typename OP::element_type::domain_type, typename OP::element_type::range_type > > makeAMG(const OP &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1175
static std::string name()
Definition: amg.hh:672
Matrix ::field_type field_type
Definition: amg.hh:626
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:1024
X Domain
The domain type.
Definition: amg.hh:87
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:681
AMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:406
AMG(std::shared_ptr< const Operator > fineOperator, const ParameterTree &configuration, const ParallelInformation &pinfo=ParallelInformation())
Constructor an AMG via ParameterTree.
Definition: amg.hh:452
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:287
SolverType
Definition: amg.hh:627
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:295
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:100
Solver< Matrix, solver > SelectedSolver
Definition: amg.hh:677
std::string operator()(const std::string &str)
Definition: amg.hh:378
S Smoother
The type of the smoother.
Definition: amg.hh:97
static std::string name()
Definition: amg.hh:649
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:279
M Operator
The matrix operator type.
Definition: amg.hh:73
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:283
SelectedSolver ::type DirectSolver
Definition: amg.hh:678
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:291
X Range
The range type.
Definition: amg.hh:89
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:406
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:1162
std::size_t levels()
Definition: amg.hh:876
InverseOperator< Vector, Vector > type
Definition: amg.hh:643
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< MatrixAdapter< M, X, Y >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1182
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:299
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< NonoverlappingSchwarzOperator< M, X, Y, C >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1226
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
static constexpr SolverType solver
Definition: amg.hh:629
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:52
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< OverlappingSchwarzOperator< M, X, Y, C >> &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1203
static constexpr bool isDirectSolver
Definition: amg.hh:679
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:221
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:383
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:84
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:91
void post(Domain &x)
Clean up.
Definition: amg.hh:1139
std::size_t maxlevels()
Definition: amg.hh:881
std::shared_ptr< Dune::Preconditioner< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL tl, const std::shared_ptr< OP > &op, const Dune::ParameterTree &config, std::enable_if_t< isValidBlockType< typename OP::matrix_type::block_type >::value, int >=0) const
Definition: amg.hh:1248
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:428
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:644
std::size_t level
The level index.
Definition: amg.hh:311
AMG(const Operator &fineOperator, const C &criterion, const SmootherArgs &smootherArgs=SmootherArgs(), const ParallelInformation &pinfo=ParallelInformation())
Construct an AMG with an inexact coarse solver based on the smoother.
Definition: amg.hh:428
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:888
Smoother SmootherType
Definition: amg.hh:275
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:82
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amg.hh:194
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:668
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:80
@ none
Definition: amg.hh:627
@ umfpack
Definition: amg.hh:627
@ superlu
Definition: amg.hh:627
@ atOnceAccu
Accumulate data to one process at once.
Definition: parameters.hh:244
@ noAccu
No data accumulution.
Definition: parameters.hh:238
@ successiveAccu
Successively accumulate to fewer processes.
Definition: parameters.hh:248
Definition: allocator.hh:11
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
ConstIterator class for sequential access.
Definition: matrix.hh:404
A generic dynamic dense matrix.
Definition: matrix.hh:561
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:589
T block_type
Export the type representing the components.
Definition: matrix.hh:568
Definition: matrixutils.hh:27
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:61
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:137
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:463
Definition: aggregates.hh:480
Definition: aggregates.hh:496
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:519
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:539
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:140
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:33
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:65
Definition: amg.hh:625
Definition: amg.hh:1169
Definition: amg.hh:1170
An overlapping Schwarz operator.
Definition: schwarz.hh:75
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:216
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:219
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:61
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:283
All parameters for AMG.
Definition: parameters.hh:393
Definition: pinfo.hh:28
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:66
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:32
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioner.hh:39
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:52
Statistics about the application of an inverse operator.
Definition: solver.hh:48
bool converged
True if convergence criterion has been met.
Definition: solver.hh:73
Categories for the solvers.
Definition: solvercategory.hh:22
Category
Definition: solvercategory.hh:23
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:34
Definition: solvercategory.hh:54
Definition: solverregistry.hh:77
Definition: solvertype.hh:16
SuperLu Solver.
Definition: superlu.hh:271
Definition: umfpack.hh:49
The UMFPack direct sparse solver.
Definition: umfpack.hh:215