dune-istl  3.0-git
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 <dune/common/exceptions.hh>
11 #include <dune/istl/solvers.hh>
13 #include <dune/istl/superlu.hh>
14 #include <dune/istl/umfpack.hh>
15 #include <dune/istl/solvertype.hh>
16 #include <dune/common/typetraits.hh>
17 #include <dune/common/exceptions.hh>
18 
19 namespace Dune
20 {
21  namespace Amg
22  {
40  template<class M, class X, class S, class P, class K, class A>
41  class KAMG;
42 
43  template<class T>
44  class KAmgTwoGrid;
45 
53  template<class M, class X, class S, class PI=SequentialInformation,
54  class A=std::allocator<X> >
55  class AMG : public Preconditioner<X,X>
56  {
57  template<class M1, class X1, class S1, class P1, class K1, class A1>
58  friend class KAMG;
59 
60  friend class KAmgTwoGrid<AMG>;
61 
62  public:
64  typedef M Operator;
71  typedef PI ParallelInformation;
76 
78  typedef X Domain;
80  typedef X Range;
88  typedef S Smoother;
89 
92 
93  enum {
95  category = S::category
96  };
97 
107  AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
108  const SmootherArgs& smootherArgs, const Parameters& parms);
109 
121  template<class C>
122  AMG(const Operator& fineOperator, const C& criterion,
123  const SmootherArgs& smootherArgs=SmootherArgs(),
124  const ParallelInformation& pinfo=ParallelInformation());
125 
129  AMG(const AMG& amg);
130 
131  ~AMG();
132 
134  void pre(Domain& x, Range& b);
135 
137  void apply(Domain& v, const Range& d);
138 
140  void post(Domain& x);
141 
146  template<class A1>
147  void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
148 
149  std::size_t levels();
150 
151  std::size_t maxlevels();
152 
162  {
163  matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
164  }
165 
170  bool usesDirectCoarseLevelSolver() const;
171 
172  private:
179  template<class C>
180  void createHierarchies(C& criterion, Operator& matrix,
181  const PI& pinfo);
188  struct LevelContext
189  {
190  typedef Smoother SmootherType;
206  typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
210  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
226  std::size_t level;
227  };
228 
229 
234  void mgc(LevelContext& levelContext);
235 
236  void additiveMgc();
237 
244  void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
245 
250  bool moveToCoarseLevel(LevelContext& levelContext);
251 
256  void initIteratorsWithFineLevel(LevelContext& levelContext);
257 
259  std::shared_ptr<OperatorHierarchy> matrices_;
261  SmootherArgs smootherArgs_;
263  std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
265  std::shared_ptr<CoarseSolver> solver_;
267  Hierarchy<Range,A>* rhs_;
269  Hierarchy<Domain,A>* lhs_;
271  Hierarchy<Domain,A>* update_;
275  typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
276  typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
278  ScalarProductPointer scalarProduct_;
280  std::size_t gamma_;
282  std::size_t preSteps_;
284  std::size_t postSteps_;
285  bool buildHierarchy_;
286  bool additive;
287  bool coarsesolverconverged;
288  std::shared_ptr<Smoother> coarseSmoother_;
290  std::size_t verbosity_;
291  };
292 
293  template<class M, class X, class S, class PI, class A>
294  inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
295  : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
296  smoothers_(amg.smoothers_), solver_(amg.solver_),
297  rhs_(), lhs_(), update_(),
298  scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
299  preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
300  buildHierarchy_(amg.buildHierarchy_),
301  additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
302  coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
303  {
304  if(amg.rhs_)
305  rhs_=new Hierarchy<Range,A>(*amg.rhs_);
306  if(amg.lhs_)
307  lhs_=new Hierarchy<Domain,A>(*amg.lhs_);
308  if(amg.update_)
309  update_=new Hierarchy<Domain,A>(*amg.update_);
310  }
311 
312  template<class M, class X, class S, class PI, class A>
313  AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
314  const SmootherArgs& smootherArgs,
315  const Parameters& parms)
316  : matrices_(&matrices), smootherArgs_(smootherArgs),
317  smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
318  rhs_(), lhs_(), update_(), scalarProduct_(0),
319  gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
320  postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
321  additive(parms.getAdditive()), coarsesolverconverged(true),
322  coarseSmoother_(), verbosity_(parms.debugLevel())
323  {
324  assert(matrices_->isBuilt());
325 
326  // build the necessary smoother hierarchies
327  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
328  }
329 
330  template<class M, class X, class S, class PI, class A>
331  template<class C>
332  AMG<M,X,S,PI,A>::AMG(const Operator& matrix,
333  const C& criterion,
334  const SmootherArgs& smootherArgs,
335  const PI& pinfo)
336  : smootherArgs_(smootherArgs),
337  smoothers_(new Hierarchy<Smoother,A>), solver_(),
338  rhs_(), lhs_(), update_(), scalarProduct_(),
339  gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
340  postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
341  additive(criterion.getAdditive()), coarsesolverconverged(true),
342  coarseSmoother_(), verbosity_(criterion.debugLevel())
343  {
344  static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
345  "Matrix and Solver must match in terms of category!");
346  // TODO: reestablish compile time checks.
347  //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
348  // "Matrix and Solver must match in terms of category!");
349  createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
350  }
351 
352 
353  template<class M, class X, class S, class PI, class A>
355  {
356  if(buildHierarchy_) {
357  if(solver_)
358  solver_.reset();
359  if(coarseSmoother_)
360  coarseSmoother_.reset();
361  }
362  if(lhs_)
363  delete lhs_;
364  lhs_=nullptr;
365  if(update_)
366  delete update_;
367  update_=nullptr;
368  if(rhs_)
369  delete rhs_;
370  rhs_=nullptr;
371  }
372 
373  template<class M, class X, class S, class PI, class A>
374  template<class C>
375  void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
376  const PI& pinfo)
377  {
378  Timer watch;
379  matrices_.reset(new OperatorHierarchy(matrix, pinfo));
380 
381  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
382 
383  // build the necessary smoother hierarchies
384  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
385 
386  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
387  // We have the carsest level. Create the coarse Solver
388  SmootherArgs sargs(smootherArgs_);
389  sargs.iterations = 1;
390 
392  cargs.setArgs(sargs);
393  if(matrices_->redistributeInformation().back().isSetup()) {
394  // Solve on the redistributed partitioning
395  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
396  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
397  }else{
398  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
399  cargs.setComm(*matrices_->parallelInformation().coarsest());
400  }
401 
402  coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
403  scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
404 
405 #if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
406 #if HAVE_SUITESPARSE_UMFPACK
407 #define DIRECTSOLVER UMFPack
408 #else
409 #define DIRECTSOLVER SuperLU
410 #endif
411  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
412  if(std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
413  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
414  || (matrices_->parallelInformation().coarsest().isRedistributed()
415  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
416  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) { // redistribute and 1 proc
417  if(matrices_->parallelInformation().coarsest().isRedistributed())
418  {
419  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
420  // We are still participating on this level
421  solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
422  else
423  solver_.reset();
424  }else
425  solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(), false, false));
426  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
427  std::cout<< "Using a direct coarse solver (" << static_cast< DIRECTSOLVER<typename M::matrix_type>* >(solver_.get())->name() << ")" << std::endl;
428  }else
429 #undef DIRECTSOLVER
430 #endif // HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK
431  {
432  if(matrices_->parallelInformation().coarsest().isRedistributed())
433  {
434  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
435  // We are still participating on this level
436  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
437  *scalarProduct_,
438  *coarseSmoother_, 1E-2, 1000, 0));
439  else
440  solver_.reset();
441  }else
442  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
443  *scalarProduct_,
444  *coarseSmoother_, 1E-2, 1000, 0));
445  }
446  }
447 
448  if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
449  std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
450  <<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
451  }
452 
453 
454  template<class M, class X, class S, class PI, class A>
455  void AMG<M,X,S,PI,A>::pre(Domain& x, Range& b)
456  {
457  // Detect Matrix rows where all offdiagonal entries are
458  // zero and set x such that A_dd*x_d=b_d
459  // Thus users can be more careless when setting up their linear
460  // systems.
461  typedef typename M::matrix_type Matrix;
462  typedef typename Matrix::ConstRowIterator RowIter;
463  typedef typename Matrix::ConstColIterator ColIter;
464  typedef typename Matrix::block_type Block;
465  Block zero;
466  zero=typename Matrix::field_type();
467 
468  const Matrix& mat=matrices_->matrices().finest()->getmat();
469  for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
470  bool isDirichlet = true;
471  bool hasDiagonal = false;
472  Block diagonal;
473  for(ColIter col=row->begin(); col!=row->end(); ++col) {
474  if(row.index()==col.index()) {
475  diagonal = *col;
476  hasDiagonal = false;
477  }else{
478  if(*col!=zero)
479  isDirichlet = false;
480  }
481  }
482  if(isDirichlet && hasDiagonal)
483  diagonal.solve(x[row.index()], b[row.index()]);
484  }
485 
486  if(smoothers_->levels()>0)
487  smoothers_->finest()->pre(x,b);
488  else
489  // No smoother to make x consistent! Do it by hand
490  matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
491  Range* copy = new Range(b);
492  if(rhs_)
493  delete rhs_;
494  rhs_ = new Hierarchy<Range,A>(copy);
495  Domain* dcopy = new Domain(x);
496  if(lhs_)
497  delete lhs_;
498  lhs_ = new Hierarchy<Domain,A>(dcopy);
499  dcopy = new Domain(x);
500  if(update_)
501  delete update_;
502  update_ = new Hierarchy<Domain,A>(dcopy);
503  matrices_->coarsenVector(*rhs_);
504  matrices_->coarsenVector(*lhs_);
505  matrices_->coarsenVector(*update_);
506 
507  // Preprocess all smoothers
508  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
509  typedef typename Hierarchy<Range,A>::Iterator RIterator;
510  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
511  Iterator coarsest = smoothers_->coarsest();
512  Iterator smoother = smoothers_->finest();
513  RIterator rhs = rhs_->finest();
514  DIterator lhs = lhs_->finest();
515  if(smoothers_->levels()>0) {
516 
517  assert(lhs_->levels()==rhs_->levels());
518  assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
519  assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
520 
521  if(smoother!=coarsest)
522  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
523  smoother->pre(*lhs,*rhs);
524  smoother->pre(*lhs,*rhs);
525  }
526 
527 
528  // The preconditioner might change x and b. So we have to
529  // copy the changes to the original vectors.
530  x = *lhs_->finest();
531  b = *rhs_->finest();
532 
533  }
534  template<class M, class X, class S, class PI, class A>
536  {
537  return matrices_->levels();
538  }
539  template<class M, class X, class S, class PI, class A>
541  {
542  return matrices_->maxlevels();
543  }
544 
546  template<class M, class X, class S, class PI, class A>
547  void AMG<M,X,S,PI,A>::apply(Domain& v, const Range& d)
548  {
549  LevelContext levelContext;
550 
551  if(additive) {
552  *(rhs_->finest())=d;
553  additiveMgc();
554  v=*lhs_->finest();
555  }else{
556  // Init all iterators for the current level
557  initIteratorsWithFineLevel(levelContext);
558 
559 
560  *levelContext.lhs = v;
561  *levelContext.rhs = d;
562  *levelContext.update=0;
563  levelContext.level=0;
564 
565  mgc(levelContext);
566 
567  if(postSteps_==0||matrices_->maxlevels()==1)
568  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
569 
570  v=*levelContext.update;
571  }
572 
573  }
574 
575  template<class M, class X, class S, class PI, class A>
576  void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
577  {
578  levelContext.smoother = smoothers_->finest();
579  levelContext.matrix = matrices_->matrices().finest();
580  levelContext.pinfo = matrices_->parallelInformation().finest();
581  levelContext.redist =
582  matrices_->redistributeInformation().begin();
583  levelContext.aggregates = matrices_->aggregatesMaps().begin();
584  levelContext.lhs = lhs_->finest();
585  levelContext.update = update_->finest();
586  levelContext.rhs = rhs_->finest();
587  }
588 
589  template<class M, class X, class S, class PI, class A>
590  bool AMG<M,X,S,PI,A>
591  ::moveToCoarseLevel(LevelContext& levelContext)
592  {
593 
594  bool processNextLevel=true;
595 
596  if(levelContext.redist->isSetup()) {
597  levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
598  levelContext.rhs.getRedistributed());
599  processNextLevel = levelContext.rhs.getRedistributed().size()>0;
600  if(processNextLevel) {
601  //restrict defect to coarse level right hand side.
602  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
603  ++levelContext.pinfo;
605  ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
606  static_cast<const Range&>(fineRhs.getRedistributed()),
607  *levelContext.pinfo);
608  }
609  }else{
610  //restrict defect to coarse level right hand side.
611  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
612  ++levelContext.pinfo;
614  ::restrictVector(*(*levelContext.aggregates),
615  *levelContext.rhs, static_cast<const Range&>(*fineRhs),
616  *levelContext.pinfo);
617  }
618 
619  if(processNextLevel) {
620  // prepare coarse system
621  ++levelContext.lhs;
622  ++levelContext.update;
623  ++levelContext.matrix;
624  ++levelContext.level;
625  ++levelContext.redist;
626 
627  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
628  // next level is not the globally coarsest one
629  ++levelContext.smoother;
630  ++levelContext.aggregates;
631  }
632  // prepare the update on the next level
633  *levelContext.update=0;
634  }
635  return processNextLevel;
636  }
637 
638  template<class M, class X, class S, class PI, class A>
639  void AMG<M,X,S,PI,A>
640  ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
641  {
642  if(processNextLevel) {
643  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
644  // previous level is not the globally coarsest one
645  --levelContext.smoother;
646  --levelContext.aggregates;
647  }
648  --levelContext.redist;
649  --levelContext.level;
650  //prolongate and add the correction (update is in coarse left hand side)
651  --levelContext.matrix;
652 
653  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
654  --levelContext.lhs;
655  --levelContext.pinfo;
656  }
657  if(levelContext.redist->isSetup()) {
658  // Need to redistribute during prolongateVector
659  levelContext.lhs.getRedistributed()=0;
661  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
662  levelContext.lhs.getRedistributed(),
663  matrices_->getProlongationDampingFactor(),
664  *levelContext.pinfo, *levelContext.redist);
665  }else{
666  *levelContext.lhs=0;
668  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
669  matrices_->getProlongationDampingFactor(),
670  *levelContext.pinfo);
671  }
672 
673 
674  if(processNextLevel) {
675  --levelContext.update;
676  --levelContext.rhs;
677  }
678 
679  *levelContext.update += *levelContext.lhs;
680  }
681 
682  template<class M, class X, class S, class PI, class A>
684  {
686  }
687 
688  template<class M, class X, class S, class PI, class A>
689  void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
690  if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
691  // Solve directly
693  res.converged=true; // If we do not compute this flag will not get updated
694  if(levelContext.redist->isSetup()) {
695  levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
696  if(levelContext.rhs.getRedistributed().size()>0) {
697  // We are still participating in the computation
698  levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
699  levelContext.rhs.getRedistributed());
700  solver_->apply(levelContext.update.getRedistributed(),
701  levelContext.rhs.getRedistributed(), res);
702  }
703  levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
704  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
705  }else{
706  levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
707  solver_->apply(*levelContext.update, *levelContext.rhs, res);
708  }
709 
710  if (!res.converged)
711  coarsesolverconverged = false;
712  }else{
713  // presmoothing
714  presmooth(levelContext, preSteps_);
715 
716 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
717  bool processNextLevel = moveToCoarseLevel(levelContext);
718 
719  if(processNextLevel) {
720  // next level
721  for(std::size_t i=0; i<gamma_; i++)
722  mgc(levelContext);
723  }
724 
725  moveToFineLevel(levelContext, processNextLevel);
726 #else
727  *lhs=0;
728 #endif
729 
730  if(levelContext.matrix == matrices_->matrices().finest()) {
731  coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
732  if(!coarsesolverconverged)
733  DUNE_THROW(MathError, "Coarse solver did not converge");
734  }
735  // postsmoothing
736  postsmooth(levelContext, postSteps_);
737 
738  }
739  }
740 
741  template<class M, class X, class S, class PI, class A>
743 
744  // restrict residual to all levels
745  typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
746  typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
747  typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
748  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
749 
750  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
751  ++pinfo;
753  ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
754  }
755 
756  // pinfo is invalid, set to coarsest level
757  //pinfo = matrices_->parallelInformation().coarsest
758  // calculate correction for all levels
759  lhs = lhs_->finest();
760  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
761 
762  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
763  // presmoothing
764  *lhs=0;
765  smoother->apply(*lhs, *rhs);
766  }
767 
768  // Coarse level solve
769 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
771  pinfo->copyOwnerToAll(*rhs, *rhs);
772  solver_->apply(*lhs, *rhs, res);
773 
774  if(!res.converged)
775  DUNE_THROW(MathError, "Coarse solver did not converge");
776 #else
777  *lhs=0;
778 #endif
779  // Prologate and add up corrections from all levels
780  --pinfo;
781  --aggregates;
782 
783  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
785  ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
786  }
787  }
788 
789 
791  template<class M, class X, class S, class PI, class A>
792  void AMG<M,X,S,PI,A>::post(Domain& x)
793  {
794  DUNE_UNUSED_PARAMETER(x);
795  // Postprocess all smoothers
796  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
797  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
798  Iterator coarsest = smoothers_->coarsest();
799  Iterator smoother = smoothers_->finest();
800  DIterator lhs = lhs_->finest();
801  if(smoothers_->levels()>0) {
802  if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
803  smoother->post(*lhs);
804  if(smoother!=coarsest)
805  for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
806  smoother->post(*lhs);
807  smoother->post(*lhs);
808  }
809  //delete &(*lhs_->finest());
810  delete lhs_;
811  lhs_=nullptr;
812  //delete &(*update_->finest());
813  delete update_;
814  update_=nullptr;
815  //delete &(*rhs_->finest());
816  delete rhs_;
817  rhs_=nullptr;
818  }
819 
820  template<class M, class X, class S, class PI, class A>
821  template<class A1>
822  void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
823  {
824  matrices_->getCoarsestAggregatesOnFinest(cont);
825  }
826 
827  } // end namespace Amg
828 } // end namespace Dune
829 
830 #endif
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:75
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:430
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:316
All parameters for AMG.
Definition: parameters.hh:390
Col col
Definition: matrixmatrix.hh:347
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:63
X Range
The range type.
Definition: amg.hh:80
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:194
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:161
A generic dynamic dense matrix.
Definition: matrix.hh:554
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:260
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:91
Statistics about the application of an inverse operator.
Definition: solver.hh:31
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:26
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:202
Implementations of the inverse operator interface.
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:55
Templates characterizing the type of a solver.
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:455
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
T block_type
Export the type representing the components.
Definition: matrix.hh:562
Definition: solvertype.hh:13
void post(Domain &x)
Clean up.
Definition: amg.hh:792
Prolongation and restriction for amg.
Two grid operator for AMG with Krylov cycle.
Definition: amg.hh:44
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:210
S Smoother
The type of the smoother.
Definition: amg.hh:88
Traits class for generically constructing non default constructable types.
Definition: novlpschwarz.hh:328
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:218
Define base class for scalar product and norm.
Definition: pinfo.hh:25
~AMG()
Definition: amg.hh:354
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:206
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:547
Matrix & mat
Definition: matrixmatrix.hh:343
Iterator finest()
Get an iterator positioned at the finest level.
Definition: hierarchy.hh:1379
Smoother SmootherType
Definition: amg.hh:190
Definition: basearray.hh:19
The solver category.
Definition: amg.hh:95
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Definition: example.cc:35
std::size_t level
The level index.
Definition: amg.hh:226
Choose the approriate scalar product for a solver category.
Definition: scalarproducts.hh:76
Definition: transfer.hh:30
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:214
std::size_t maxlevels()
Definition: amg.hh:540
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:198
Classes for the generic construction and application of the smoothers.
ConstIterator class for sequential access.
Definition: matrix.hh:397
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:822
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:408
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:257
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:73
X Domain
The domain type.
Definition: amg.hh:78
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:222
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:82
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:313
M Operator
The matrix operator type.
Definition: amg.hh:64
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:683
Classes for using UMFPack with ISTL matrices.
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1385
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:574
an algebraic multigrid method using a Krylov-cycle.
Definition: amg.hh:41
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
Classes for using SuperLU with ISTL matrices.
std::size_t levels()
Definition: amg.hh:535
Provides a classes representing the hierarchies in AMG.
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:71