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