dune-istl  2.2.1
amg.hh
Go to the documentation of this file.
1 // $Id: amg.hh 1762 2013-01-19 11:21:36Z mblatt $
2 #ifndef DUNE_AMG_AMG_HH
3 #define DUNE_AMG_AMG_HH
4 
5 #include<memory>
6 #include<dune/common/exceptions.hh>
10 #include<dune/istl/solvers.hh>
12 #include<dune/istl/superlu.hh>
14 #include<dune/common/typetraits.hh>
15 #include<dune/common/exceptions.hh>
16 
17 namespace Dune
18 {
19  namespace Amg
20  {
38  template<class M, class X, class S, class P, class K, class A>
39  class KAMG;
40 
41  template<class T>
42  class KAmgTwoGrid;
43 
51  template<class M, class X, class S, class PI=SequentialInformation,
52  class A=std::allocator<X> >
53  class AMG : public Preconditioner<X,X>
54  {
55  template<class M1, class X1, class S1, class P1, class K1, class A1>
56  friend class KAMG;
57 
58  friend class KAmgTwoGrid<AMG>;
59 
60  public:
62  typedef M Operator;
69  typedef PI ParallelInformation;
74 
76  typedef X Domain;
78  typedef X Range;
86  typedef S Smoother;
87 
90 
91  enum {
93  category = S::category
94  };
95 
111  AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
112  const SmootherArgs& smootherArgs, std::size_t gamma,
113  std::size_t preSmoothingSteps,
114  std::size_t postSmoothingSteps,
115  bool additive=false) DUNE_DEPRECATED;
116 
126  AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
127  const SmootherArgs& smootherArgs, const Parameters& parms);
128 
147  template<class C>
148  AMG(const Operator& fineOperator, const C& criterion,
149  const SmootherArgs& smootherArgs, std::size_t gamma,
150  std::size_t preSmoothingSteps,
151  std::size_t postSmoothingSteps,
152  bool additive=false,
153  const ParallelInformation& pinfo=ParallelInformation()) DUNE_DEPRECATED;
154 
166  template<class C>
167  AMG(const Operator& fineOperator, const C& criterion,
168  const SmootherArgs& smootherArgs,
170 
171  ~AMG();
172 
174  void pre(Domain& x, Range& b);
175 
177  void apply(Domain& v, const Range& d);
178 
180  void post(Domain& x);
181 
186  template<class A1>
187  void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
188 
189  std::size_t levels();
190 
191  std::size_t maxlevels();
192 
202  {
203  matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
204  }
205 
210  bool usesDirectCoarseLevelSolver() const;
211 
212  private:
214  void mgc();
215 
216  typename Hierarchy<Smoother,A>::Iterator smoother;
219  typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
220  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
221  typename Hierarchy<Domain,A>::Iterator lhs;
222  typename Hierarchy<Domain,A>::Iterator update;
223  typename Hierarchy<Range,A>::Iterator rhs;
224 
225  void additiveMgc();
226 
228  void presmooth();
229 
231  void postsmooth();
232 
236  void moveToFineLevel(bool processedFineLevel);
237 
239  bool moveToCoarseLevel();
240 
242  void initIteratorsWithFineLevel();
243 
245  OperatorHierarchy* matrices_;
247  SmootherArgs smootherArgs_;
249  Hierarchy<Smoother,A> smoothers_;
251  CoarseSolver* solver_;
253  Hierarchy<Range,A>* rhs_;
255  Hierarchy<Domain,A>* lhs_;
257  Hierarchy<Domain,A>* update_;
261  typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
263  ScalarProduct* scalarProduct_;
265  std::size_t gamma_;
267  std::size_t preSteps_;
269  std::size_t postSteps_;
270  std::size_t level;
271  bool buildHierarchy_;
272  bool additive;
273  bool coarsesolverconverged;
274  Smoother *coarseSmoother_;
276  std::size_t verbosity_;
277  };
278 
279  template<class M, class X, class S, class PI, class A>
280  AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
281  const SmootherArgs& smootherArgs,
282  std::size_t gamma, std::size_t preSmoothingSteps,
283  std::size_t postSmoothingSteps, bool additive_)
284  : matrices_(&matrices), smootherArgs_(smootherArgs),
285  smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
286  gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
287  additive(additive_), coarsesolverconverged(true),
288  coarseSmoother_(), verbosity_(2)
289  {
290  assert(matrices_->isBuilt());
291 
292  // build the necessary smoother hierarchies
293  matrices_->coarsenSmoother(smoothers_, smootherArgs_);
294  }
295 
296  template<class M, class X, class S, class PI, class A>
297  AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
298  const SmootherArgs& smootherArgs,
299  const Parameters& parms)
300  : matrices_(&matrices), smootherArgs_(smootherArgs),
301  smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
302  gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
303  postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
304  additive(parms.getAdditive()), coarsesolverconverged(true),
305  coarseSmoother_(), verbosity_(parms.debugLevel())
306  {
307  assert(matrices_->isBuilt());
308 
309  // build the necessary smoother hierarchies
310  matrices_->coarsenSmoother(smoothers_, smootherArgs_);
311  }
312 
313  template<class M, class X, class S, class PI, class A>
314  template<class C>
316  const C& criterion,
317  const SmootherArgs& smootherArgs,
318  std::size_t gamma, std::size_t preSmoothingSteps,
319  std::size_t postSmoothingSteps,
320  bool additive_,
321  const PI& pinfo)
322  : smootherArgs_(smootherArgs),
323  smoothers_(), solver_(), scalarProduct_(0), gamma_(gamma),
324  preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
325  additive(additive_), coarsesolverconverged(true),
326  coarseSmoother_(), verbosity_(criterion.debugLevel())
327  {
328  dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
329  "Matrix and Solver must match in terms of category!");
330  // TODO: reestablish compile time checks.
331  //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
332  // "Matrix and Solver must match in terms of category!");
333  Timer watch;
334  matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
335 
336  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
337 
338  // build the necessary smoother hierarchies
339  matrices_->coarsenSmoother(smoothers_, smootherArgs_);
340 
341  if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
342  std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
343  }
344 
345  template<class M, class X, class S, class PI, class A>
346  template<class C>
348  const C& criterion,
349  const SmootherArgs& smootherArgs,
350  const PI& pinfo)
351  : smootherArgs_(smootherArgs),
352  smoothers_(), solver_(), scalarProduct_(0),
353  gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
354  postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
355  additive(criterion.getAdditive()), coarsesolverconverged(true),
356  coarseSmoother_(), verbosity_(criterion.debugLevel())
357  {
358  dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
359  "Matrix and Solver must match in terms of category!");
360  // TODO: reestablish compile time checks.
361  //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
362  // "Matrix and Solver must match in terms of category!");
363  Timer watch;
364  matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
365 
366  matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
367 
368  // build the necessary smoother hierarchies
369  matrices_->coarsenSmoother(smoothers_, smootherArgs_);
370 
371  if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
372  std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
373  }
374 
375  template<class M, class X, class S, class PI, class A>
377  {
378  if(buildHierarchy_){
379  delete matrices_;
380  }
381  if(scalarProduct_)
382  delete scalarProduct_;
383  }
384 
385 
386  template<class M, class X, class S, class PI, class A>
388  {
389  // Detect Matrix rows where all offdiagonal entries are
390  // zero and set x such that A_dd*x_d=b_d
391  // Thus users can be more careless when setting up their linear
392  // systems.
393  typedef typename M::matrix_type Matrix;
394  typedef typename Matrix::ConstRowIterator RowIter;
395  typedef typename Matrix::ConstColIterator ColIter;
396  typedef typename Matrix::block_type Block;
397  Block zero;
398  zero=typename Matrix::field_type();
399 
400  const Matrix& mat=matrices_->matrices().finest()->getmat();
401  for(RowIter row=mat.begin(); row!=mat.end(); ++row){
402  bool isDirichlet = true;
403  bool hasDiagonal = false;
404  Block diagonal;
405  for(ColIter col=row->begin(); col!=row->end(); ++col){
406  if(row.index()==col.index()){
407  diagonal = *col;
408  hasDiagonal = false;
409  }else{
410  if(*col!=zero)
411  isDirichlet = false;
412  }
413  }
414  if(isDirichlet && hasDiagonal)
415  diagonal.solve(x[row.index()], b[row.index()]);
416  }
417 
418  if(smoothers_.levels()>0)
419  smoothers_.finest()->pre(x,b);
420  else
421  // No smoother to make x consistent! Do it by hand
422  matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
423  Range* copy = new Range(b);
424  rhs_ = new Hierarchy<Range,A>(*copy);
425  Domain* dcopy = new Domain(x);
426  lhs_ = new Hierarchy<Domain,A>(*dcopy);
427  dcopy = new Domain(x);
428  update_ = new Hierarchy<Domain,A>(*dcopy);
429  matrices_->coarsenVector(*rhs_);
430  matrices_->coarsenVector(*lhs_);
431  matrices_->coarsenVector(*update_);
432 
433  // Preprocess all smoothers
434  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
435  typedef typename Hierarchy<Range,A>::Iterator RIterator;
436  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
437  Iterator coarsest = smoothers_.coarsest();
438  Iterator smoother = smoothers_.finest();
439  RIterator rhs = rhs_->finest();
440  DIterator lhs = lhs_->finest();
441  if(smoothers_.levels()>0){
442 
443  assert(lhs_->levels()==rhs_->levels());
444  assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
445  assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
446 
447  if(smoother!=coarsest)
448  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
449  smoother->pre(*lhs,*rhs);
450  smoother->pre(*lhs,*rhs);
451  }
452 
453 
454  // The preconditioner might change x and b. So we have to
455  // copy the changes to the original vectors.
456  x = *lhs_->finest();
457  b = *rhs_->finest();
458 
459  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){
460  // We have the carsest level. Create the coarse Solver
461  SmootherArgs sargs(smootherArgs_);
462  sargs.iterations = 1;
463 
465  cargs.setArgs(sargs);
466  if(matrices_->redistributeInformation().back().isSetup()){
467  // Solve on the redistributed partitioning
468  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
469  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
470 
471  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
472  scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
473  }else{
474  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
475  cargs.setComm(*matrices_->parallelInformation().coarsest());
476 
477  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
478  scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
479  }
480 #if HAVE_SUPERLU
481  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
482  if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode
483  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
484  || (matrices_->parallelInformation().coarsest().isRedistributed()
485  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
486  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){ // redistribute and 1 proc
487  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
488  std::cout<<"Using superlu"<<std::endl;
489  if(matrices_->parallelInformation().coarsest().isRedistributed())
490  {
491  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
492  // We are still participating on this level
493  solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat());
494  else
495  solver_ = 0;
496  }else
497  solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat());
498  }else
499 #endif
500  {
501  if(matrices_->parallelInformation().coarsest().isRedistributed())
502  {
503  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
504  // We are still participating on this level
505  solver_ = new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
506  *scalarProduct_,
507  *coarseSmoother_, 1E-2, 10000, 0);
508  else
509  solver_ = 0;
510  }else
511  solver_ = new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
512  *scalarProduct_,
513  *coarseSmoother_, 1E-2, 1000, 0);
514  }
515  }
516  }
517  template<class M, class X, class S, class PI, class A>
519  {
520  return matrices_->levels();
521  }
522  template<class M, class X, class S, class PI, class A>
524  {
525  return matrices_->maxlevels();
526  }
527 
529  template<class M, class X, class S, class PI, class A>
531  {
532  if(additive){
533  *(rhs_->finest())=d;
534  additiveMgc();
535  v=*lhs_->finest();
536  }else{
537  // Init all iterators for the current level
538  initIteratorsWithFineLevel();
539 
540 
541  *lhs = v;
542  *rhs = d;
543  *update=0;
544  level=0;
545 
546  mgc();
547 
548  if(postSteps_==0||matrices_->maxlevels()==1)
549  pinfo->copyOwnerToAll(*update, *update);
550 
551  v=*update;
552  }
553 
554  }
555 
556  template<class M, class X, class S, class PI, class A>
558  {
559  smoother = smoothers_.finest();
560  matrix = matrices_->matrices().finest();
561  pinfo = matrices_->parallelInformation().finest();
562  redist =
563  matrices_->redistributeInformation().begin();
564  aggregates = matrices_->aggregatesMaps().begin();
565  lhs = lhs_->finest();
566  update = update_->finest();
567  rhs = rhs_->finest();
568  }
569 
570  template<class M, class X, class S, class PI, class A>
571  bool AMG<M,X,S,PI,A>
572  ::moveToCoarseLevel()
573  {
574 
575  bool processNextLevel=true;
576 
577  if(redist->isSetup()){
578  redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed());
579  processNextLevel =rhs.getRedistributed().size()>0;
580  if(processNextLevel){
581  //restrict defect to coarse level right hand side.
582  typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
583  ++pinfo;
585  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo);
586  }
587  }else{
588  //restrict defect to coarse level right hand side.
589  typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
590  ++pinfo;
592  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
593  }
594 
595  if(processNextLevel){
596  // prepare coarse system
597  ++lhs;
598  ++update;
599  ++matrix;
600  ++level;
601  ++redist;
602 
603  if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
604  // next level is not the globally coarsest one
605  ++smoother;
606  ++aggregates;
607  }
608  // prepare the update on the next level
609  *update=0;
610  }
611  return processNextLevel;
612  }
613 
614  template<class M, class X, class S, class PI, class A>
615  void AMG<M,X,S,PI,A>
616  ::moveToFineLevel(bool processNextLevel)
617  {
618  if(processNextLevel){
619  if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
620  // previous level is not the globally coarsest one
621  --smoother;
622  --aggregates;
623  }
624  --redist;
625  --level;
626  //prolongate and add the correction (update is in coarse left hand side)
627  --matrix;
628 
629  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
630  --lhs;
631  --pinfo;
632  }
633  if(redist->isSetup()){
634  // Need to redistribute during prolongate
635  lhs.getRedistributed()=0;
637  ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(),
638  *pinfo, *redist);
639  }else{
640  *lhs=0;
642  ::prolongate(*(*aggregates), *update, *lhs,
643  matrices_->getProlongationDampingFactor(), *pinfo);
644  }
645 
646 
647  if(processNextLevel){
648  --update;
649  --rhs;
650  }
651 
652  *update += *lhs;
653  }
654 
655 
656  template<class M, class X, class S, class PI, class A>
657  void AMG<M,X,S,PI,A>
658  ::presmooth()
659  {
660 
661  for(std::size_t i=0; i < preSteps_; ++i){
662  *lhs=0;
663  SmootherApplier<S>::preSmooth(*smoother, *lhs, *rhs);
664  // Accumulate update
665  *update += *lhs;
666 
667  // update defect
668  matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
669  pinfo->project(*rhs);
670  }
671  }
672 
673  template<class M, class X, class S, class PI, class A>
674  void AMG<M,X,S,PI,A>
675  ::postsmooth()
676  {
677 
678  for(std::size_t i=0; i < postSteps_; ++i){
679  // update defect
680  matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
681  *lhs=0;
682  pinfo->project(*rhs);
683  SmootherApplier<S>::postSmooth(*smoother, *lhs, *rhs);
684  // Accumulate update
685  *update += *lhs;
686  }
687  }
688 
689 
690  template<class M, class X, class S, class PI, class A>
692  {
694  }
695 
696  template<class M, class X, class S, class PI, class A>
697  void AMG<M,X,S,PI,A>::mgc(){
698  if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){
699  // Solve directly
701  res.converged=true; // If we do not compute this flag will not get updated
702  if(redist->isSetup()){
703  redist->redistribute(*rhs, rhs.getRedistributed());
704  if(rhs.getRedistributed().size()>0){
705  // We are still participating in the computation
706  pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed());
707  solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res);
708  }
709  redist->redistributeBackward(*update, update.getRedistributed());
710  pinfo->copyOwnerToAll(*update, *update);
711  }else{
712  pinfo->copyOwnerToAll(*rhs, *rhs);
713  solver_->apply(*update, *rhs, res);
714  }
715 
716  if (!res.converged)
717  coarsesolverconverged = false;
718  }else{
719  // presmoothing
720  presmooth();
721 
722 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
723  bool processNextLevel = moveToCoarseLevel();
724 
725  if(processNextLevel){
726  // next level
727  for(std::size_t i=0; i<gamma_; i++)
728  mgc();
729  }
730 
731  moveToFineLevel(processNextLevel);
732 #else
733  *lhs=0;
734 #endif
735 
736  if(matrix == matrices_->matrices().finest()){
737  coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
738  if(!coarsesolverconverged)
739  DUNE_THROW(MathError, "Coarse solver did not converge");
740  }
741  // postsmoothing
742  postsmooth();
743 
744  }
745  }
746 
747  template<class M, class X, class S, class PI, class A>
748  void AMG<M,X,S,PI,A>::additiveMgc(){
749 
750  // restrict residual to all levels
751  typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
752  typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
753  typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
754  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
755 
756  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates){
757  ++pinfo;
759  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
760  }
761 
762  // pinfo is invalid, set to coarsest level
763  //pinfo = matrices_->parallelInformation().coarsest
764  // calculate correction for all levels
765  lhs = lhs_->finest();
766  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest();
767 
768  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){
769  // presmoothing
770  *lhs=0;
771  smoother->apply(*lhs, *rhs);
772  }
773 
774  // Coarse level solve
775 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
776  InverseOperatorResult res;
777  pinfo->copyOwnerToAll(*rhs, *rhs);
778  solver_->apply(*lhs, *rhs, res);
779 
780  if(!res.converged)
781  DUNE_THROW(MathError, "Coarse solver did not converge");
782 #else
783  *lhs=0;
784 #endif
785  // Prologate and add up corrections from all levels
786  --pinfo;
787  --aggregates;
788 
789  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo){
791  ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo);
792  }
793  }
794 
795 
797  template<class M, class X, class S, class PI, class A>
799  {
800  if(buildHierarchy_){
801  if(solver_)
802  delete solver_;
803  if(coarseSmoother_)
805  }
806 
807  // Postprocess all smoothers
808  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
809  typedef typename Hierarchy<Range,A>::Iterator RIterator;
810  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
811  Iterator coarsest = smoothers_.coarsest();
812  Iterator smoother = smoothers_.finest();
813  DIterator lhs = lhs_->finest();
814  if(smoothers_.levels()>0){
815  if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
816  smoother->post(*lhs);
817  if(smoother!=coarsest)
818  for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
819  smoother->post(*lhs);
820  smoother->post(*lhs);
821  }
822 
823  delete &(*lhs_->finest());
824  delete lhs_;
825  delete &(*update_->finest());
826  delete update_;
827  delete &(*rhs_->finest());
828  delete rhs_;
829  }
830 
831  template<class M, class X, class S, class PI, class A>
832  template<class A1>
833  void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
834  {
835  matrices_->getCoarsestAggregatesOnFinest(cont);
836  }
837 
838  } // end namespace Amg
839 }// end namespace Dune
840 
841 #endif