My Project
amgcpr.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_CPR_HH
4 #define DUNE_AMG_AMG_CPR_HH
5 
6 // NOTE: This file is a modified version of dune/istl/paamg/amg.hh from
7 // dune-istl release 2.6.0. Modifications have been kept as minimal as possible.
8 
9 #include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
10 
11 #include <dune/common/exceptions.hh>
12 #include <dune/istl/paamg/smoother.hh>
13 #include <dune/istl/paamg/transfer.hh>
14 #include <dune/istl/paamg/hierarchy.hh>
15 #include <dune/istl/solvers.hh>
16 #include <dune/istl/scalarproducts.hh>
17 #include <dune/istl/superlu.hh>
18 #include <dune/istl/umfpack.hh>
19 #include <dune/istl/solvertype.hh>
20 #include <dune/common/typetraits.hh>
21 #include <dune/common/exceptions.hh>
22 
23 #include <memory>
24 
25 namespace Dune
26 {
27  namespace Amg
28  {
29 
30 
31 #if HAVE_MPI
32  template<class M, class T>
33  void redistributeMatrixAmg(M&, M&, SequentialInformation&, SequentialInformation&, T&)
34  {
35  // noop
36  }
37 
38  template<class M, class PI>
39  typename std::enable_if<!std::is_same<PI,SequentialInformation>::value,void>::type
40  redistributeMatrixAmg(M& mat, M& matRedist, PI& info, PI& infoRedist,
41  Dune::RedistributeInformation<PI>& redistInfo)
42  {
43  info.buildGlobalLookup(mat.N());
44  redistributeMatrixEntries(mat, matRedist, info, infoRedist, redistInfo);
45  info.freeGlobalLookup();
46  }
47 #endif
48 
66  template<class M, class X, class S, class P, class K, class A>
67  class KAMG;
68 
69  template<class T>
70  class KAmgTwoGrid;
71 
82  template<class M, class X, class S, class PI=SequentialInformation,
83  class A=std::allocator<X> >
84  class AMGCPR : public PreconditionerWithUpdate<X,X>
85  {
86  template<class M1, class X1, class S1, class P1, class K1, class A1>
87  friend class KAMG;
88 
89  friend class KAmgTwoGrid<AMGCPR>;
90 
91  public:
93  typedef M Operator;
102  typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
104  typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
105 
107  typedef X Domain;
109  typedef X Range;
111  typedef InverseOperator<X,X> CoarseSolver;
117  typedef S Smoother;
118 
120  typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
121 
131  AMGCPR(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
132  const SmootherArgs& smootherArgs, const Parameters& parms);
133 
145  template<class C>
146  AMGCPR(const Operator& fineOperator, const C& criterion,
147  const SmootherArgs& smootherArgs=SmootherArgs(),
149 
153  AMGCPR(const AMGCPR& amg);
154 
155  ~AMGCPR();
156 
158  void pre(Domain& x, Range& b);
159 
161  void apply(Domain& v, const Range& d);
162 
164  virtual SolverCategory::Category category() const
165  {
166  return category_;
167  }
168 
170  void post(Domain& x);
171 
176  template<class A1>
177  void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
178 
179  std::size_t levels();
180 
181  std::size_t maxlevels();
182 
192  {
193  auto copyFlags = NegateSet<typename PI::OwnerSet>();
194  const auto& matrices = matrices_->matrices();
195  const auto& aggregatesMapHierarchy = matrices_->aggregatesMaps();
196  const auto& infoHierarchy = matrices_->parallelInformation();
197  const auto& redistInfoHierarchy = matrices_->redistributeInformation();
198  BaseGalerkinProduct productBuilder;
199  auto aggregatesMap = aggregatesMapHierarchy.begin();
200  auto info = infoHierarchy.finest();
201  auto redistInfo = redistInfoHierarchy.begin();
202  auto matrix = matrices.finest();
203  auto coarsestMatrix = matrices.coarsest();
204 
205  using Matrix = typename M::matrix_type;
206 
207 #if HAVE_MPI
208  if(matrix.isRedistributed()) {
209  redistributeMatrixAmg(const_cast<Matrix&>(matrix->getmat()),
210  const_cast<Matrix&>(matrix.getRedistributed().getmat()),
211  const_cast<PI&>(*info), const_cast<PI&>(info.getRedistributed()),
212  const_cast<Dune::RedistributeInformation<PI>&>(*redistInfo));
213  }
214 #endif
215 
216  for(; matrix!=coarsestMatrix; ++aggregatesMap) {
217  const Matrix& fine = (matrix.isRedistributed() ? matrix.getRedistributed() : *matrix).getmat();
218  ++matrix;
219  ++info;
220  ++redistInfo;
221  productBuilder.calculate(fine, *(*aggregatesMap), const_cast<Matrix&>(matrix->getmat()), *info, copyFlags);
222 #if HAVE_MPI
223  if(matrix.isRedistributed()) {
224  redistributeMatrixAmg(const_cast<Matrix&>(matrix->getmat()),
225  const_cast<Matrix&>(matrix.getRedistributed().getmat()),
226  const_cast<PI&>(*info), const_cast<PI&>(info.getRedistributed()),
227  const_cast<Dune::RedistributeInformation<PI>&>(*redistInfo));
228  }
229 #endif
230  }
231  }
232 
236  template<class C>
237  void updateSolver(C& criterion, const Operator& /* matrix */, const PI& pinfo);
238 
242  virtual void update();
243 
248  bool usesDirectCoarseLevelSolver() const;
249 
250  private:
257 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7 )
258  template<class C>
259  void createHierarchies(C& criterion, Operator& matrix,
260  const PI& pinfo)
261  {
262  // create shared_ptr with empty deleter
263  std::shared_ptr< Operator > op( &matrix, []( Operator* ) {});
264  std::shared_ptr< PI > pifo( const_cast< PI* > (&pinfo), []( PI * ) {});
265  createHierarchies( criterion, op, pifo);
266  }
267 
268  template<class C>
269  void createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
270  std::shared_ptr< PI > pinfo );
271 
272 #else
273  template<class C>
274  void createHierarchies(C& criterion, Operator& matrix,
275  const PI& pinfo);
276 #endif
277 
278  void setupCoarseSolver();
279 
286  struct LevelContext
287  {
288  typedef Smoother SmootherType;
292  typename Hierarchy<Smoother,A>::Iterator smoother;
296  typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
300  typename ParallelInformationHierarchy::Iterator pinfo;
304  typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
308  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
312  typename Hierarchy<Domain,A>::Iterator lhs;
316  typename Hierarchy<Domain,A>::Iterator update;
320  typename Hierarchy<Range,A>::Iterator rhs;
324  std::size_t level;
325  };
326 
327 
332  void mgc(LevelContext& levelContext);
333 
334  void additiveMgc();
335 
342  void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
343 
348  bool moveToCoarseLevel(LevelContext& levelContext);
349 
354  void initIteratorsWithFineLevel(LevelContext& levelContext);
355 
357  std::shared_ptr<OperatorHierarchy> matrices_;
359  SmootherArgs smootherArgs_;
361  std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
363  std::shared_ptr<CoarseSolver> solver_;
365  std::shared_ptr< Hierarchy<Range,A> > rhs_;
367  std::shared_ptr< Hierarchy<Domain,A> > lhs_;
369  std::shared_ptr< Hierarchy<Domain,A> > update_;
371  using ScalarProduct = Dune::ScalarProduct<X>;
373  std::shared_ptr<ScalarProduct> scalarProduct_;
375  std::size_t gamma_;
377  std::size_t preSteps_;
379  std::size_t postSteps_;
380  bool buildHierarchy_;
381  bool additive;
382  bool coarsesolverconverged;
383  std::shared_ptr<Smoother> coarseSmoother_;
385  SolverCategory::Category category_;
387  std::size_t verbosity_;
388  };
389 
390  template<class M, class X, class S, class PI, class A>
392  : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
393  smoothers_(amg.smoothers_), solver_(amg.solver_),
394  rhs_(), lhs_(), update_(),
395  scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
396  preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
397  buildHierarchy_(amg.buildHierarchy_),
398  additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
399  coarseSmoother_(amg.coarseSmoother_),
400  category_(amg.category_),
401  verbosity_(amg.verbosity_)
402  {
403  if(amg.rhs_)
404  rhs_.reset( new Hierarchy<Range,A>(*amg.rhs_) );
405  if(amg.lhs_)
406  lhs_.reset( new Hierarchy<Domain,A>(*amg.lhs_) );
407  if(amg.update_)
408  update_.reset( new Hierarchy<Domain,A>(*amg.update_) );
409  }
410 
411  template<class M, class X, class S, class PI, class A>
413  const SmootherArgs& smootherArgs,
414  const Parameters& parms)
415  : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
416  smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
417  rhs_(), lhs_(), update_(), scalarProduct_(0),
418  gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
419  postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
420  additive(parms.getAdditive()), coarsesolverconverged(true),
421  coarseSmoother_(),
422 // #warning should category be retrieved from matrices?
423  category_(SolverCategory::category(*smoothers_->coarsest())),
424  verbosity_(parms.debugLevel())
425  {
426  assert(matrices_->isBuilt());
427 
428  // build the necessary smoother hierarchies
429  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
430  }
431 
432  template<class M, class X, class S, class PI, class A>
433  template<class C>
435  const C& criterion,
436  const SmootherArgs& smootherArgs,
437  const PI& pinfo)
438  : smootherArgs_(smootherArgs),
439  smoothers_(new Hierarchy<Smoother,A>), solver_(),
440  rhs_(), lhs_(), update_(), scalarProduct_(),
441  gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
442  postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
443  additive(criterion.getAdditive()), coarsesolverconverged(true),
444  coarseSmoother_(),
445  category_(SolverCategory::category(pinfo)),
446  verbosity_(criterion.debugLevel())
447  {
448  if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
449  DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
450  createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
451  }
452 
453  template<class M, class X, class S, class PI, class A>
455  {
456  if(buildHierarchy_) {
457  if(solver_)
458  solver_.reset();
459  if(coarseSmoother_)
460  coarseSmoother_.reset();
461  }
462  }
463 
464  template<class M, class X, class S, class PI, class A>
465  template<class C>
466  void AMGCPR<M,X,S,PI,A>::updateSolver(C& /* criterion */, const Operator& /* matrix */, const PI& /* pinfo */)
467  {
468  update();
469  }
470 
471  template<class M, class X, class S, class PI, class A>
473  {
474  Timer watch;
475  smoothers_.reset(new Hierarchy<Smoother,A>);
476  solver_.reset();
477  coarseSmoother_.reset();
478  scalarProduct_.reset();
479  buildHierarchy_= true;
480  coarsesolverconverged = true;
481  smoothers_.reset(new Hierarchy<Smoother,A>);
482  recalculateHierarchy();
483  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
484  setupCoarseSolver();
485  if (verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) {
486  std::cout << "Recalculating galerkin and coarse smoothers "<< matrices_->maxlevels() << " levels "
487  << watch.elapsed() << " seconds." << std::endl;
488  }
489  }
490 
491  template<class M, class X, class S, class PI, class A>
492  template<class C>
493 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
494  void AMGCPR<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr< Operator > matrix,
495  std::shared_ptr< PI > pinfo )
496 #else
497  void AMGCPR<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix, const PI& pinfo )
498 #endif
499  {
500  Timer watch;
501  matrices_.reset(new OperatorHierarchy(matrix, pinfo));
502 
503  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
504 
505  // build the necessary smoother hierarchies
506  matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
507  setupCoarseSolver();
508  if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
509  std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
510  <<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
511  }
512 
513  template<class M, class X, class S, class PI, class A>
514  void AMGCPR<M,X,S,PI,A>::setupCoarseSolver()
515  {
516  // test whether we should solve on the coarse level. That is the case if we
517  // have that level and if there was a redistribution on this level then our
518  // communicator has to be valid (size()>0) as the smoother might try to communicate
519  // in the constructor.
520  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
521  && ( ! matrices_->redistributeInformation().back().isSetup() ||
522  matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
523  {
524  // We have the carsest level. Create the coarse Solver
525  SmootherArgs sargs(smootherArgs_);
526  sargs.iterations = 1;
527 
528  typename ConstructionTraits<Smoother>::Arguments cargs;
529  cargs.setArgs(sargs);
530  if(matrices_->redistributeInformation().back().isSetup()) {
531  // Solve on the redistributed partitioning
532  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
533  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
534  }else{
535  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
536  cargs.setComm(*matrices_->parallelInformation().coarsest());
537  }
538 
539 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
540  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
541 #else
542  coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
543 #endif
544 
545  scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
546 
547 
548  typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
549 
550  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
551  if( SolverSelector::isDirectSolver &&
552  (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
553  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
554  || (matrices_->parallelInformation().coarsest().isRedistributed()
555  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
556  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
557  )
558  { // redistribute and 1 proc
559  if(matrices_->parallelInformation().coarsest().isRedistributed())
560  {
561  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
562  {
563  // We are still participating on this level
564  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
565  }
566  else
567  solver_.reset();
568  }
569  else
570  {
571  solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
572  }
573  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
574  std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
575  }
576  else
577  {
578  if(matrices_->parallelInformation().coarsest().isRedistributed())
579  {
580  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
581  // We are still participating on this level
582 
583  // we have to allocate these types using the rebound allocator
584  // in order to ensure that we fulfill the alignement requirements
585  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
586  // Cast needed for Dune <=2.5
587  reinterpret_cast<typename
588  std::conditional<std::is_same<PI,SequentialInformation>::value,
589  Dune::SeqScalarProduct<X>,
590  Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
591  *coarseSmoother_, 1E-2, 1000, 0));
592  else
593  solver_.reset();
594  }else
595  {
596  solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
597  // Cast needed for Dune <=2.5
598  reinterpret_cast<typename
599  std::conditional<std::is_same<PI,SequentialInformation>::value,
600  Dune::SeqScalarProduct<X>,
601  Dune::OverlappingSchwarzScalarProduct<X,PI> >::type&>(*scalarProduct_),
602  *coarseSmoother_, 1E-2, 1000, 0));
603  // // we have to allocate these types using the rebound allocator
604  // // in order to ensure that we fulfill the alignement requirements
605  // using Alloc = typename A::template rebind<BiCGSTABSolver<X>>::other;
606  // Alloc alloc;
607  // auto p = alloc.allocate(1);
608  // alloc.construct(p,
609  // const_cast<M&>(*matrices_->matrices().coarsest()),
610  // *scalarProduct_,
611  // *coarseSmoother_, 1E-2, 1000, 0);
612  // solver_.reset(p,[](BiCGSTABSolver<X>* p){
613  // Alloc alloc;
614  // alloc.destroy(p);
615  // alloc.deallocate(p,1);
616  // });
617  }
618  }
619  }
620  }
621 
622  template<class M, class X, class S, class PI, class A>
624  {
625  // Detect Matrix rows where all offdiagonal entries are
626  // zero and set x such that A_dd*x_d=b_d
627  // Thus users can be more careless when setting up their linear
628  // systems.
629  typedef typename M::matrix_type Matrix;
630  typedef typename Matrix::ConstRowIterator RowIter;
631  typedef typename Matrix::ConstColIterator ColIter;
632  typedef typename Matrix::block_type Block;
633  Block zero;
634  zero=typename Matrix::field_type();
635 
636  const Matrix& mat=matrices_->matrices().finest()->getmat();
637  for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
638  bool isDirichlet = true;
639  bool hasDiagonal = false;
640  Block diagonal;
641  for(ColIter col=row->begin(); col!=row->end(); ++col) {
642  if(row.index()==col.index()) {
643  diagonal = *col;
644  hasDiagonal = false;
645  }else{
646  if(*col!=zero)
647  isDirichlet = false;
648  }
649  }
650  if(isDirichlet && hasDiagonal)
651  diagonal.solve(x[row.index()], b[row.index()]);
652  }
653 
654  if(smoothers_->levels()>0)
655  smoothers_->finest()->pre(x,b);
656  else
657  // No smoother to make x consistent! Do it by hand
658  matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
659 
660 
661 #if DUNE_VERSION_NEWER( DUNE_ISTL, 2, 7)
662  typedef std::shared_ptr< Range > RangePtr ;
663  typedef std::shared_ptr< Domain > DomainPtr;
664 #else
665  typedef Range* RangePtr;
666  typedef Domain* DomainPtr;
667 #endif
668 
669  // Hierarchy takes ownership of pointers
670  RangePtr copy( new Range(b) );
671  rhs_.reset( new Hierarchy<Range,A>(copy) );
672  DomainPtr dcopy( new Domain(x) );
673  lhs_.reset( new Hierarchy<Domain,A>(dcopy) );
674  DomainPtr dcopy2( new Domain(x) );
675  update_.reset( new Hierarchy<Domain,A>(dcopy2) );
676 
677  matrices_->coarsenVector(*rhs_);
678  matrices_->coarsenVector(*lhs_);
679  matrices_->coarsenVector(*update_);
680 
681  // Preprocess all smoothers
682  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
683  typedef typename Hierarchy<Range,A>::Iterator RIterator;
684  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
685  Iterator coarsest = smoothers_->coarsest();
686  Iterator smoother = smoothers_->finest();
687  RIterator rhs = rhs_->finest();
688  DIterator lhs = lhs_->finest();
689  if(smoothers_->levels()>0) {
690 
691  assert(lhs_->levels()==rhs_->levels());
692  assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
693  assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
694 
695  if(smoother!=coarsest)
696  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
697  smoother->pre(*lhs,*rhs);
698  smoother->pre(*lhs,*rhs);
699  }
700 
701 
702  // The preconditioner might change x and b. So we have to
703  // copy the changes to the original vectors.
704  x = *lhs_->finest();
705  b = *rhs_->finest();
706 
707  }
708  template<class M, class X, class S, class PI, class A>
709  std::size_t AMGCPR<M,X,S,PI,A>::levels()
710  {
711  return matrices_->levels();
712  }
713  template<class M, class X, class S, class PI, class A>
714  std::size_t AMGCPR<M,X,S,PI,A>::maxlevels()
715  {
716  return matrices_->maxlevels();
717  }
718 
720  template<class M, class X, class S, class PI, class A>
722  {
723  LevelContext levelContext;
724 
725  if(additive) {
726  *(rhs_->finest())=d;
727  additiveMgc();
728  v=*lhs_->finest();
729  }else{
730  // Init all iterators for the current level
731  initIteratorsWithFineLevel(levelContext);
732 
733 
734  *levelContext.lhs = v;
735  *levelContext.rhs = d;
736  *levelContext.update=0;
737  levelContext.level=0;
738 
739  mgc(levelContext);
740 
741  if(postSteps_==0||matrices_->maxlevels()==1)
742  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
743 
744  v=*levelContext.update;
745  }
746 
747  }
748 
749  template<class M, class X, class S, class PI, class A>
750  void AMGCPR<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
751  {
752  levelContext.smoother = smoothers_->finest();
753  levelContext.matrix = matrices_->matrices().finest();
754  levelContext.pinfo = matrices_->parallelInformation().finest();
755  levelContext.redist =
756  matrices_->redistributeInformation().begin();
757  levelContext.aggregates = matrices_->aggregatesMaps().begin();
758  levelContext.lhs = lhs_->finest();
759  levelContext.update = update_->finest();
760  levelContext.rhs = rhs_->finest();
761  }
762 
763  template<class M, class X, class S, class PI, class A>
764  bool AMGCPR<M,X,S,PI,A>
765  ::moveToCoarseLevel(LevelContext& levelContext)
766  {
767 
768  bool processNextLevel=true;
769 
770  if(levelContext.redist->isSetup()) {
771  levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
772  levelContext.rhs.getRedistributed());
773  processNextLevel = levelContext.rhs.getRedistributed().size()>0;
774  if(processNextLevel) {
775  //restrict defect to coarse level right hand side.
776  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
777  ++levelContext.pinfo;
778  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
779  ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
780  static_cast<const Range&>(fineRhs.getRedistributed()),
781  *levelContext.pinfo);
782  }
783  }else{
784  //restrict defect to coarse level right hand side.
785  typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
786  ++levelContext.pinfo;
787  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
788  ::restrictVector(*(*levelContext.aggregates),
789  *levelContext.rhs, static_cast<const Range&>(*fineRhs),
790  *levelContext.pinfo);
791  }
792 
793  if(processNextLevel) {
794  // prepare coarse system
795  ++levelContext.lhs;
796  ++levelContext.update;
797  ++levelContext.matrix;
798  ++levelContext.level;
799  ++levelContext.redist;
800 
801  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
802  // next level is not the globally coarsest one
803  ++levelContext.smoother;
804  ++levelContext.aggregates;
805  }
806  // prepare the update on the next level
807  *levelContext.update=0;
808  }
809  return processNextLevel;
810  }
811 
812  template<class M, class X, class S, class PI, class A>
813  void AMGCPR<M,X,S,PI,A>
814  ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
815  {
816  if(processNextLevel) {
817  if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
818  // previous level is not the globally coarsest one
819  --levelContext.smoother;
820  --levelContext.aggregates;
821  }
822  --levelContext.redist;
823  --levelContext.level;
824  //prolongate and add the correction (update is in coarse left hand side)
825  --levelContext.matrix;
826 
827  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
828  --levelContext.lhs;
829  --levelContext.pinfo;
830  }
831  if(levelContext.redist->isSetup()) {
832  // Need to redistribute during prolongateVector
833  levelContext.lhs.getRedistributed()=0;
834  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
835  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
836  levelContext.lhs.getRedistributed(),
837  matrices_->getProlongationDampingFactor(),
838  *levelContext.pinfo, *levelContext.redist);
839  }else{
840  *levelContext.lhs=0;
841  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
842  ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
843  matrices_->getProlongationDampingFactor(),
844  *levelContext.pinfo);
845  }
846 
847 
848  if(processNextLevel) {
849  --levelContext.update;
850  --levelContext.rhs;
851  }
852 
853  *levelContext.update += *levelContext.lhs;
854  }
855 
856  template<class M, class X, class S, class PI, class A>
858  {
859  return IsDirectSolver< CoarseSolver>::value;
860  }
861 
862  template<class M, class X, class S, class PI, class A>
863  void AMGCPR<M,X,S,PI,A>::mgc(LevelContext& levelContext){
864  if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
865  // Solve directly
866  InverseOperatorResult res;
867  res.converged=true; // If we do not compute this flag will not get updated
868  if(levelContext.redist->isSetup()) {
869  levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
870  if(levelContext.rhs.getRedistributed().size()>0) {
871  // We are still participating in the computation
872  levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
873  levelContext.rhs.getRedistributed());
874  solver_->apply(levelContext.update.getRedistributed(),
875  levelContext.rhs.getRedistributed(), res);
876  }
877  levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
878  levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
879  }else{
880  levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
881  solver_->apply(*levelContext.update, *levelContext.rhs, res);
882  }
883 
884  if (!res.converged)
885  coarsesolverconverged = false;
886  }else{
887  // presmoothing
888  presmooth(levelContext, preSteps_);
889 
890 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
891  bool processNextLevel = moveToCoarseLevel(levelContext);
892 
893  if(processNextLevel) {
894  // next level
895  for(std::size_t i=0; i<gamma_; i++)
896  mgc(levelContext);
897  }
898 
899  moveToFineLevel(levelContext, processNextLevel);
900 #else
901  *lhs=0;
902 #endif
903 
904  if(levelContext.matrix == matrices_->matrices().finest()) {
905  coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
906  if(!coarsesolverconverged){
907  //DUNE_THROW(MathError, "Coarse solver did not converge");
908  }
909  }
910  // postsmoothing
911  postsmooth(levelContext, postSteps_);
912 
913  }
914  }
915 
916  template<class M, class X, class S, class PI, class A>
917  void AMGCPR<M,X,S,PI,A>::additiveMgc(){
918 
919  // restrict residual to all levels
920  typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
921  typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
922  typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
923  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
924 
925  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
926  ++pinfo;
927  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
928  ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
929  }
930 
931  // pinfo is invalid, set to coarsest level
932  //pinfo = matrices_->parallelInformation().coarsest
933  // calculate correction for all levels
934  lhs = lhs_->finest();
935  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
936 
937  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
938  // presmoothing
939  *lhs=0;
940  smoother->apply(*lhs, *rhs);
941  }
942 
943  // Coarse level solve
944 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
945  InverseOperatorResult res;
946  pinfo->copyOwnerToAll(*rhs, *rhs);
947  solver_->apply(*lhs, *rhs, res);
948 
949  if(!res.converged)
950  DUNE_THROW(MathError, "Coarse solver did not converge");
951 #else
952  *lhs=0;
953 #endif
954  // Prologate and add up corrections from all levels
955  --pinfo;
956  --aggregates;
957 
958  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
959  Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
960  ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
961  }
962  }
963 
964 
966  template<class M, class X, class S, class PI, class A>
968  {
969  DUNE_UNUSED_PARAMETER(x);
970  // Postprocess all smoothers
971  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
972  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
973  Iterator coarsest = smoothers_->coarsest();
974  Iterator smoother = smoothers_->finest();
975  DIterator lhs = lhs_->finest();
976  if(smoothers_->levels()>0) {
977  if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
978  smoother->post(*lhs);
979  if(smoother!=coarsest)
980  for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
981  smoother->post(*lhs);
982  smoother->post(*lhs);
983  }
984  //delete &(*lhs_->finest());
985  lhs_.reset();
986  //delete &(*update_->finest());
987  update_.reset();
988  //delete &(*rhs_->finest());
989  rhs_.reset();
990  }
991 
992  template<class M, class X, class S, class PI, class A>
993  template<class A1>
994  void AMGCPR<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
995  {
996  matrices_->getCoarsestAggregatesOnFinest(cont);
997  }
998 
999  } // end namespace Amg
1000 } // end namespace Dune
1001 
1002 #endif
Parallel algebraic multigrid based on agglomeration.
Definition: amgcpr.hh:85
Definition: amgcpr.hh:67
Definition: amgcpr.hh:70
Interface class adding the update() method to the preconditioner interface.
Definition: PreconditionerWithUpdate.hpp:32
M Operator
The matrix operator type.
Definition: amgcpr.hh:93
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amgcpr.hh:292
AMGCPR(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amgcpr.hh:412
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amgcpr.hh:316
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amgcpr.hh:320
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amgcpr.hh:102
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amgcpr.hh:300
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amgcpr.hh:164
void apply(Domain &v, const Range &d)
Definition: amgcpr.hh:721
void pre(Domain &x, Range &b)
Definition: amgcpr.hh:623
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amgcpr.hh:296
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amgcpr.hh:994
void post(Domain &x)
Definition: amgcpr.hh:967
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amgcpr.hh:104
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amgcpr.hh:191
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amgcpr.hh:120
X Range
The range type.
Definition: amgcpr.hh:109
X Domain
The domain type.
Definition: amgcpr.hh:107
S Smoother
The type of the smoother.
Definition: amgcpr.hh:117
void updateSolver(C &criterion, const Operator &, const PI &pinfo)
Update the coarse solver and the hierarchies.
Definition: amgcpr.hh:466
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amgcpr.hh:308
virtual void update()
Update the coarse solver and the hierarchies.
Definition: amgcpr.hh:472
PI ParallelInformation
The type of the parallel information.
Definition: amgcpr.hh:100
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amgcpr.hh:111
std::size_t level
The level index.
Definition: amgcpr.hh:324
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amgcpr.hh:312
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amgcpr.hh:857
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amgcpr.hh:304