dune-istl  2.2.0
amg.hh
Go to the documentation of this file.
1 // $Id: amg.hh 1547 2012-04-19 17:19:40Z 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  if(smoothers_.levels()>0)
390  smoothers_.finest()->pre(x,b);
391  else
392  // No smoother to make x consistent! Do it by hand
393  matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
394  Range* copy = new Range(b);
395  rhs_ = new Hierarchy<Range,A>(*copy);
396  Domain* dcopy = new Domain(x);
397  lhs_ = new Hierarchy<Domain,A>(*dcopy);
398  dcopy = new Domain(x);
399  update_ = new Hierarchy<Domain,A>(*dcopy);
400  matrices_->coarsenVector(*rhs_);
401  matrices_->coarsenVector(*lhs_);
402  matrices_->coarsenVector(*update_);
403 
404  // Preprocess all smoothers
405  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
406  typedef typename Hierarchy<Range,A>::Iterator RIterator;
407  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
408  Iterator coarsest = smoothers_.coarsest();
409  Iterator smoother = smoothers_.finest();
410  RIterator rhs = rhs_->finest();
411  DIterator lhs = lhs_->finest();
412  if(smoothers_.levels()>0){
413 
414  assert(lhs_->levels()==rhs_->levels());
415  assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
416  assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
417 
418  if(smoother!=coarsest)
419  for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
420  smoother->pre(*lhs,*rhs);
421  smoother->pre(*lhs,*rhs);
422  }
423 
424 
425  // The preconditioner might change x and b. So we have to
426  // copy the changes to the original vectors.
427  x = *lhs_->finest();
428  b = *rhs_->finest();
429 
430  if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){
431  // We have the carsest level. Create the coarse Solver
432  SmootherArgs sargs(smootherArgs_);
433  sargs.iterations = 1;
434 
436  cargs.setArgs(sargs);
437  if(matrices_->redistributeInformation().back().isSetup()){
438  // Solve on the redistributed partitioning
439  cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
440  cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
441 
442  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
443  scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
444  }else{
445  cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
446  cargs.setComm(*matrices_->parallelInformation().coarsest());
447 
448  coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
449  scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
450  }
451 #if HAVE_SUPERLU
452  // Use superlu if we are purely sequential or with only one processor on the coarsest level.
453  if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode
454  || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
455  || (matrices_->parallelInformation().coarsest().isRedistributed()
456  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
457  && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){ // redistribute and 1 proc
458  if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
459  std::cout<<"Using superlu"<<std::endl;
460  if(matrices_->parallelInformation().coarsest().isRedistributed())
461  {
462  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
463  // We are still participating on this level
464  solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat());
465  else
466  solver_ = 0;
467  }else
468  solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat());
469  }else
470 #endif
471  {
472  if(matrices_->parallelInformation().coarsest().isRedistributed())
473  {
474  if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
475  // We are still participating on this level
476  solver_ = new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
477  *scalarProduct_,
478  *coarseSmoother_, 1E-2, 10000, 0);
479  else
480  solver_ = 0;
481  }else
482  solver_ = new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
483  *scalarProduct_,
484  *coarseSmoother_, 1E-2, 1000, 0);
485  }
486  }
487  }
488  template<class M, class X, class S, class PI, class A>
490  {
491  return matrices_->levels();
492  }
493  template<class M, class X, class S, class PI, class A>
495  {
496  return matrices_->maxlevels();
497  }
498 
500  template<class M, class X, class S, class PI, class A>
502  {
503  if(additive){
504  *(rhs_->finest())=d;
505  additiveMgc();
506  v=*lhs_->finest();
507  }else{
508  // Init all iterators for the current level
509  initIteratorsWithFineLevel();
510 
511 
512  *lhs = v;
513  *rhs = d;
514  *update=0;
515  level=0;
516 
517  mgc();
518 
519  if(postSteps_==0||matrices_->maxlevels()==1)
520  pinfo->copyOwnerToAll(*update, *update);
521 
522  v=*update;
523  }
524 
525  }
526 
527  template<class M, class X, class S, class PI, class A>
529  {
530  smoother = smoothers_.finest();
531  matrix = matrices_->matrices().finest();
532  pinfo = matrices_->parallelInformation().finest();
533  redist =
534  matrices_->redistributeInformation().begin();
535  aggregates = matrices_->aggregatesMaps().begin();
536  lhs = lhs_->finest();
537  update = update_->finest();
538  rhs = rhs_->finest();
539  }
540 
541  template<class M, class X, class S, class PI, class A>
542  bool AMG<M,X,S,PI,A>
543  ::moveToCoarseLevel()
544  {
545 
546  bool processNextLevel=true;
547 
548  if(redist->isSetup()){
549  redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed());
550  processNextLevel =rhs.getRedistributed().size()>0;
551  if(processNextLevel){
552  //restrict defect to coarse level right hand side.
553  typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
554  ++pinfo;
556  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo);
557  }
558  }else{
559  //restrict defect to coarse level right hand side.
560  typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
561  ++pinfo;
563  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
564  }
565 
566  if(processNextLevel){
567  // prepare coarse system
568  ++lhs;
569  ++update;
570  ++matrix;
571  ++level;
572  ++redist;
573 
574  if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
575  // next level is not the globally coarsest one
576  ++smoother;
577  ++aggregates;
578  }
579  // prepare the update on the next level
580  *update=0;
581  }
582  return processNextLevel;
583  }
584 
585  template<class M, class X, class S, class PI, class A>
586  void AMG<M,X,S,PI,A>
587  ::moveToFineLevel(bool processNextLevel)
588  {
589  if(processNextLevel){
590  if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
591  // previous level is not the globally coarsest one
592  --smoother;
593  --aggregates;
594  }
595  --redist;
596  --level;
597  //prolongate and add the correction (update is in coarse left hand side)
598  --matrix;
599 
600  //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
601  --lhs;
602  --pinfo;
603  }
604  if(redist->isSetup()){
605  // Need to redistribute during prolongate
606  lhs.getRedistributed()=0;
608  ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(),
609  *pinfo, *redist);
610  }else{
611  *lhs=0;
613  ::prolongate(*(*aggregates), *update, *lhs,
614  matrices_->getProlongationDampingFactor(), *pinfo);
615  }
616 
617 
618  if(processNextLevel){
619  --update;
620  --rhs;
621  }
622 
623  *update += *lhs;
624  }
625 
626 
627  template<class M, class X, class S, class PI, class A>
628  void AMG<M,X,S,PI,A>
629  ::presmooth()
630  {
631 
632  for(std::size_t i=0; i < preSteps_; ++i){
633  *lhs=0;
634  SmootherApplier<S>::preSmooth(*smoother, *lhs, *rhs);
635  // Accumulate update
636  *update += *lhs;
637 
638  // update defect
639  matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
640  pinfo->project(*rhs);
641  }
642  }
643 
644  template<class M, class X, class S, class PI, class A>
645  void AMG<M,X,S,PI,A>
646  ::postsmooth()
647  {
648 
649  for(std::size_t i=0; i < postSteps_; ++i){
650  // update defect
651  matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
652  *lhs=0;
653  pinfo->project(*rhs);
654  SmootherApplier<S>::postSmooth(*smoother, *lhs, *rhs);
655  // Accumulate update
656  *update += *lhs;
657  }
658  }
659 
660 
661  template<class M, class X, class S, class PI, class A>
663  {
665  }
666 
667  template<class M, class X, class S, class PI, class A>
668  void AMG<M,X,S,PI,A>::mgc(){
669  if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){
670  // Solve directly
672  res.converged=true; // If we do not compute this flag will not get updated
673  if(redist->isSetup()){
674  redist->redistribute(*rhs, rhs.getRedistributed());
675  if(rhs.getRedistributed().size()>0){
676  // We are still participating in the computation
677  pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed());
678  solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res);
679  }
680  redist->redistributeBackward(*update, update.getRedistributed());
681  pinfo->copyOwnerToAll(*update, *update);
682  }else{
683  pinfo->copyOwnerToAll(*rhs, *rhs);
684  solver_->apply(*update, *rhs, res);
685  }
686 
687  if (!res.converged)
688  coarsesolverconverged = false;
689  }else{
690  // presmoothing
691  presmooth();
692 
693 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
694  bool processNextLevel = moveToCoarseLevel();
695 
696  if(processNextLevel){
697  // next level
698  for(std::size_t i=0; i<gamma_; i++)
699  mgc();
700  }
701 
702  moveToFineLevel(processNextLevel);
703 #else
704  *lhs=0;
705 #endif
706 
707  if(matrix == matrices_->matrices().finest()){
708  coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
709  if(!coarsesolverconverged)
710  DUNE_THROW(MathError, "Coarse solver did not converge");
711  }
712  // postsmoothing
713  postsmooth();
714 
715  }
716  }
717 
718  template<class M, class X, class S, class PI, class A>
719  void AMG<M,X,S,PI,A>::additiveMgc(){
720 
721  // restrict residual to all levels
722  typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
723  typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
724  typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
725  typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
726 
727  for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates){
728  ++pinfo;
730  ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
731  }
732 
733  // pinfo is invalid, set to coarsest level
734  //pinfo = matrices_->parallelInformation().coarsest
735  // calculate correction for all levels
736  lhs = lhs_->finest();
737  typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest();
738 
739  for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){
740  // presmoothing
741  *lhs=0;
742  smoother->apply(*lhs, *rhs);
743  }
744 
745  // Coarse level solve
746 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
747  InverseOperatorResult res;
748  pinfo->copyOwnerToAll(*rhs, *rhs);
749  solver_->apply(*lhs, *rhs, res);
750 
751  if(!res.converged)
752  DUNE_THROW(MathError, "Coarse solver did not converge");
753 #else
754  *lhs=0;
755 #endif
756  // Prologate and add up corrections from all levels
757  --pinfo;
758  --aggregates;
759 
760  for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo){
762  ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo);
763  }
764  }
765 
766 
768  template<class M, class X, class S, class PI, class A>
770  {
771  if(buildHierarchy_){
772  if(solver_)
773  delete solver_;
774  if(coarseSmoother_)
776  }
777 
778  // Postprocess all smoothers
779  typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
780  typedef typename Hierarchy<Range,A>::Iterator RIterator;
781  typedef typename Hierarchy<Domain,A>::Iterator DIterator;
782  Iterator coarsest = smoothers_.coarsest();
783  Iterator smoother = smoothers_.finest();
784  DIterator lhs = lhs_->finest();
785  if(smoothers_.levels()>0){
786  if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
787  smoother->post(*lhs);
788  if(smoother!=coarsest)
789  for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
790  smoother->post(*lhs);
791  smoother->post(*lhs);
792  }
793 
794  delete &(*lhs_->finest());
795  delete lhs_;
796  delete &(*update_->finest());
797  delete update_;
798  delete &(*rhs_->finest());
799  delete rhs_;
800  }
801 
802  template<class M, class X, class S, class PI, class A>
803  template<class A1>
804  void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
805  {
806  matrices_->getCoarsestAggregatesOnFinest(cont);
807  }
808 
809  } // end namespace Amg
810 }// end namespace Dune
811 
812 #endif