DUNE PDELab (git)

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