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