Dune Core Modules (2.4.1)

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
114 AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
115 const SmootherArgs& smootherArgs, std::size_t gamma,
116 std::size_t preSmoothingSteps,
117 std::size_t postSmoothingSteps,
118 bool additive=false) DUNE_DEPRECATED;
119
129 AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
130 const SmootherArgs& smootherArgs, const Parameters& parms);
131
151 template<class C>
152 AMG(const Operator& fineOperator, const C& criterion,
153 const SmootherArgs& smootherArgs, std::size_t gamma,
154 std::size_t preSmoothingSteps,
155 std::size_t postSmoothingSteps,
156 bool additive=false,
158
170 template<class C>
171 AMG(const Operator& fineOperator, const C& criterion,
172 const SmootherArgs& smootherArgs=SmootherArgs(),
174
178 AMG(const AMG& amg);
179
180 ~AMG();
181
183 void pre(Domain& x, Range& b);
184
186 void apply(Domain& v, const Range& d);
187
189 void post(Domain& x);
190
195 template<class A1>
196 void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
197
198 std::size_t levels();
199
200 std::size_t maxlevels();
201
211 {
212 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
213 }
214
220
221 private:
228 template<class C>
229 void createHierarchies(C& criterion, Operator& matrix,
230 const PI& pinfo);
237 struct LevelContext
238 {
239 typedef Smoother SmootherType;
255 typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
259 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
275 std::size_t level;
276 };
277
278
283 void mgc(LevelContext& levelContext);
284
285 void additiveMgc();
286
293 void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
294
299 bool moveToCoarseLevel(LevelContext& levelContext);
300
305 void initIteratorsWithFineLevel(LevelContext& levelContext);
306
308 std::shared_ptr<OperatorHierarchy> matrices_;
310 SmootherArgs smootherArgs_;
312 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
314 std::shared_ptr<CoarseSolver> solver_;
316 Hierarchy<Range,A>* rhs_;
320 Hierarchy<Domain,A>* update_;
324 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
325 typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
327 ScalarProductPointer scalarProduct_;
329 std::size_t gamma_;
331 std::size_t preSteps_;
333 std::size_t postSteps_;
334 bool buildHierarchy_;
335 bool additive;
336 bool coarsesolverconverged;
337 std::shared_ptr<Smoother> coarseSmoother_;
339 std::size_t verbosity_;
340 };
341
342 template<class M, class X, class S, class PI, class A>
343 inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
344 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
345 smoothers_(amg.smoothers_), solver_(amg.solver_),
346 rhs_(), lhs_(), update_(),
347 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
348 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
349 buildHierarchy_(amg.buildHierarchy_),
350 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
351 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
352 {
353 if(amg.rhs_)
354 rhs_=new Hierarchy<Range,A>(*amg.rhs_);
355 if(amg.lhs_)
356 lhs_=new Hierarchy<Domain,A>(*amg.lhs_);
357 if(amg.update_)
358 update_=new Hierarchy<Domain,A>(*amg.update_);
359 }
360
361 template<class M, class X, class S, class PI, class A>
363 const SmootherArgs& smootherArgs,
364 std::size_t gamma, std::size_t preSmoothingSteps,
365 std::size_t postSmoothingSteps, bool additive_)
366 : matrices_(&matrices), smootherArgs_(smootherArgs),
367 smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
368 rhs_(), lhs_(), update_(), scalarProduct_(0),
369 gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
370 additive(additive_), coarsesolverconverged(true),
371 coarseSmoother_(), verbosity_(2)
372 {
373 assert(matrices_->isBuilt());
374
375 // build the necessary smoother hierarchies
376 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
377 }
378
379 template<class M, class X, class S, class PI, class A>
381 const SmootherArgs& smootherArgs,
382 const Parameters& parms)
383 : matrices_(&matrices), smootherArgs_(smootherArgs),
384 smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
385 rhs_(), lhs_(), update_(), scalarProduct_(0),
386 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
387 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
388 additive(parms.getAdditive()), coarsesolverconverged(true),
389 coarseSmoother_(), verbosity_(parms.debugLevel())
390 {
391 assert(matrices_->isBuilt());
392
393 // build the necessary smoother hierarchies
394 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
395 }
396
397 template<class M, class X, class S, class PI, class A>
398 template<class C>
400 const C& criterion,
401 const SmootherArgs& smootherArgs,
402 std::size_t gamma, std::size_t preSmoothingSteps,
403 std::size_t postSmoothingSteps,
404 bool additive_,
405 const PI& pinfo)
406 : smootherArgs_(smootherArgs),
407 smoothers_(new Hierarchy<Smoother,A>), solver_(),
408 rhs_(), lhs_(), update_(), scalarProduct_(), gamma_(gamma),
409 preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
410 additive(additive_), coarsesolverconverged(true),
411 coarseSmoother_(), verbosity_(criterion.debugLevel())
412 {
413 static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
414 "Matrix and Solver must match in terms of category!");
415 // TODO: reestablish compile time checks.
416 //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
417 // "Matrix and Solver must match in terms of category!");
418 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
419 }
420
421 template<class M, class X, class S, class PI, class A>
422 template<class C>
424 const C& criterion,
425 const SmootherArgs& smootherArgs,
426 const PI& pinfo)
427 : smootherArgs_(smootherArgs),
428 smoothers_(new Hierarchy<Smoother,A>), solver_(),
429 rhs_(), lhs_(), update_(), scalarProduct_(),
430 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
431 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
432 additive(criterion.getAdditive()), coarsesolverconverged(true),
433 coarseSmoother_(), verbosity_(criterion.debugLevel())
434 {
435 static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
436 "Matrix and Solver must match in terms of category!");
437 // TODO: reestablish compile time checks.
438 //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
439 // "Matrix and Solver must match in terms of category!");
440 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
441 }
442
443
444 template<class M, class X, class S, class PI, class A>
446 {
447 if(buildHierarchy_) {
448 if(solver_)
449 solver_.reset();
450 if(coarseSmoother_)
451 coarseSmoother_.reset();
452 }
453 if(lhs_)
454 delete lhs_;
455 lhs_=nullptr;
456 if(update_)
457 delete update_;
458 update_=nullptr;
459 if(rhs_)
460 delete rhs_;
461 rhs_=nullptr;
462 }
463
464 template<class M, class X, class S, class PI, class A>
465 template<class C>
466 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
467 const PI& pinfo)
468 {
469 Timer watch;
470 matrices_.reset(new OperatorHierarchy(matrix, pinfo));
471
472 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
473
474 // build the necessary smoother hierarchies
475 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
476
477 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
478 // We have the carsest level. Create the coarse Solver
479 SmootherArgs sargs(smootherArgs_);
480 sargs.iterations = 1;
481
483 cargs.setArgs(sargs);
484 if(matrices_->redistributeInformation().back().isSetup()) {
485 // Solve on the redistributed partitioning
486 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
487 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
488 }else{
489 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
490 cargs.setComm(*matrices_->parallelInformation().coarsest());
491 }
492
493 coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
494 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
495
496#if HAVE_SUPERLU || HAVE_UMFPACK
497#if HAVE_UMFPACK
498#define DIRECTSOLVER UMFPack
499#else
500#define DIRECTSOLVER SuperLU
501#endif
502 // Use superlu if we are purely sequential or with only one processor on the coarsest level.
503 if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode
504 || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
505 || (matrices_->parallelInformation().coarsest().isRedistributed()
506 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
507 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) { // redistribute and 1 proc
508 if(matrices_->parallelInformation().coarsest().isRedistributed())
509 {
510 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
511 // We are still participating on this level
512 solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
513 else
514 solver_.reset();
515 }else
516 solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(), false, false));
517 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
518 std::cout<< "Using a direct coarse solver (" << static_cast< DIRECTSOLVER<typename M::matrix_type>* >(solver_.get())->name() << ")" << std::endl;
519 }else
520#undef DIRECTSOLVER
521#endif
522 {
523 if(matrices_->parallelInformation().coarsest().isRedistributed())
524 {
525 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
526 // We are still participating on this level
527 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
528 *scalarProduct_,
529 *coarseSmoother_, 1E-2, 1000, 0));
530 else
531 solver_.reset();
532 }else
533 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
534 *scalarProduct_,
535 *coarseSmoother_, 1E-2, 1000, 0));
536 }
537 }
538
539 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
540 std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
541 <<"(inclusive coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
542 }
543
544
545 template<class M, class X, class S, class PI, class A>
547 {
548 // Detect Matrix rows where all offdiagonal entries are
549 // zero and set x such that A_dd*x_d=b_d
550 // Thus users can be more careless when setting up their linear
551 // systems.
552 typedef typename M::matrix_type Matrix;
553 typedef typename Matrix::ConstRowIterator RowIter;
554 typedef typename Matrix::ConstColIterator ColIter;
555 typedef typename Matrix::block_type Block;
556 Block zero;
557 zero=typename Matrix::field_type();
558
559 const Matrix& mat=matrices_->matrices().finest()->getmat();
560 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
561 bool isDirichlet = true;
562 bool hasDiagonal = false;
563 Block diagonal;
564 for(ColIter col=row->begin(); col!=row->end(); ++col) {
565 if(row.index()==col.index()) {
566 diagonal = *col;
567 hasDiagonal = false;
568 }else{
569 if(*col!=zero)
570 isDirichlet = false;
571 }
572 }
573 if(isDirichlet && hasDiagonal)
574 diagonal.solve(x[row.index()], b[row.index()]);
575 }
576
577 if(smoothers_->levels()>0)
578 smoothers_->finest()->pre(x,b);
579 else
580 // No smoother to make x consistent! Do it by hand
581 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
582 Range* copy = new Range(b);
583 if(rhs_)
584 delete rhs_;
585 rhs_ = new Hierarchy<Range,A>(copy);
586 Domain* dcopy = new Domain(x);
587 if(lhs_)
588 delete lhs_;
589 lhs_ = new Hierarchy<Domain,A>(dcopy);
590 dcopy = new Domain(x);
591 if(update_)
592 delete update_;
593 update_ = new Hierarchy<Domain,A>(dcopy);
594 matrices_->coarsenVector(*rhs_);
595 matrices_->coarsenVector(*lhs_);
596 matrices_->coarsenVector(*update_);
597
598 // Preprocess all smoothers
599 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
600 typedef typename Hierarchy<Range,A>::Iterator RIterator;
601 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
602 Iterator coarsest = smoothers_->coarsest();
603 Iterator smoother = smoothers_->finest();
604 RIterator rhs = rhs_->finest();
605 DIterator lhs = lhs_->finest();
606 if(smoothers_->levels()>0) {
607
608 assert(lhs_->levels()==rhs_->levels());
609 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
610 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
611
612 if(smoother!=coarsest)
613 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
614 smoother->pre(*lhs,*rhs);
615 smoother->pre(*lhs,*rhs);
616 }
617
618
619 // The preconditioner might change x and b. So we have to
620 // copy the changes to the original vectors.
621 x = *lhs_->finest();
622 b = *rhs_->finest();
623
624 }
625 template<class M, class X, class S, class PI, class A>
626 std::size_t AMG<M,X,S,PI,A>::levels()
627 {
628 return matrices_->levels();
629 }
630 template<class M, class X, class S, class PI, class A>
631 std::size_t AMG<M,X,S,PI,A>::maxlevels()
632 {
633 return matrices_->maxlevels();
634 }
635
637 template<class M, class X, class S, class PI, class A>
639 {
640 LevelContext levelContext;
641
642 if(additive) {
643 *(rhs_->finest())=d;
644 additiveMgc();
645 v=*lhs_->finest();
646 }else{
647 // Init all iterators for the current level
648 initIteratorsWithFineLevel(levelContext);
649
650
651 *levelContext.lhs = v;
652 *levelContext.rhs = d;
653 *levelContext.update=0;
654 levelContext.level=0;
655
656 mgc(levelContext);
657
658 if(postSteps_==0||matrices_->maxlevels()==1)
659 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
660
661 v=*levelContext.update;
662 }
663
664 }
665
666 template<class M, class X, class S, class PI, class A>
667 void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
668 {
669 levelContext.smoother = smoothers_->finest();
670 levelContext.matrix = matrices_->matrices().finest();
671 levelContext.pinfo = matrices_->parallelInformation().finest();
672 levelContext.redist =
673 matrices_->redistributeInformation().begin();
674 levelContext.aggregates = matrices_->aggregatesMaps().begin();
675 levelContext.lhs = lhs_->finest();
676 levelContext.update = update_->finest();
677 levelContext.rhs = rhs_->finest();
678 }
679
680 template<class M, class X, class S, class PI, class A>
681 bool AMG<M,X,S,PI,A>
682 ::moveToCoarseLevel(LevelContext& levelContext)
683 {
684
685 bool processNextLevel=true;
686
687 if(levelContext.redist->isSetup()) {
688 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
689 levelContext.rhs.getRedistributed());
690 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
691 if(processNextLevel) {
692 //restrict defect to coarse level right hand side.
693 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
694 ++levelContext.pinfo;
695 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
696 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
697 static_cast<const Range&>(fineRhs.getRedistributed()),
698 *levelContext.pinfo);
699 }
700 }else{
701 //restrict defect to coarse level right hand side.
702 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
703 ++levelContext.pinfo;
704 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
705 ::restrictVector(*(*levelContext.aggregates),
706 *levelContext.rhs, static_cast<const Range&>(*fineRhs),
707 *levelContext.pinfo);
708 }
709
710 if(processNextLevel) {
711 // prepare coarse system
712 ++levelContext.lhs;
713 ++levelContext.update;
714 ++levelContext.matrix;
715 ++levelContext.level;
716 ++levelContext.redist;
717
718 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
719 // next level is not the globally coarsest one
720 ++levelContext.smoother;
721 ++levelContext.aggregates;
722 }
723 // prepare the update on the next level
724 *levelContext.update=0;
725 }
726 return processNextLevel;
727 }
728
729 template<class M, class X, class S, class PI, class A>
730 void AMG<M,X,S,PI,A>
731 ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
732 {
733 if(processNextLevel) {
734 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
735 // previous level is not the globally coarsest one
736 --levelContext.smoother;
737 --levelContext.aggregates;
738 }
739 --levelContext.redist;
740 --levelContext.level;
741 //prolongate and add the correction (update is in coarse left hand side)
742 --levelContext.matrix;
743
744 //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
745 --levelContext.lhs;
746 --levelContext.pinfo;
747 }
748 if(levelContext.redist->isSetup()) {
749 // Need to redistribute during prolongateVector
750 levelContext.lhs.getRedistributed()=0;
751 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
752 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
753 levelContext.lhs.getRedistributed(),
754 matrices_->getProlongationDampingFactor(),
755 *levelContext.pinfo, *levelContext.redist);
756 }else{
757 *levelContext.lhs=0;
758 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
759 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
760 matrices_->getProlongationDampingFactor(),
761 *levelContext.pinfo);
762 }
763
764
765 if(processNextLevel) {
766 --levelContext.update;
767 --levelContext.rhs;
768 }
769
770 *levelContext.update += *levelContext.lhs;
771 }
772
773 template<class M, class X, class S, class PI, class A>
775 {
776 return IsDirectSolver< CoarseSolver>::value;
777 }
778
779 template<class M, class X, class S, class PI, class A>
780 void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
781 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
782 // Solve directly
784 res.converged=true; // If we do not compute this flag will not get updated
785 if(levelContext.redist->isSetup()) {
786 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
787 if(levelContext.rhs.getRedistributed().size()>0) {
788 // We are still participating in the computation
789 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
790 levelContext.rhs.getRedistributed());
791 solver_->apply(levelContext.update.getRedistributed(),
792 levelContext.rhs.getRedistributed(), res);
793 }
794 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
795 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
796 }else{
797 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
798 solver_->apply(*levelContext.update, *levelContext.rhs, res);
799 }
800
801 if (!res.converged)
802 coarsesolverconverged = false;
803 }else{
804 // presmoothing
805 presmooth(levelContext, preSteps_);
806
807#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
808 bool processNextLevel = moveToCoarseLevel(levelContext);
809
810 if(processNextLevel) {
811 // next level
812 for(std::size_t i=0; i<gamma_; i++)
813 mgc(levelContext);
814 }
815
816 moveToFineLevel(levelContext, processNextLevel);
817#else
818 *lhs=0;
819#endif
820
821 if(levelContext.matrix == matrices_->matrices().finest()) {
822 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
823 if(!coarsesolverconverged)
824 DUNE_THROW(MathError, "Coarse solver did not converge");
825 }
826 // postsmoothing
827 postsmooth(levelContext, postSteps_);
828
829 }
830 }
831
832 template<class M, class X, class S, class PI, class A>
833 void AMG<M,X,S,PI,A>::additiveMgc(){
834
835 // restrict residual to all levels
836 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
837 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
838 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
839 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
840
841 for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
842 ++pinfo;
843 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
844 ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
845 }
846
847 // pinfo is invalid, set to coarsest level
848 //pinfo = matrices_->parallelInformation().coarsest
849 // calculate correction for all levels
850 lhs = lhs_->finest();
851 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
852
853 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
854 // presmoothing
855 *lhs=0;
856 smoother->apply(*lhs, *rhs);
857 }
858
859 // Coarse level solve
860#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
861 InverseOperatorResult res;
862 pinfo->copyOwnerToAll(*rhs, *rhs);
863 solver_->apply(*lhs, *rhs, res);
864
865 if(!res.converged)
866 DUNE_THROW(MathError, "Coarse solver did not converge");
867#else
868 *lhs=0;
869#endif
870 // Prologate and add up corrections from all levels
871 --pinfo;
872 --aggregates;
873
874 for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
875 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
876 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
877 }
878 }
879
880
882 template<class M, class X, class S, class PI, class A>
884 {
886 // Postprocess all smoothers
887 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
888 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
889 Iterator coarsest = smoothers_->coarsest();
890 Iterator smoother = smoothers_->finest();
891 DIterator lhs = lhs_->finest();
892 if(smoothers_->levels()>0) {
893 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
894 smoother->post(*lhs);
895 if(smoother!=coarsest)
896 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
897 smoother->post(*lhs);
898 smoother->post(*lhs);
899 }
900 //delete &(*lhs_->finest());
901 delete lhs_;
902 lhs_=nullptr;
903 //delete &(*update_->finest());
904 delete update_;
905 update_=nullptr;
906 //delete &(*rhs_->finest());
907 delete rhs_;
908 rhs_=nullptr;
909 }
910
911 template<class M, class X, class S, class PI, class A>
912 template<class A1>
913 void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
914 {
915 matrices_->getCoarsestAggregatesOnFinest(cont);
916 }
917
918 } // end namespace Amg
919} // end namespace Dune
920
921#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
A generic dynamic dense matrix.
Definition: matrix.hh:25
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:29
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:53
T block_type
Export the type representing the components.
Definition: matrix.hh:32
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
ConstIterator class for sequential access.
Definition: vbvector.hh:647
A few common exception classes.
#define DUNE_DEPRECATED
Mark some entity as deprecated.
Definition: deprecated.hh:84
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:546
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:267
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:271
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:774
AMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, bool additive=false) DUNE_DEPRECATED
Construct a new amg with a specific coarse solver.
Definition: amg.hh:362
X Domain
The domain type.
Definition: amg.hh:78
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:251
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:259
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:243
M Operator
The matrix operator type.
Definition: amg.hh:64
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:247
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:255
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:913
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:263
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:210
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:883
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:275
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:638
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:10
STL namespace.
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 intentional 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 13, 23:29, 2024)