dune-istl 2.8.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>
7#include <sstream>
8#include <dune/common/exceptions.hh>
12#include <dune/istl/solvers.hh>
14#include <dune/istl/superlu.hh>
15#include <dune/istl/umfpack.hh>
17#include <dune/common/typetraits.hh>
18#include <dune/common/exceptions.hh>
19#include <dune/common/scalarvectorview.hh>
20#include <dune/common/scalarmatrixview.hh>
21#include <dune/common/parametertree.hh>
22
23namespace Dune
24{
25 namespace Amg
26 {
44 template<class M, class X, class S, class P, class K, class A>
45 class KAMG;
46
47 template<class T>
48 class KAmgTwoGrid;
49
60 template<class M, class X, class S, class PI=SequentialInformation,
61 class A=std::allocator<X> >
62 class AMG : public Preconditioner<X,X>
63 {
64 template<class M1, class X1, class S1, class P1, class K1, class A1>
65 friend class KAMG;
66
67 friend class KAmgTwoGrid<AMG>;
68
69 public:
71 typedef M Operator;
83
85 typedef X Domain;
87 typedef X Range;
95 typedef S Smoother;
96
99
109 AMG(OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
110 const SmootherArgs& smootherArgs, const Parameters& parms);
111
123 template<class C>
124 AMG(const Operator& fineOperator, const C& criterion,
125 const SmootherArgs& smootherArgs=SmootherArgs(),
127
178 AMG(std::shared_ptr<const Operator> fineOperator, const ParameterTree& configuration, const ParallelInformation& pinfo=ParallelInformation());
179
183 AMG(const AMG& amg);
184
186 void pre(Domain& x, Range& b);
187
189 void apply(Domain& v, const Range& d);
190
193 {
194 return category_;
195 }
196
198 void post(Domain& x);
199
204 template<class A1>
205 void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
206
207 std::size_t levels();
208
209 std::size_t maxlevels();
210
220 {
221 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
222 }
223
229
230 private:
231 /*
232 * @brief Helper function to create hierarchies with parameter tree.
233 *
234 * Will create the coarsen criterion with the norm and create the
235 * Hierarchies
236 * \tparam Norm Type of the norm to use.
237 */
238 template<class Norm>
239 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
240 const PI& pinfo, const Norm&,
241 const ParameterTree& configuration,
242 std::true_type compiles = std::true_type());
243 template<class Norm>
244 void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
245 const PI& pinfo, const Norm&,
246 const ParameterTree& configuration,
247 std::false_type);
252 template<class C>
253 void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
254 const PI& pinfo, const ParameterTree& configuration);
261 template<class C>
262 void createHierarchies(C& criterion,
263 const std::shared_ptr<const Operator>& matrixptr,
264 const PI& pinfo);
271 struct LevelContext
272 {
289 typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
293 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
309 std::size_t level;
310 };
311
312
317 void mgc(LevelContext& levelContext);
318
319 void additiveMgc();
320
327 void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
328
333 bool moveToCoarseLevel(LevelContext& levelContext);
334
339 void initIteratorsWithFineLevel(LevelContext& levelContext);
340
342 std::shared_ptr<OperatorHierarchy> matrices_;
344 SmootherArgs smootherArgs_;
346 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
348 std::shared_ptr<CoarseSolver> solver_;
350 std::shared_ptr<Hierarchy<Range,A>> rhs_;
352 std::shared_ptr<Hierarchy<Domain,A>> lhs_;
354 std::shared_ptr<Hierarchy<Domain,A>> update_;
358 std::shared_ptr<ScalarProduct> scalarProduct_;
360 std::size_t gamma_;
362 std::size_t preSteps_;
364 std::size_t postSteps_;
365 bool buildHierarchy_;
366 bool additive;
367 bool coarsesolverconverged;
368 std::shared_ptr<Smoother> coarseSmoother_;
370 SolverCategory::Category category_;
372 std::size_t verbosity_;
373
374 struct ToLower
375 {
376 std::string operator()(const std::string& str)
377 {
378 std::stringstream retval;
379 std::ostream_iterator<char> out(retval);
380 std::transform(str.begin(), str.end(), out,
381 [](char c){
382 return std::tolower(c, std::locale::classic());
383 });
384 return retval.str();
385 }
386 };
387 };
388
389 template<class M, class X, class S, class PI, class A>
390 inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
391 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
392 smoothers_(amg.smoothers_), solver_(amg.solver_),
393 rhs_(), lhs_(), update_(),
394 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
395 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
396 buildHierarchy_(amg.buildHierarchy_),
397 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
398 coarseSmoother_(amg.coarseSmoother_),
399 category_(amg.category_),
400 verbosity_(amg.verbosity_)
401 {}
402
403 template<class M, class X, class S, class PI, class A>
405 const SmootherArgs& smootherArgs,
406 const Parameters& parms)
407 : matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
408 smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
409 rhs_(), lhs_(), update_(), scalarProduct_(0),
410 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
411 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
412 additive(parms.getAdditive()), coarsesolverconverged(true),
413 coarseSmoother_(),
414// #warning should category be retrieved from matrices?
415 category_(SolverCategory::category(*smoothers_->coarsest())),
416 verbosity_(parms.debugLevel())
417 {
418 assert(matrices_->isBuilt());
419
420 // build the necessary smoother hierarchies
421 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
422 }
423
424 template<class M, class X, class S, class PI, class A>
425 template<class C>
427 const C& criterion,
428 const SmootherArgs& smootherArgs,
429 const PI& pinfo)
430 : smootherArgs_(smootherArgs),
431 smoothers_(new Hierarchy<Smoother,A>), solver_(),
432 rhs_(), lhs_(), update_(), scalarProduct_(),
433 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
434 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
435 additive(criterion.getAdditive()), coarsesolverconverged(true),
436 coarseSmoother_(),
437 category_(SolverCategory::category(pinfo)),
438 verbosity_(criterion.debugLevel())
439 {
441 DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
442 // TODO: reestablish compile time checks.
443 //static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
444 // "Matrix and Solver must match in terms of category!");
445 auto matrixptr = stackobject_to_shared_ptr(matrix);
446 createHierarchies(criterion, matrixptr, pinfo);
447 }
448
449 template<class M, class X, class S, class PI, class A>
450 AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
451 const ParameterTree& configuration,
452 const ParallelInformation& pinfo) :
453 smoothers_(new Hierarchy<Smoother,A>),
454 solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
455 coarsesolverconverged(true), coarseSmoother_(),
456 category_(SolverCategory::category(pinfo))
457 {
458
459 if (configuration.hasKey ("smootherIterations"))
460 smootherArgs_.iterations = configuration.get<int>("smootherIterations");
461
462 if (configuration.hasKey ("smootherRelaxation"))
463 smootherArgs_.relaxationFactor = configuration.get<typename SmootherArgs::RelaxationFactor>("smootherRelaxation");
464
465 auto normName = ToLower()(configuration.get("strengthMeasure", "diagonal"));
466 auto index = configuration.get<int>("diagonalRowIndex", 0);
467
468 if ( normName == "diagonal")
469 {
470 using field_type = typename M::field_type;
471 using real_type = typename FieldTraits<field_type>::real_type;
472 std::is_convertible<field_type, real_type> compiles;
473
474 switch (index)
475 {
476 case 0:
477 createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<0>(), configuration, compiles);
478 break;
479 case 1:
480 createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<1>(), configuration, compiles);
481 break;
482 case 2:
483 createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<2>(), configuration, compiles);
484 break;
485 case 3:
486 createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<3>(), configuration, compiles);
487 break;
488 case 4:
489 createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<4>(), configuration, compiles);
490 break;
491 default:
492 DUNE_THROW(InvalidStateException, "Currently strengthIndex>4 is not supported.");
493 }
494 }
495 else if (normName == "rowsum")
496 createCriterionAndHierarchies(matrixptr, pinfo, RowSum(), configuration);
497 else if (normName == "frobenius")
498 createCriterionAndHierarchies(matrixptr, pinfo, FrobeniusNorm(), configuration);
499 else if (normName == "one")
500 createCriterionAndHierarchies(matrixptr, pinfo, AlwaysOneNorm(), configuration);
501 else
502 DUNE_THROW(Dune::NotImplemented, "Wrong config file: strengthMeasure "<<normName<<" is not supported by AMG");
503 }
504
505 template<class M, class X, class S, class PI, class A>
506 template<class Norm>
507 void AMG<M,X,S,PI,A>::createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const Norm&, const ParameterTree& configuration, std::false_type)
508 {
509 DUNE_THROW(InvalidStateException, "Strength of connection measure does not support this type ("
510 << className<typename M::field_type>() << ") as it is lacking a conversion to"
511 << className<typename FieldTraits<typename M::field_type>::real_type>() << ".");
512 }
513
514 template<class M, class X, class S, class PI, class A>
515 template<class Norm>
516 void AMG<M,X,S,PI,A>::createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const Norm&, const ParameterTree& configuration, std::true_type)
517 {
518 if (configuration.get<bool>("criterionSymmetric", true))
519 {
520 using Criterion = Dune::Amg::CoarsenCriterion<
522 Criterion criterion;
523 createHierarchies(criterion, matrixptr, pinfo, configuration);
524 }
525 else
526 {
527 using Criterion = Dune::Amg::CoarsenCriterion<
529 Criterion criterion;
530 createHierarchies(criterion, matrixptr, pinfo, configuration);
531 }
532 }
533
534 template<class M, class X, class S, class PI, class A>
535 template<class C>
536 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const ParameterTree& configuration)
537 {
538 if (configuration.hasKey ("maxLevel"))
539 criterion.setMaxLevel(configuration.get<int>("maxLevel"));
540
541 if (configuration.hasKey ("minCoarseningRate"))
542 criterion.setMinCoarsenRate(configuration.get<int>("minCoarseningRate"));
543
544 if (configuration.hasKey ("coarsenTarget"))
545 criterion.setCoarsenTarget (configuration.get<int>("coarsenTarget"));
546
547 if (configuration.hasKey ("accumulationMode"))
548 {
549 std::string mode = ToLower()(configuration.get<std::string>("accumulationMode"));
550 if ( mode == "none")
551 criterion.setAccumulate(AccumulationMode::noAccu);
552 else if ( mode == "atonce" )
553 criterion.setAccumulate(AccumulationMode::atOnceAccu);
554 else if ( mode == "successive")
555 criterion.setCoarsenTarget (AccumulationMode::successiveAccu);
556 else
557 DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
558 << mode <<".");
559 }
560
561 if (configuration.hasKey ("prolongationDampingFactor"))
562 criterion.setProlongationDampingFactor (configuration.get<double>("prolongationDampingFactor"));
563
564 if (configuration.hasKey("defaultAggregationSizeMode"))
565 {
566 auto mode = ToLower()(configuration.get<std::string>("defaultAggregationSizeMode"));
567 auto dim = configuration.get<std::size_t>("defaultAggregationDimension");
568 std::size_t maxDistance = 2;
569 if (configuration.hasKey("MaxAggregateDistance"))
570 maxDistance = configuration.get<std::size_t>("maxAggregateDistance");
571 if (mode == "isotropic")
572 criterion.setDefaultValuesIsotropic(dim, maxDistance);
573 else if(mode == "anisotropic")
574 criterion.setDefaultValuesAnisotropic(dim, maxDistance);
575 else
576 DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
577 << mode <<".");
578 }
579
580 if (configuration.hasKey("maxAggregateDistance"))
581 criterion.setMaxDistance(configuration.get<std::size_t>("maxAggregateDistance"));
582
583 if (configuration.hasKey("minAggregateSize"))
584 criterion.setMinAggregateSize(configuration.get<std::size_t>("minAggregateSize"));
585
586 if (configuration.hasKey("maxAggregateSize"))
587 criterion.setMaxAggregateSize(configuration.get<std::size_t>("maxAggregateSize"));
588
589 if (configuration.hasKey("maxAggregateConnectivity"))
590 criterion.setMaxConnectivity(configuration.get<std::size_t>("maxAggregateConnectivity"));
591
592 if (configuration.hasKey ("alpha"))
593 criterion.setAlpha (configuration.get<double> ("alpha"));
594
595 if (configuration.hasKey ("beta"))
596 criterion.setBeta (configuration.get<double> ("beta"));
597
598 if (configuration.hasKey ("gamma"))
599 criterion.setGamma (configuration.get<std::size_t> ("gamma"));
600 gamma_ = criterion.getGamma();
601
602 if (configuration.hasKey ("additive"))
603 criterion.setAdditive (configuration.get<bool>("additive"));
604 additive = criterion.getAdditive();
605
606 if (configuration.hasKey ("preSteps"))
607 criterion.setNoPreSmoothSteps (configuration.get<std::size_t> ("preSteps"));
608 preSteps_ = criterion.getNoPreSmoothSteps ();
609
610 if (configuration.hasKey ("postSteps"))
611 criterion.setNoPostSmoothSteps (configuration.get<std::size_t> ("preSteps"));
612 postSteps_ = criterion.getNoPostSmoothSteps ();
613
614 verbosity_ = configuration.get("verbosity", 0);
615 criterion.setDebugLevel (verbosity_);
616
617 createHierarchies(criterion, matrixptr, pinfo);
618 }
619
620 template <class Matrix,
621 class Vector>
623 {
626
627 static constexpr SolverType solver =
628#if DISABLE_AMG_DIRECTSOLVER
629 none;
630#elif HAVE_SUITESPARSE_UMFPACK
632#elif HAVE_SUPERLU
633 superlu ;
634#else
635 none;
636#endif
637
638 template <class M, SolverType>
639 struct Solver
640 {
642 static type* create(const M& mat, bool verbose, bool reusevector )
643 {
644 DUNE_THROW(NotImplemented,"DirectSolver not selected");
645 return nullptr;
646 }
647 static std::string name () { return "None"; }
648 };
649#if HAVE_SUITESPARSE_UMFPACK
650 template <class M>
651 struct Solver< M, umfpack >
652 {
653 typedef UMFPack< M > type;
654 static type* create(const M& mat, bool verbose, bool reusevector )
655 {
656 return new type(mat, verbose, reusevector );
657 }
658 static std::string name () { return "UMFPack"; }
659 };
660#endif
661#if HAVE_SUPERLU
662 template <class M>
663 struct Solver< M, superlu >
664 {
666 static type* create(const M& mat, bool verbose, bool reusevector )
667 {
668 return new type(mat, verbose, reusevector );
669 }
670 static std::string name () { return "SuperLU"; }
671 };
672#endif
673
674 // define direct solver type to be used
677 static constexpr bool isDirectSolver = solver != none;
678 static std::string name() { return SelectedSolver :: name (); }
679 static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
680 {
681 return SelectedSolver :: create( mat, verbose, reusevector );
682 }
683 };
684
685 template<class M, class X, class S, class PI, class A>
686 template<class C>
687 void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
688 const std::shared_ptr<const Operator>& matrixptr,
689 const PI& pinfo)
690 {
691 Timer watch;
692 matrices_ = std::make_shared<OperatorHierarchy>(
693 std::const_pointer_cast<Operator>(matrixptr),
694 stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
695
696 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
697
698 // build the necessary smoother hierarchies
699 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
700
701 // test whether we should solve on the coarse level. That is the case if we
702 // have that level and if there was a redistribution on this level then our
703 // communicator has to be valid (size()>0) as the smoother might try to communicate
704 // in the constructor.
705 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
706 && ( ! matrices_->redistributeInformation().back().isSetup() ||
707 matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
708 {
709 // We have the carsest level. Create the coarse Solver
710 SmootherArgs sargs(smootherArgs_);
711 sargs.iterations = 1;
712
714 cargs.setArgs(sargs);
715 if(matrices_->redistributeInformation().back().isSetup()) {
716 // Solve on the redistributed partitioning
717 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
718 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
719 }else{
720 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
721 cargs.setComm(*matrices_->parallelInformation().coarsest());
722 }
723
724 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
725 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
726
727 typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
728
729 // Use superlu if we are purely sequential or with only one processor on the coarsest level.
730 if( SolverSelector::isDirectSolver &&
731 (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
732 || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
733 || (matrices_->parallelInformation().coarsest().isRedistributed()
734 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
735 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
736 )
737 { // redistribute and 1 proc
738 if(matrices_->parallelInformation().coarsest().isRedistributed())
739 {
740 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
741 {
742 // We are still participating on this level
743 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
744 }
745 else
746 solver_.reset();
747 }
748 else
749 {
750 solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
751 }
752 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
753 std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
754 }
755 else
756 {
757 if(matrices_->parallelInformation().coarsest().isRedistributed())
758 {
759 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
760 // We are still participating on this level
761
762 // we have to allocate these types using the rebound allocator
763 // in order to ensure that we fulfill the alignment requirements
764 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
765 *scalarProduct_,
766 *coarseSmoother_, 1E-2, 1000, 0));
767 else
768 solver_.reset();
769 }else
770 {
771 solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
772 *scalarProduct_,
773 *coarseSmoother_, 1E-2, 1000, 0));
774 // // we have to allocate these types using the rebound allocator
775 // // in order to ensure that we fulfill the alignment requirements
776 // using Alloc = typename std::allocator_traits<A>::template rebind_alloc<BiCGSTABSolver<X>>;
777 // Alloc alloc;
778 // auto p = alloc.allocate(1);
779 // std::allocator_traits<Alloc>::construct(alloc, p,
780 // const_cast<M&>(*matrices_->matrices().coarsest()),
781 // *scalarProduct_,
782 // *coarseSmoother_, 1E-2, 1000, 0);
783 // solver_.reset(p,[](BiCGSTABSolver<X>* p){
784 // Alloc alloc;
785 // std::allocator_traits<Alloc>::destroy(alloc, p);
786 // alloc.deallocate(p,1);
787 // });
788 }
789 }
790 }
791
792 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
793 std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
794 <<"(including coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
795 }
796
797
798 template<class M, class X, class S, class PI, class A>
800 {
801 // Detect Matrix rows where all offdiagonal entries are
802 // zero and set x such that A_dd*x_d=b_d
803 // Thus users can be more careless when setting up their linear
804 // systems.
805 typedef typename M::matrix_type Matrix;
806 typedef typename Matrix::ConstRowIterator RowIter;
807 typedef typename Matrix::ConstColIterator ColIter;
808 typedef typename Matrix::block_type Block;
809 Block zero;
810 zero=typename Matrix::field_type();
811
812 const Matrix& mat=matrices_->matrices().finest()->getmat();
813 for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
814 bool isDirichlet = true;
815 bool hasDiagonal = false;
816 Block diagonal{};
817 for(ColIter col=row->begin(); col!=row->end(); ++col) {
818 if(row.index()==col.index()) {
819 diagonal = *col;
820 hasDiagonal = true;
821 }else{
822 if(*col!=zero)
823 isDirichlet = false;
824 }
825 }
826 if(isDirichlet && hasDiagonal)
827 {
828 auto&& xEntry = Impl::asVector(x[row.index()]);
829 auto&& bEntry = Impl::asVector(b[row.index()]);
830 Impl::asMatrix(diagonal).solve(xEntry, bEntry);
831 }
832 }
833
834 if(smoothers_->levels()>0)
835 smoothers_->finest()->pre(x,b);
836 else
837 // No smoother to make x consistent! Do it by hand
838 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
839 rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
840 lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
841 update_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
842 matrices_->coarsenVector(*rhs_);
843 matrices_->coarsenVector(*lhs_);
844 matrices_->coarsenVector(*update_);
845
846 // Preprocess all smoothers
847 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
848 typedef typename Hierarchy<Range,A>::Iterator RIterator;
849 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
850 Iterator coarsest = smoothers_->coarsest();
851 Iterator smoother = smoothers_->finest();
852 RIterator rhs = rhs_->finest();
853 DIterator lhs = lhs_->finest();
854 if(smoothers_->levels()>1) {
855
856 assert(lhs_->levels()==rhs_->levels());
857 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
858 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
859
860 if(smoother!=coarsest)
861 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
862 smoother->pre(*lhs,*rhs);
863 smoother->pre(*lhs,*rhs);
864 }
865
866
867 // The preconditioner might change x and b. So we have to
868 // copy the changes to the original vectors.
869 x = *lhs_->finest();
870 b = *rhs_->finest();
871
872 }
873 template<class M, class X, class S, class PI, class A>
875 {
876 return matrices_->levels();
877 }
878 template<class M, class X, class S, class PI, class A>
880 {
881 return matrices_->maxlevels();
882 }
883
885 template<class M, class X, class S, class PI, class A>
887 {
888 LevelContext levelContext;
889
890 if(additive) {
891 *(rhs_->finest())=d;
892 additiveMgc();
893 v=*lhs_->finest();
894 }else{
895 // Init all iterators for the current level
896 initIteratorsWithFineLevel(levelContext);
897
898
899 *levelContext.lhs = v;
900 *levelContext.rhs = d;
901 *levelContext.update=0;
902 levelContext.level=0;
903
904 mgc(levelContext);
905
906 if(postSteps_==0||matrices_->maxlevels()==1)
907 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
908
909 v=*levelContext.update;
910 }
911
912 }
913
914 template<class M, class X, class S, class PI, class A>
915 void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
916 {
917 levelContext.smoother = smoothers_->finest();
918 levelContext.matrix = matrices_->matrices().finest();
919 levelContext.pinfo = matrices_->parallelInformation().finest();
920 levelContext.redist =
921 matrices_->redistributeInformation().begin();
922 levelContext.aggregates = matrices_->aggregatesMaps().begin();
923 levelContext.lhs = lhs_->finest();
924 levelContext.update = update_->finest();
925 levelContext.rhs = rhs_->finest();
926 }
927
928 template<class M, class X, class S, class PI, class A>
929 bool AMG<M,X,S,PI,A>
930 ::moveToCoarseLevel(LevelContext& levelContext)
931 {
932
933 bool processNextLevel=true;
934
935 if(levelContext.redist->isSetup()) {
936 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
937 levelContext.rhs.getRedistributed());
938 processNextLevel = levelContext.rhs.getRedistributed().size()>0;
939 if(processNextLevel) {
940 //restrict defect to coarse level right hand side.
941 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
942 ++levelContext.pinfo;
944 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
945 static_cast<const Range&>(fineRhs.getRedistributed()),
946 *levelContext.pinfo);
947 }
948 }else{
949 //restrict defect to coarse level right hand side.
950 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
951 ++levelContext.pinfo;
953 ::restrictVector(*(*levelContext.aggregates),
954 *levelContext.rhs, static_cast<const Range&>(*fineRhs),
955 *levelContext.pinfo);
956 }
957
958 if(processNextLevel) {
959 // prepare coarse system
960 ++levelContext.lhs;
961 ++levelContext.update;
962 ++levelContext.matrix;
963 ++levelContext.level;
964 ++levelContext.redist;
965
966 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
967 // next level is not the globally coarsest one
968 ++levelContext.smoother;
969 ++levelContext.aggregates;
970 }
971 // prepare the update on the next level
972 *levelContext.update=0;
973 }
974 return processNextLevel;
975 }
976
977 template<class M, class X, class S, class PI, class A>
978 void AMG<M,X,S,PI,A>
979 ::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
980 {
981 if(processNextLevel) {
982 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
983 // previous level is not the globally coarsest one
984 --levelContext.smoother;
985 --levelContext.aggregates;
986 }
987 --levelContext.redist;
988 --levelContext.level;
989 //prolongate and add the correction (update is in coarse left hand side)
990 --levelContext.matrix;
991
992 //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
993 --levelContext.lhs;
994 --levelContext.pinfo;
995 }
996 if(levelContext.redist->isSetup()) {
997 // Need to redistribute during prolongateVector
998 levelContext.lhs.getRedistributed()=0;
1000 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1001 levelContext.lhs.getRedistributed(),
1002 matrices_->getProlongationDampingFactor(),
1003 *levelContext.pinfo, *levelContext.redist);
1004 }else{
1005 *levelContext.lhs=0;
1007 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
1008 matrices_->getProlongationDampingFactor(),
1009 *levelContext.pinfo);
1010 }
1011
1012
1013 if(processNextLevel) {
1014 --levelContext.update;
1015 --levelContext.rhs;
1016 }
1017
1018 *levelContext.update += *levelContext.lhs;
1019 }
1020
1021 template<class M, class X, class S, class PI, class A>
1023 {
1025 }
1026
1027 template<class M, class X, class S, class PI, class A>
1028 void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
1029 if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
1030 // Solve directly
1032 res.converged=true; // If we do not compute this flag will not get updated
1033 if(levelContext.redist->isSetup()) {
1034 levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
1035 if(levelContext.rhs.getRedistributed().size()>0) {
1036 // We are still participating in the computation
1037 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
1038 levelContext.rhs.getRedistributed());
1039 solver_->apply(levelContext.update.getRedistributed(),
1040 levelContext.rhs.getRedistributed(), res);
1041 }
1042 levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
1043 levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
1044 }else{
1045 levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
1046 solver_->apply(*levelContext.update, *levelContext.rhs, res);
1047 }
1048
1049 if (!res.converged)
1050 coarsesolverconverged = false;
1051 }else{
1052 // presmoothing
1053 presmooth(levelContext, preSteps_);
1054
1055#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1056 bool processNextLevel = moveToCoarseLevel(levelContext);
1057
1058 if(processNextLevel) {
1059 // next level
1060 for(std::size_t i=0; i<gamma_; i++)
1061 mgc(levelContext);
1062 }
1063
1064 moveToFineLevel(levelContext, processNextLevel);
1065#else
1066 *lhs=0;
1067#endif
1068
1069 if(levelContext.matrix == matrices_->matrices().finest()) {
1070 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
1071 if(!coarsesolverconverged)
1072 DUNE_THROW(MathError, "Coarse solver did not converge");
1073 }
1074 // postsmoothing
1075 postsmooth(levelContext, postSteps_);
1076
1077 }
1078 }
1079
1080 template<class M, class X, class S, class PI, class A>
1081 void AMG<M,X,S,PI,A>::additiveMgc(){
1082
1083 // restrict residual to all levels
1084 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
1085 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
1086 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
1087 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
1088
1089 for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
1090 ++pinfo;
1092 ::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
1093 }
1094
1095 // pinfo is invalid, set to coarsest level
1096 //pinfo = matrices_->parallelInformation().coarsest
1097 // calculate correction for all levels
1098 lhs = lhs_->finest();
1099 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
1100
1101 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
1102 // presmoothing
1103 *lhs=0;
1104 smoother->apply(*lhs, *rhs);
1105 }
1106
1107 // Coarse level solve
1108#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
1109 InverseOperatorResult res;
1110 pinfo->copyOwnerToAll(*rhs, *rhs);
1111 solver_->apply(*lhs, *rhs, res);
1112
1113 if(!res.converged)
1114 DUNE_THROW(MathError, "Coarse solver did not converge");
1115#else
1116 *lhs=0;
1117#endif
1118 // Prologate and add up corrections from all levels
1119 --pinfo;
1120 --aggregates;
1121
1122 for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
1124 ::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
1125 }
1126 }
1127
1128
1130 template<class M, class X, class S, class PI, class A>
1131 void AMG<M,X,S,PI,A>::post([[maybe_unused]] Domain& x)
1132 {
1133 // Postprocess all smoothers
1134 typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
1135 typedef typename Hierarchy<Domain,A>::Iterator DIterator;
1136 Iterator coarsest = smoothers_->coarsest();
1137 Iterator smoother = smoothers_->finest();
1138 DIterator lhs = lhs_->finest();
1139 if(smoothers_->levels()>0) {
1140 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
1141 smoother->post(*lhs);
1142 if(smoother!=coarsest)
1143 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
1144 smoother->post(*lhs);
1145 smoother->post(*lhs);
1146 }
1147 lhs_ = nullptr;
1148 update_ = nullptr;
1149 rhs_ = nullptr;
1150 }
1151
1152 template<class M, class X, class S, class PI, class A>
1153 template<class A1>
1154 void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
1155 {
1156 matrices_->getCoarsestAggregatesOnFinest(cont);
1157 }
1158
1159 } // end namespace Amg
1160
1162 template<class> struct isValidBlockType : std::false_type{};
1163 template<class T, int n, int m> struct isValidBlockType<FieldMatrix<T,n,m>> : std::true_type{};
1164
1165 template<class OP>
1166 std::shared_ptr<Dune::Preconditioner<typename OP::element_type::domain_type, typename OP::element_type::range_type> >
1167 makeAMG(const OP& op, const std::string& smoother, const Dune::ParameterTree& config) const
1168 {
1169 DUNE_THROW(Dune::Exception, "Operator type not supported by AMG");
1170 }
1171
1172 template<class M, class X, class Y>
1173 std::shared_ptr<Dune::Preconditioner<X,Y> >
1174 makeAMG(const std::shared_ptr<MatrixAdapter<M,X,Y>>& op, const std::string& smoother,
1175 const Dune::ParameterTree& config) const
1176 {
1177 using OP = MatrixAdapter<M,X,Y>;
1178
1179 if(smoother == "ssor")
1180 return std::make_shared<Amg::AMG<OP, X, SeqSSOR<M,X,Y>>>(op, config);
1181 if(smoother == "sor")
1182 return std::make_shared<Amg::AMG<OP, X, SeqSOR<M,X,Y>>>(op, config);
1183 if(smoother == "jac")
1184 return std::make_shared<Amg::AMG<OP, X, SeqJac<M,X,Y>>>(op, config);
1185 if(smoother == "gs")
1186 return std::make_shared<Amg::AMG<OP, X, SeqGS<M,X,Y>>>(op, config);
1187 if(smoother == "ilu")
1188 return std::make_shared<Amg::AMG<OP, X, SeqILU<M,X,Y>>>(op, config);
1189 else
1190 DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1191 }
1192
1193 template<class M, class X, class Y, class C>
1194 std::shared_ptr<Dune::Preconditioner<X,Y> >
1195 makeAMG(const std::shared_ptr<OverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1196 const Dune::ParameterTree& config) const
1197 {
1199
1200 auto cop = std::static_pointer_cast<const OP>(op);
1201
1202 if(smoother == "ssor")
1203 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1204 if(smoother == "sor")
1205 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
1206 if(smoother == "jac")
1207 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqJac<M,X,Y>>,C>>(cop, config, op->getCommunication());
1208 if(smoother == "gs")
1209 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqGS<M,X,Y>>,C>>(cop, config, op->getCommunication());
1210 if(smoother == "ilu")
1211 return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqILU<M,X,Y>>,C>>(cop, config, op->getCommunication());
1212 else
1213 DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1214 }
1215
1216 template<class M, class X, class Y, class C>
1217 std::shared_ptr<Dune::Preconditioner<X,Y> >
1218 makeAMG(const std::shared_ptr<NonoverlappingSchwarzOperator<M,X,Y,C>>& op, const std::string& smoother,
1219 const Dune::ParameterTree& config) const
1220 {
1222
1223 if(smoother == "ssor")
1224 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1225 if(smoother == "sor")
1226 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
1227 if(smoother == "jac")
1228 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqJac<M,X,Y>>,C>>(op, config, op->getCommunication());
1229 if(smoother == "gs")
1230 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqGS<M,X,Y>>,C>>(op, config, op->getCommunication());
1231 if(smoother == "ilu")
1232 return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqILU<M,X,Y>>,C>>(op, config, op->getCommunication());
1233 else
1234 DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
1235 }
1236
1237 template<typename TL, typename OP>
1238 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1239 typename Dune::TypeListElement<2, TL>::type>>
1240 operator() (TL tl, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
1242 {
1243 using field_type = typename OP::matrix_type::field_type;
1244 using real_type = typename FieldTraits<field_type>::real_type;
1245 if (!std::is_convertible<field_type, real_type>())
1246 DUNE_THROW(UnsupportedType, "AMG needs field_type(" <<
1247 className<field_type>() <<
1248 ") to be convertible to its real_type (" <<
1249 className<real_type>() <<
1250 ").");
1251 using D = typename Dune::TypeListElement<1, decltype(tl)>::type;
1252 using R = typename Dune::TypeListElement<2, decltype(tl)>::type;
1253 std::shared_ptr<Preconditioner<D,R>> amg;
1254 std::string smoother = config.get("smoother", "ssor");
1255 return makeAMG(op, smoother, config);
1256 }
1257
1258 template<typename TL, typename OP>
1259 std::shared_ptr<Dune::Preconditioner<typename Dune::TypeListElement<1, TL>::type,
1260 typename Dune::TypeListElement<2, TL>::type>>
1261 operator() (TL /*tl*/, const std::shared_ptr<OP>& /*mat*/, const Dune::ParameterTree& /*config*/,
1263 {
1264 DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
1265 }
1266 };
1267
1269} // end namespace Dune
1270
1271#endif
Provides a classes representing the hierarchies in AMG.
Classes for the generic construction and application of the smoothers.
Prolongation and restriction for amg.
Define base class for scalar product and norm.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Classes for using SuperLU with ISTL matrices.
Classes for using UMFPack with ISTL matrices.
Col col
Definition: matrixmatrix.hh:349
Matrix & mat
Definition: matrixmatrix.hh:345
AMG(const AMG &amg)
Copy constructor.
Definition: amg.hh:390
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: amg.hh:799
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition: amg.hh:679
static std::string name()
Definition: amg.hh:678
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition: amg.hh:301
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: amg.hh:305
static std::string name()
Definition: amg.hh:670
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: amg.hh:1022
X Domain
The domain type.
Definition: amg.hh:85
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:642
AMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition: amg.hh:404
AMG(std::shared_ptr< const Operator > fineOperator, const ParameterTree &configuration, const ParallelInformation &pinfo=ParallelInformation())
Constructor an AMG via ParameterTree.
Definition: amg.hh:450
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: amg.hh:285
SolverType
Definition: amg.hh:625
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: amg.hh:293
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition: amg.hh:98
Solver< Matrix, solver > SelectedSolver
Definition: amg.hh:675
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< MatrixAdapter< M, X, Y > > &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1174
std::string operator()(const std::string &str)
Definition: amg.hh:376
std::shared_ptr< Dune::Preconditioner< typename OP::element_type::domain_type, typename OP::element_type::range_type > > makeAMG(const OP &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1167
S Smoother
The type of the smoother.
Definition: amg.hh:95
static std::string name()
Definition: amg.hh:647
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition: amg.hh:277
M Operator
The matrix operator type.
Definition: amg.hh:71
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: amg.hh:281
static type * create(const M &mat, bool verbose, bool reusevector)
Definition: amg.hh:666
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: amg.hh:289
X Range
The range type.
Definition: amg.hh:87
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:404
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: amg.hh:1154
std::size_t levels()
Definition: amg.hh:874
InverseOperator< Vector, Vector > type
Definition: amg.hh:641
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< OverlappingSchwarzOperator< M, X, Y, C > > &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1195
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: amg.hh:297
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:42
static constexpr SolverType solver
Definition: amg.hh:627
static std::shared_ptr< T > construct(Arguments &args)
Construct an object with the specified arguments.
Definition: construction.hh:50
static constexpr bool isDirectSolver
Definition: amg.hh:677
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: amg.hh:219
Iterator coarsest()
Get an iterator positioned at the coarsest level.
Definition: hierarchy.hh:381
Matrix::field_type field_type
Definition: amg.hh:624
SelectedSolver::type DirectSolver
Definition: amg.hh:676
std::shared_ptr< Dune::Preconditioner< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL tl, const std::shared_ptr< OP > &op, const Dune::ParameterTree &config, std::enable_if_t< isValidBlockType< typename OP::matrix_type::block_type >::value, int >=0) const
Definition: amg.hh:1240
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: amg.hh:82
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: amg.hh:89
void post(Domain &x)
Clean up.
Definition: amg.hh:1131
std::size_t maxlevels()
Definition: amg.hh:879
std::shared_ptr< Dune::Preconditioner< X, Y > > makeAMG(const std::shared_ptr< NonoverlappingSchwarzOperator< M, X, Y, C > > &op, const std::string &smoother, const Dune::ParameterTree &config) const
Definition: amg.hh:1218
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:426
std::size_t level
The level index.
Definition: amg.hh:309
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:426
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: amg.hh:886
Smoother SmootherType
Definition: amg.hh:273
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: amg.hh:80
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: amg.hh:192
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: amg.hh:78
@ none
Definition: amg.hh:625
@ umfpack
Definition: amg.hh:625
@ superlu
Definition: amg.hh:625
@ atOnceAccu
Accumulate data to one process at once.
Definition: parameters.hh:242
@ noAccu
No data accumulution.
Definition: parameters.hh:236
@ successiveAccu
Successively accumulate to fewer processes.
Definition: parameters.hh:246
Definition: allocator.hh:9
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator())
ConstIterator class for sequential access.
Definition: matrix.hh:402
A generic dynamic dense matrix.
Definition: matrix.hh:559
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:563
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:587
T block_type
Export the type representing the components.
Definition: matrix.hh:566
Definition: matrixutils.hh:25
A nonoverlapping operator with communication object.
Definition: novlpschwarz.hh:62
Adapter to turn a matrix into a linear operator.
Definition: operators.hh:135
Functor using the row sum (infinity) norm to determine strong couplings.
Definition: aggregates.hh:461
Definition: aggregates.hh:478
Definition: aggregates.hh:494
Criterion taking advantage of symmetric matrices.
Definition: aggregates.hh:517
Criterion suitable for unsymmetric matrices.
Definition: aggregates.hh:537
an algebraic multigrid method using a Krylov-cycle.
Definition: kamg.hh:138
Two grid operator for AMG with Krylov cycle.
Definition: kamg.hh:31
Parallel algebraic multigrid based on agglomeration.
Definition: amg.hh:63
Definition: amg.hh:623
Definition: amg.hh:1161
Definition: amg.hh:1162
An overlapping Schwarz operator.
Definition: schwarz.hh:76
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:214
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:217
The hierarchies build by the coarsening process.
Definition: matrixhierarchy.hh:59
The criterion describing the stop criteria for the coarsening process.
Definition: matrixhierarchy.hh:281
All parameters for AMG.
Definition: parameters.hh:391
Definition: pinfo.hh:26
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:64
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
X::field_type field_type
The field type of the preconditioner.
Definition: preconditioner.hh:37
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:50
Statistics about the application of an inverse operator.
Definition: solver.hh:46
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Categories for the solvers.
Definition: solvercategory.hh:20
Category
Definition: solvercategory.hh:21
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:32
Definition: solvercategory.hh:52
Definition: solverregistry.hh:75
Definition: solvertype.hh:14
SuperLu Solver.
Definition: superlu.hh:269
Definition: umfpack.hh:47
The UMFPack direct sparse solver.
Definition: umfpack.hh:213