Dune Core Modules (2.5.0)

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>
11#include <dune/istl/solvers.hh>
13#include <dune/istl/superlu.hh>
14#include <dune/istl/umfpack.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,
108 const SmootherArgs& smootherArgs, const Parameters& parms);
109
121 template<class C>
122 AMG(const Operator& fineOperator, const C& criterion,
123 const SmootherArgs& smootherArgs=SmootherArgs(),
125
129 AMG(const AMG& amg);
130
131 ~AMG();
132
134 void pre(Domain& x, Range& b);
135
137 void apply(Domain& v, const Range& d);
138
140 void post(Domain& x);
141
146 template<class A1>
147 void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
148
149 std::size_t levels();
150
151 std::size_t maxlevels();
152
162 {
163 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
164 }
165
171
172 private:
179 template<class C>
180 void createHierarchies(C& criterion, Operator& matrix,
181 const PI& pinfo);
188 struct LevelContext
189 {
190 typedef Smoother SmootherType;
206 typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
210 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
226 std::size_t level;
227 };
228
229
234 void mgc(LevelContext& levelContext);
235
236 void additiveMgc();
237
244 void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
245
250 bool moveToCoarseLevel(LevelContext& levelContext);
251
256 void initIteratorsWithFineLevel(LevelContext& levelContext);
257
259 std::shared_ptr<OperatorHierarchy> matrices_;
261 SmootherArgs smootherArgs_;
263 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
265 std::shared_ptr<CoarseSolver> solver_;
267 Hierarchy<Range,A>* rhs_;
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>
314 const SmootherArgs& smootherArgs,
315 const Parameters& parms)
316 : matrices_(&matrices), smootherArgs_(smootherArgs),
317 smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
318 rhs_(), lhs_(), update_(), scalarProduct_(0),
319 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
320 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
321 additive(parms.getAdditive()), coarsesolverconverged(true),
322 coarseSmoother_(), verbosity_(parms.debugLevel())
323 {
324 assert(matrices_->isBuilt());
325
326 // build the necessary smoother hierarchies
327 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
328 }
329
330 template<class M, class X, class S, class PI, class A>
331 template<class C>
333 const C& criterion,
334 const SmootherArgs& smootherArgs,
335 const PI& pinfo)
336 : smootherArgs_(smootherArgs),
337 smoothers_(new Hierarchy<Smoother,A>), solver_(),
338 rhs_(), lhs_(), update_(), scalarProduct_(),
339 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
340 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
341 additive(criterion.getAdditive()), coarsesolverconverged(true),
342 coarseSmoother_(), verbosity_(criterion.debugLevel())
343 {
344 static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
345 "Matrix and Solver must match in terms of category!");
346 // TODO: reestablish compile time checks.
347 //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
348 // "Matrix and Solver must match in terms of category!");
349 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
350 }
351
352
353 template<class M, class X, class S, class PI, class A>
355 {
356 if(buildHierarchy_) {
357 if(solver_)
358 solver_.reset();
359 if(coarseSmoother_)
360 coarseSmoother_.reset();
361 }
362 if(lhs_)
363 delete lhs_;
364 lhs_=nullptr;
365 if(update_)
366 delete update_;
367 update_=nullptr;
368 if(rhs_)
369 delete rhs_;
370 rhs_=nullptr;
371 }
372
373 template <class Matrix,
374 class Vector>
375 struct DirectSolverSelector
376 {
377 typedef typename Matrix :: field_type field_type;
378 enum SolverType { umfpack, superlu, none };
379
380 static constexpr SolverType solver =
381#if HAVE_SUITESPARSE_UMFPACK
382 UMFPackMethodChooser< field_type > :: valid ? umfpack : none ;
383#elif HAVE_SUPERLU
384 superlu ;
385#else
386 none;
387#endif
388
389 template <class M, SolverType>
390 struct Solver
391 {
392 typedef InverseOperator<Vector,Vector> type;
393 static type* create(const M& mat, bool verbose, bool reusevector )
394 {
395 DUNE_THROW(NotImplemented,"DirectSolver not selected");
396 return nullptr;
397 }
398 static std::string name () { return "None"; }
399 };
400#if HAVE_SUITESPARSE_UMFPACK
401 template <class M>
402 struct Solver< M, umfpack >
403 {
404 typedef UMFPack< M > type;
405 static type* create(const M& mat, bool verbose, bool reusevector )
406 {
407 return new type(mat, verbose, reusevector );
408 }
409 static std::string name () { return "UMFPack"; }
410 };
411#endif
412#if HAVE_SUPERLU
413 template <class M>
414 struct Solver< M, superlu >
415 {
416 typedef SuperLU< M > type;
417 static type* create(const M& mat, bool verbose, bool reusevector )
418 {
419 return new type(mat, verbose, reusevector );
420 }
421 static std::string name () { return "SuperLU"; }
422 };
423#endif
424
425 // define direct solver type to be used
426 typedef Solver< Matrix, solver > SelectedSolver ;
427 typedef typename SelectedSolver :: type DirectSolver;
428 static constexpr bool isDirectSolver = solver != none;
429 static std::string name() { return SelectedSolver :: name (); }
430 static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
431 {
432 return SelectedSolver :: create( mat, verbose, reusevector );
433 }
434 };
435
436 template<class M, class X, class S, class PI, class A>
437 template<class C>
438 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
439 const PI& pinfo)
440 {
441 Timer watch;
442 matrices_.reset(new OperatorHierarchy(matrix, pinfo));
443
444 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
445
446 // build the necessary smoother hierarchies
447 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
448
449 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels())
450 {
451 // We have the carsest level. Create the coarse Solver
452 SmootherArgs sargs(smootherArgs_);
453 sargs.iterations = 1;
454
456 cargs.setArgs(sargs);
457 if(matrices_->redistributeInformation().back().isSetup()) {
458 // Solve on the redistributed partitioning
459 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
460 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
461 }else{
462 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
463 cargs.setComm(*matrices_->parallelInformation().coarsest());
464 }
465
466 coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
467 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
468
469 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
470
471 // Use superlu if we are purely sequential or with only one processor on the coarsest level.
472 if( SolverSelector::isDirectSolver &&
473 (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
474 || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
475 || (matrices_->parallelInformation().coarsest().isRedistributed()
476 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
477 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
478 )
479 { // redistribute and 1 proc
480 if(matrices_->parallelInformation().coarsest().isRedistributed())
481 {
482 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
483 {
484 // We are still participating on this level
485 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
486 }
487 else
488 solver_.reset();
489 }
490 else
491 {
492 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
493 }
494 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
495 std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
496 }
497 else
498 {
499 if(matrices_->parallelInformation().coarsest().isRedistributed())
500 {
501 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
502 // We are still participating on this level
503 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
504 *scalarProduct_,
505 *coarseSmoother_, 1E-2, 1000, 0));
506 else
507 solver_.reset();
508 }else
509 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
510 *scalarProduct_,
511 *coarseSmoother_, 1E-2, 1000, 0));
512 }
513 }
514
515 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
516 std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
517 <<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
518 }
519
520
521 template<class M, class X, class S, class PI, class A>
523 {
524 // Detect Matrix rows where all offdiagonal entries are
525 // zero and set x such that A_dd*x_d=b_d
526 // Thus users can be more careless when setting up their linear
527 // systems.
528 typedef typename M::matrix_type Matrix;
529 typedef typename Matrix::ConstRowIterator RowIter;
530 typedef typename Matrix::ConstColIterator ColIter;
531 typedef typename Matrix::block_type Block;
532 Block zero;
533 zero=typename Matrix::field_type();
534
535 const Matrix& mat=matrices_->matrices().finest()->getmat();
536 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
537 bool isDirichlet = true;
538 bool hasDiagonal = false;
539 Block diagonal;
540 for(ColIter col=row->begin(); col!=row->end(); ++col) {
541 if(row.index()==col.index()) {
542 diagonal = *col;
543 hasDiagonal = false;
544 }else{
545 if(*col!=zero)
546 isDirichlet = false;
547 }
548 }
549 if(isDirichlet && hasDiagonal)
550 diagonal.solve(x[row.index()], b[row.index()]);
551 }
552
553 if(smoothers_->levels()>0)
554 smoothers_->finest()->pre(x,b);
555 else
556 // No smoother to make x consistent! Do it by hand
557 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
558 Range* copy = new Range(b);
559 if(rhs_)
560 delete rhs_;
561 rhs_ = new Hierarchy<Range,A>(copy);
562 Domain* dcopy = new Domain(x);
563 if(lhs_)
564 delete lhs_;
565 lhs_ = new Hierarchy<Domain,A>(dcopy);
566 dcopy = new Domain(x);
567 if(update_)
568 delete update_;
569 update_ = new Hierarchy<Domain,A>(dcopy);
570 matrices_->coarsenVector(*rhs_);
571 matrices_->coarsenVector(*lhs_);
572 matrices_->coarsenVector(*update_);
573
574 // Preprocess all smoothers
575 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
576 typedef typename Hierarchy<Range,A>::Iterator RIterator;
577 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
578 Iterator coarsest = smoothers_->coarsest();
579 Iterator smoother = smoothers_->finest();
580 RIterator rhs = rhs_->finest();
581 DIterator lhs = lhs_->finest();
582 if(smoothers_->levels()>0) {
583
584 assert(lhs_->levels()==rhs_->levels());
585 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
586 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
587
588 if(smoother!=coarsest)
589 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
590 smoother->pre(*lhs,*rhs);
591 smoother->pre(*lhs,*rhs);
592 }
593
594
595 // The preconditioner might change x and b. So we have to
596 // copy the changes to the original vectors.
597 x = *lhs_->finest();
598 b = *rhs_->finest();
599
600 }
601 template<class M, class X, class S, class PI, class A>
602 std::size_t AMG<M,X,S,PI,A>::levels()
603 {
604 return matrices_->levels();
605 }
606 template<class M, class X, class S, class PI, class A>
607 std::size_t AMG<M,X,S,PI,A>::maxlevels()
608 {
609 return matrices_->maxlevels();
610 }
611
613 template<class M, class X, class S, class PI, class A>
615 {
616 LevelContext levelContext;
617
618 if(additive) {
619 *(rhs_->finest())=d;
620 additiveMgc();
621 v=*lhs_->finest();
622 }else{
623 // Init all iterators for the current level
624 initIteratorsWithFineLevel(levelContext);
625
626
627 *levelContext.lhs = v;
628 *levelContext.rhs = d;
629 *levelContext.update=0;
630 levelContext.level=0;
631
632 mgc(levelContext);
633
634 if(postSteps_==0||matrices_->maxlevels()==1)
635 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
636
637 v=*levelContext.update;
638 }
639
640 }
641
642 template<class M, class X, class S, class PI, class A>
643 void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
644 {
645 levelContext.smoother = smoothers_->finest();
646 levelContext.matrix = matrices_->matrices().finest();
647 levelContext.pinfo = matrices_->parallelInformation().finest();
648 levelContext.redist =
649 matrices_->redistributeInformation().begin();
650 levelContext.aggregates = matrices_->aggregatesMaps().begin();
651 levelContext.lhs = lhs_->finest();
652 levelContext.update = update_->finest();
653 levelContext.rhs = rhs_->finest();
654 }
655
656 template<class M, class X, class S, class PI, class A>
657 bool AMG<M,X,S,PI,A>
658 ::moveToCoarseLevel(LevelContext& levelContext)
659 {
660
661 bool processNextLevel=true;
662
663 if(levelContext.redist->isSetup()) {
664 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
665 levelContext.rhs.getRedistributed());
666 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
667 if(processNextLevel) {
668 //restrict defect to coarse level right hand side.
669 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
670 ++levelContext.pinfo;
671 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
672 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
673 static_cast<const Range&>(fineRhs.getRedistributed()),
674 *levelContext.pinfo);
675 }
676 }else{
677 //restrict defect to coarse level right hand side.
678 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
679 ++levelContext.pinfo;
680 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
681 ::restrictVector(*(*levelContext.aggregates),
682 *levelContext.rhs, static_cast<const Range&>(*fineRhs),
683 *levelContext.pinfo);
684 }
685
686 if(processNextLevel) {
687 // prepare coarse system
688 ++levelContext.lhs;
689 ++levelContext.update;
690 ++levelContext.matrix;
691 ++levelContext.level;
692 ++levelContext.redist;
693
694 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
695 // next level is not the globally coarsest one
696 ++levelContext.smoother;
697 ++levelContext.aggregates;
698 }
699 // prepare the update on the next level
700 *levelContext.update=0;
701 }
702 return processNextLevel;
703 }
704
705 template<class M, class X, class S, class PI, class A>
706 void AMG<M,X,S,PI,A>
707 ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
708 {
709 if(processNextLevel) {
710 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
711 // previous level is not the globally coarsest one
712 --levelContext.smoother;
713 --levelContext.aggregates;
714 }
715 --levelContext.redist;
716 --levelContext.level;
717 //prolongate and add the correction (update is in coarse left hand side)
718 --levelContext.matrix;
719
720 //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
721 --levelContext.lhs;
722 --levelContext.pinfo;
723 }
724 if(levelContext.redist->isSetup()) {
725 // Need to redistribute during prolongateVector
726 levelContext.lhs.getRedistributed()=0;
727 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
728 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
729 levelContext.lhs.getRedistributed(),
730 matrices_->getProlongationDampingFactor(),
731 *levelContext.pinfo, *levelContext.redist);
732 }else{
733 *levelContext.lhs=0;
734 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
735 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
736 matrices_->getProlongationDampingFactor(),
737 *levelContext.pinfo);
738 }
739
740
741 if(processNextLevel) {
742 --levelContext.update;
743 --levelContext.rhs;
744 }
745
746 *levelContext.update += *levelContext.lhs;
747 }
748
749 template<class M, class X, class S, class PI, class A>
751 {
752 return IsDirectSolver< CoarseSolver>::value;
753 }
754
755 template<class M, class X, class S, class PI, class A>
756 void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
757 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
758 // Solve directly
760 res.converged=true; // If we do not compute this flag will not get updated
761 if(levelContext.redist->isSetup()) {
762 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
763 if(levelContext.rhs.getRedistributed().size()>0) {
764 // We are still participating in the computation
765 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
766 levelContext.rhs.getRedistributed());
767 solver_->apply(levelContext.update.getRedistributed(),
768 levelContext.rhs.getRedistributed(), res);
769 }
770 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
771 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
772 }else{
773 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
774 solver_->apply(*levelContext.update, *levelContext.rhs, res);
775 }
776
777 if (!res.converged)
778 coarsesolverconverged = false;
779 }else{
780 // presmoothing
781 presmooth(levelContext, preSteps_);
782
783#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
784 bool processNextLevel = moveToCoarseLevel(levelContext);
785
786 if(processNextLevel) {
787 // next level
788 for(std::size_t i=0; i<gamma_; i++)
789 mgc(levelContext);
790 }
791
792 moveToFineLevel(levelContext, processNextLevel);
793#else
794 *lhs=0;
795#endif
796
797 if(levelContext.matrix == matrices_->matrices().finest()) {
798 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
799 if(!coarsesolverconverged)
800 DUNE_THROW(MathError, "Coarse solver did not converge");
801 }
802 // postsmoothing
803 postsmooth(levelContext, postSteps_);
804
805 }
806 }
807
808 template<class M, class X, class S, class PI, class A>
809 void AMG<M,X,S,PI,A>::additiveMgc(){
810
811 // restrict residual to all levels
812 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
813 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
814 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
815 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
816
817 for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
818 ++pinfo;
819 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
820 ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
821 }
822
823 // pinfo is invalid, set to coarsest level
824 //pinfo = matrices_->parallelInformation().coarsest
825 // calculate correction for all levels
826 lhs = lhs_->finest();
827 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
828
829 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
830 // presmoothing
831 *lhs=0;
832 smoother->apply(*lhs, *rhs);
833 }
834
835 // Coarse level solve
836#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
837 InverseOperatorResult res;
838 pinfo->copyOwnerToAll(*rhs, *rhs);
839 solver_->apply(*lhs, *rhs, res);
840
841 if(!res.converged)
842 DUNE_THROW(MathError, "Coarse solver did not converge");
843#else
844 *lhs=0;
845#endif
846 // Prologate and add up corrections from all levels
847 --pinfo;
848 --aggregates;
849
850 for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
851 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
852 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
853 }
854 }
855
856
858 template<class M, class X, class S, class PI, class A>
860 {
862 // Postprocess all smoothers
863 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
864 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
865 Iterator coarsest = smoothers_->coarsest();
866 Iterator smoother = smoothers_->finest();
867 DIterator lhs = lhs_->finest();
868 if(smoothers_->levels()>0) {
869 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
870 smoother->post(*lhs);
871 if(smoother!=coarsest)
872 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
873 smoother->post(*lhs);
874 smoother->post(*lhs);
875 }
876 //delete &(*lhs_->finest());
877 delete lhs_;
878 lhs_=nullptr;
879 //delete &(*update_->finest());
880 delete update_;
881 update_=nullptr;
882 //delete &(*rhs_->finest());
883 delete rhs_;
884 rhs_=nullptr;
885 }
886
887 template<class M, class X, class S, class PI, class A>
888 template<class A1>
889 void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
890 {
891 matrices_->getCoarsestAggregatesOnFinest(cont);
892 }
893
894 } // end namespace Amg
895} // end namespace Dune
896
897#endif
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:56
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:257
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:260
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:141
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:31
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:317
All parameters for AMG.
Definition: parameters.hh:391
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
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
T block_type
Export the type representing the components.
Definition: matrix.hh:562
The negation of a set. An item is contained in the set if and only if it is not contained in the nega...
Definition: enumset.hh:95
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:26
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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:522
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:750
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:889
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:214
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:161
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:1385
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:859
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:614
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
Provides a classes representing the hierarchies in AMG.
Dune namespace.
Definition: alignment.hh:11
Define base class for scalar product and norm.
Classes for the generic construction and application of the smoothers.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
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
Choose the approriate scalar product for a solver category.
Definition: scalarproducts.hh:77
Classes for using SuperLU with ISTL matrices.
Prolongation and restriction for amg.
Traits for type conversions and type information.
Classes for using UMFPack with ISTL matrices.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)