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>
16#include <dune/common/typetraits.hh>
17#include <dune/common/exceptions.hh>
18
19namespace 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;
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,
109
121 template<class C>
125
129 AMG(const AMG& amg);
130
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
171
172 private:
179 template<class C>
180 void createHierarchies(C& criterion, Operator& matrix,
181 const PI& pinfo);
188 struct LevelContext
189 {
198 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
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_;
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>
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>
333 const C& criterion,
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
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>
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 }
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;
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>
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>
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;
604 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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;
613 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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;
660 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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;
667 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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>
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>
742 void AMG<M,X,S,PI,A>::additiveMgc(){
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;
752 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
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
770 InverseOperatorResult res;
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) {
784 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
785 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
786 }
787 }
788
789
791 template<class M, class X, class S, class PI, class A>
793 {
795 // Postprocess all smoothers
796 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
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
Classes for using UMFPack with ISTL matrices.
Classes for using SuperLU with ISTL matrices.
Define base class for scalar product and norm.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Prolongation and restriction for amg.
Provides a classes representing the hierarchies in AMG.
Classes for the generic construction and application of the smoothers.
Col col
Definition matrixmatrix.hh:347
Matrix & mat
Definition matrixmatrix.hh:343
AMG(const AMG &amg)
Copy constructor.
Definition amg.hh:294
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
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition amg.hh:455
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition amg.hh:218
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition amg.hh:222
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition amg.hh:683
X Domain
The domain type.
Definition amg.hh:78
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition amg.hh:202
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition amg.hh:210
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition amg.hh:91
S Smoother
The type of the smoother.
Definition amg.hh:88
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition amg.hh:194
M Operator
The matrix operator type.
Definition amg.hh:64
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition amg.hh:198
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition amg.hh:206
X Range
The range type.
Definition amg.hh:80
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition smoother.hh:408
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition amg.hh:822
std::size_t levels()
Definition amg.hh:535
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition amg.hh:214
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition amg.hh:161
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition amg.hh:75
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition amg.hh:82
void post(Domain &x)
Clean up.
Definition amg.hh:792
std::size_t maxlevels()
Definition amg.hh:540
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition smoother.hh:430
static T * construct(Arguments &args)
Construct an object with the specified arguments.
Definition construction.hh:52
std::size_t level
The level index.
Definition amg.hh:226
AMG(const Operator &fineOperator, const C &criterion, const SmootherArgs &smootherArgs=SmootherArgs(), const ParallelInformation &pinfo=ParallelInformation())
Construct an AMG with an inexact coarse solver based on the smoother.
Definition amg.hh:332
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition amg.hh:547
Smoother SmootherType
Definition amg.hh:190
~AMG()
Definition amg.hh:354
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition amg.hh:73
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition amg.hh:71
@ category
The solver category.
Definition amg.hh:95
Definition basearray.hh:19
Iterator implementation class
Definition basearray.hh:84
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
ConstIterator class for sequential access.
Definition matrix.hh:398
A generic dynamic dense matrix.
Definition matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition matrix.hh:559
T block_type
Export the type representing the components.
Definition matrix.hh:562
an algebraic multigrid method using a Krylov-cycle.
Definition kamg.hh:141
Two grid operator for AMG with Krylov cycle.
Definition kamg.hh:31
Parallel algebraic multigrid based on agglomeration.
Definition amg.hh:56
The hierarchies build by the coarsening process.
Definition hierarchy.hh:317
All parameters for AMG.
Definition parameters.hh:391
Definition pinfo.hh:26
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:26
Statistics about the application of an inverse operator.
Definition solver.hh:32
bool converged
True if convergence criterion has been met.
Definition solver.hh:56
Definition solvertype.hh:14
Definition example.cc:36