Loading...
Searching...
No Matches
amg.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_AMG_AMG_HH
#define DUNE_AMG_AMG_HH
#include <memory>
#include <sstream>
#include <dune/common/exceptions.hh>
#include <dune/istl/paamg/smoother.hh>
#include <dune/istl/paamg/transfer.hh>
#include <dune/istl/paamg/matrixhierarchy.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/superlu.hh>
#include <dune/istl/umfpack.hh>
#include <dune/istl/solvertype.hh>
#include <dune/istl/solverregistry.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include <dune/common/parametertree.hh>
namespace Dune
{
namespace Amg
{
template<class M, class X, class S, class P, class K, class A>
class KAMG;
template<class T>
class KAmgTwoGrid;
template<class M, class X, class S, class PI=SequentialInformation,
class A=std::allocator<X> >
class AMG : public Preconditioner<X,X>
{
template<class M1, class X1, class S1, class P1, class K1, class A1>
friend class KAMG;
friend class KAmgTwoGrid<AMG>;
public:
typedef M Operator;
typedef PI ParallelInformation;
typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
typedef X Domain;
typedef X Range;
typedef InverseOperator<X,X> CoarseSolver;
typedef S Smoother;
typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
AMG(OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs, const Parameters& parms);
template<class C>
AMG(const Operator& fineOperator, const C& criterion,
const SmootherArgs& smootherArgs=SmootherArgs(),
const ParallelInformation& pinfo=ParallelInformation());
AMG(std::shared_ptr<const Operator> fineOperator, const ParameterTree& configuration, const ParallelInformation& pinfo=ParallelInformation());
AMG(const AMG& amg);
void pre(Domain& x, Range& b);
void apply(Domain& v, const Range& d);
virtual SolverCategory::Category category() const
{
return category_;
}
void post(Domain& x);
template<class A1>
void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont);
std::size_t levels();
std::size_t maxlevels();
void recalculateHierarchy()
{
matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
}
bool usesDirectCoarseLevelSolver() const;
private:
/*
* @brief Helper function to create hierarchies with parameter tree.
*
* Will create the coarsen criterion with the norm and create the
* Hierarchies
* \tparam Norm Type of the norm to use.
*/
template<class Norm>
void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
const PI& pinfo, const Norm&,
const ParameterTree& configuration,
std::true_type compiles = std::true_type());
template<class Norm>
void createCriterionAndHierarchies(std::shared_ptr<const Operator> matrixptr,
const PI& pinfo, const Norm&,
const ParameterTree& configuration,
std::false_type);
template<class C>
void createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr,
const PI& pinfo, const ParameterTree& configuration);
template<class C>
void createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo);
struct LevelContext
{
typedef Smoother SmootherType;
typename Hierarchy<Smoother,A>::Iterator smoother;
typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix;
typename ParallelInformationHierarchy::Iterator pinfo;
typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
typename Hierarchy<Domain,A>::Iterator lhs;
typename Hierarchy<Domain,A>::Iterator update;
typename Hierarchy<Range,A>::Iterator rhs;
std::size_t level;
};
void mgc(LevelContext& levelContext);
void additiveMgc();
void moveToFineLevel(LevelContext& levelContext,bool processedFineLevel);
bool moveToCoarseLevel(LevelContext& levelContext);
void initIteratorsWithFineLevel(LevelContext& levelContext);
std::shared_ptr<OperatorHierarchy> matrices_;
SmootherArgs smootherArgs_;
std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
std::shared_ptr<CoarseSolver> solver_;
std::shared_ptr<Hierarchy<Range,A>> rhs_;
std::shared_ptr<Hierarchy<Domain,A>> lhs_;
std::shared_ptr<Hierarchy<Domain,A>> update_;
using ScalarProduct = Dune::ScalarProduct<X>;
std::shared_ptr<ScalarProduct> scalarProduct_;
std::size_t gamma_;
std::size_t preSteps_;
std::size_t postSteps_;
bool buildHierarchy_;
bool additive;
bool coarsesolverconverged;
std::shared_ptr<Smoother> coarseSmoother_;
SolverCategory::Category category_;
std::size_t verbosity_;
struct ToLower
{
std::string operator()(const std::string& str)
{
std::stringstream retval;
std::ostream_iterator<char> out(retval);
std::transform(str.begin(), str.end(), out,
[](char c){
return std::tolower(c, std::locale::classic());
});
return retval.str();
}
};
};
template<class M, class X, class S, class PI, class A>
inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
: matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
smoothers_(amg.smoothers_), solver_(amg.solver_),
rhs_(), lhs_(), update_(),
scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
buildHierarchy_(amg.buildHierarchy_),
additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
coarseSmoother_(amg.coarseSmoother_),
category_(amg.category_),
verbosity_(amg.verbosity_)
{}
template<class M, class X, class S, class PI, class A>
AMG<M,X,S,PI,A>::AMG(OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs,
const Parameters& parms)
: matrices_(stackobject_to_shared_ptr(matrices)), smootherArgs_(smootherArgs),
smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
rhs_(), lhs_(), update_(), scalarProduct_(0),
gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
additive(parms.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(),
// #warning should category be retrieved from matrices?
category_(SolverCategory::category(*smoothers_->coarsest())),
verbosity_(parms.debugLevel())
{
assert(matrices_->isBuilt());
// build the necessary smoother hierarchies
matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
}
template<class M, class X, class S, class PI, class A>
template<class C>
AMG<M,X,S,PI,A>::AMG(const Operator& matrix,
const C& criterion,
const SmootherArgs& smootherArgs,
const PI& pinfo)
: smootherArgs_(smootherArgs),
smoothers_(new Hierarchy<Smoother,A>), solver_(),
rhs_(), lhs_(), update_(), scalarProduct_(),
gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
additive(criterion.getAdditive()), coarsesolverconverged(true),
coarseSmoother_(),
category_(SolverCategory::category(pinfo)),
verbosity_(criterion.debugLevel())
{
if(SolverCategory::category(matrix) != SolverCategory::category(pinfo))
DUNE_THROW(InvalidSolverCategory, "Matrix and Communication must have the same SolverCategory!");
// TODO: reestablish compile time checks.
//static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
// "Matrix and Solver must match in terms of category!");
auto matrixptr = stackobject_to_shared_ptr(matrix);
createHierarchies(criterion, matrixptr, pinfo);
}
template<class M, class X, class S, class PI, class A>
AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
const ParameterTree& configuration,
const ParallelInformation& pinfo) :
smoothers_(new Hierarchy<Smoother,A>),
solver_(), rhs_(), lhs_(), update_(), scalarProduct_(), buildHierarchy_(true),
coarsesolverconverged(true), coarseSmoother_(),
category_(SolverCategory::category(pinfo))
{
if (configuration.hasKey ("smootherIterations"))
smootherArgs_.iterations = configuration.get<int>("smootherIterations");
if (configuration.hasKey ("smootherRelaxation"))
smootherArgs_.relaxationFactor = configuration.get<typename SmootherArgs::RelaxationFactor>("smootherRelaxation");
auto normName = ToLower()(configuration.get("strengthMeasure", "diagonal"));
auto index = configuration.get<int>("diagonalRowIndex", 0);
if ( normName == "diagonal")
{
using field_type = typename M::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
std::is_convertible<field_type, real_type> compiles;
switch (index)
{
case 0:
createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<0>(), configuration, compiles);
break;
case 1:
createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<1>(), configuration, compiles);
break;
case 2:
createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<2>(), configuration, compiles);
break;
case 3:
createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<3>(), configuration, compiles);
break;
case 4:
createCriterionAndHierarchies(matrixptr, pinfo, Diagonal<4>(), configuration, compiles);
break;
default:
DUNE_THROW(InvalidStateException, "Currently strengthIndex>4 is not supported.");
}
}
else if (normName == "rowsum")
createCriterionAndHierarchies(matrixptr, pinfo, RowSum(), configuration);
else if (normName == "frobenius")
createCriterionAndHierarchies(matrixptr, pinfo, FrobeniusNorm(), configuration);
else if (normName == "one")
createCriterionAndHierarchies(matrixptr, pinfo, AlwaysOneNorm(), configuration);
else
DUNE_THROW(Dune::NotImplemented, "Wrong config file: strengthMeasure "<<normName<<" is not supported by AMG");
}
template<class M, class X, class S, class PI, class A>
template<class Norm>
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)
{
DUNE_THROW(InvalidStateException, "Strength of connection measure does not support this type ("
<< className<typename M::field_type>() << ") as it is lacking a conversion to"
<< className<typename FieldTraits<typename M::field_type>::real_type>() << ".");
}
template<class M, class X, class S, class PI, class A>
template<class Norm>
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)
{
if (configuration.get<bool>("criterionSymmetric", true))
{
using Criterion = Dune::Amg::CoarsenCriterion<
Dune::Amg::SymmetricCriterion<typename M::matrix_type,Norm> >;
Criterion criterion;
createHierarchies(criterion, matrixptr, pinfo, configuration);
}
else
{
using Criterion = Dune::Amg::CoarsenCriterion<
Dune::Amg::UnSymmetricCriterion<typename M::matrix_type,Norm> >;
Criterion criterion;
createHierarchies(criterion, matrixptr, pinfo, configuration);
}
}
template<class M, class X, class S, class PI, class A>
template<class C>
void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, std::shared_ptr<const Operator> matrixptr, const PI& pinfo, const ParameterTree& configuration)
{
if (configuration.hasKey ("maxLevel"))
criterion.setMaxLevel(configuration.get<int>("maxLevel"));
if (configuration.hasKey ("minCoarseningRate"))
criterion.setMinCoarsenRate(configuration.get<int>("minCoarseningRate"));
if (configuration.hasKey ("coarsenTarget"))
criterion.setCoarsenTarget (configuration.get<int>("coarsenTarget"));
if (configuration.hasKey ("accumulationMode"))
{
std::string mode = ToLower()(configuration.get<std::string>("accumulationMode"));
if ( mode == "none")
criterion.setAccumulate(AccumulationMode::noAccu);
else if ( mode == "atonce" )
criterion.setAccumulate(AccumulationMode::atOnceAccu);
else if ( mode == "successive")
criterion.setCoarsenTarget (AccumulationMode::successiveAccu);
else
DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
<< mode <<".");
}
if (configuration.hasKey ("prolongationDampingFactor"))
criterion.setProlongationDampingFactor (configuration.get<double>("prolongationDampingFactor"));
if (configuration.hasKey("defaultAggregationSizeMode"))
{
auto mode = ToLower()(configuration.get<std::string>("defaultAggregationSizeMode"));
auto dim = configuration.get<std::size_t>("defaultAggregationDimension");
std::size_t maxDistance = 2;
if (configuration.hasKey("MaxAggregateDistance"))
maxDistance = configuration.get<std::size_t>("maxAggregateDistance");
if (mode == "isotropic")
criterion.setDefaultValuesIsotropic(dim, maxDistance);
else if(mode == "anisotropic")
criterion.setDefaultValuesAnisotropic(dim, maxDistance);
else
DUNE_THROW(InvalidSolverFactoryConfiguration, "Parameter accumulationMode does not allow value "
<< mode <<".");
}
if (configuration.hasKey("maxAggregateDistance"))
criterion.setMaxDistance(configuration.get<std::size_t>("maxAggregateDistance"));
if (configuration.hasKey("minAggregateSize"))
criterion.setMinAggregateSize(configuration.get<std::size_t>("minAggregateSize"));
if (configuration.hasKey("maxAggregateSize"))
criterion.setMaxAggregateSize(configuration.get<std::size_t>("maxAggregateSize"));
if (configuration.hasKey("maxAggregateConnectivity"))
criterion.setMaxConnectivity(configuration.get<std::size_t>("maxAggregateConnectivity"));
if (configuration.hasKey ("alpha"))
criterion.setAlpha (configuration.get<double> ("alpha"));
if (configuration.hasKey ("beta"))
criterion.setBeta (configuration.get<double> ("beta"));
if (configuration.hasKey ("gamma"))
criterion.setGamma (configuration.get<std::size_t> ("gamma"));
gamma_ = criterion.getGamma();
if (configuration.hasKey ("additive"))
criterion.setAdditive (configuration.get<bool>("additive"));
additive = criterion.getAdditive();
if (configuration.hasKey ("preSteps"))
criterion.setNoPreSmoothSteps (configuration.get<std::size_t> ("preSteps"));
preSteps_ = criterion.getNoPreSmoothSteps ();
if (configuration.hasKey ("postSteps"))
criterion.setNoPostSmoothSteps (configuration.get<std::size_t> ("postSteps"));
postSteps_ = criterion.getNoPostSmoothSteps ();
verbosity_ = configuration.get("verbosity", 0);
criterion.setDebugLevel (verbosity_);
createHierarchies(criterion, matrixptr, pinfo);
}
template <class Matrix,
class Vector>
struct DirectSolverSelector
{
typedef typename Matrix :: field_type field_type;
enum SolverType { umfpack, superlu, none };
static constexpr SolverType solver =
#if DISABLE_AMG_DIRECTSOLVER
none;
#elif HAVE_SUITESPARSE_UMFPACK
UMFPackMethodChooser< field_type > :: valid ? umfpack : none ;
#elif HAVE_SUPERLU
superlu ;
#else
none;
#endif
template <class M, SolverType>
struct Solver
{
typedef InverseOperator<Vector,Vector> type;
static type* create(const M& mat, bool verbose, bool reusevector )
{
DUNE_THROW(NotImplemented,"DirectSolver not selected");
return nullptr;
}
static std::string name () { return "None"; }
};
#if HAVE_SUITESPARSE_UMFPACK
template <class M>
struct Solver< M, umfpack >
{
typedef UMFPack< M > type;
static type* create(const M& mat, bool verbose, bool reusevector )
{
return new type(mat, verbose, reusevector );
}
static std::string name () { return "UMFPack"; }
};
#endif
#if HAVE_SUPERLU
template <class M>
struct Solver< M, superlu >
{
typedef SuperLU< M > type;
static type* create(const M& mat, bool verbose, bool reusevector )
{
return new type(mat, verbose, reusevector );
}
static std::string name () { return "SuperLU"; }
};
#endif
// define direct solver type to be used
typedef Solver< Matrix, solver > SelectedSolver ;
typedef typename SelectedSolver :: type DirectSolver;
static constexpr bool isDirectSolver = solver != none;
static std::string name() { return SelectedSolver :: name (); }
static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
{
return SelectedSolver :: create( mat, verbose, reusevector );
}
};
template<class M, class X, class S, class PI, class A>
template<class C>
void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo)
{
Timer watch;
matrices_ = std::make_shared<OperatorHierarchy>(
std::const_pointer_cast<Operator>(matrixptr),
stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
// build the necessary smoother hierarchies
matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
// test whether we should solve on the coarse level. That is the case if we
// have that level and if there was a redistribution on this level then our
// communicator has to be valid (size()>0) as the smoother might try to communicate
// in the constructor.
if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()
&& ( ! matrices_->redistributeInformation().back().isSetup() ||
matrices_->parallelInformation().coarsest().getRedistributed().communicator().size() ) )
{
// We have the carsest level. Create the coarse Solver
SmootherArgs sargs(smootherArgs_);
sargs.iterations = 1;
typename ConstructionTraits<Smoother>::Arguments cargs;
cargs.setArgs(sargs);
if(matrices_->redistributeInformation().back().isSetup()) {
// Solve on the redistributed partitioning
cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
}else{
cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
cargs.setComm(*matrices_->parallelInformation().coarsest());
}
coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
scalarProduct_ = createScalarProduct<X>(cargs.getComm(),category());
typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
// Use superlu if we are purely sequential or with only one processor on the coarsest level.
if( SolverSelector::isDirectSolver &&
(std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
|| matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
|| (matrices_->parallelInformation().coarsest().isRedistributed()
&& matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
&& matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
)
{ // redistribute and 1 proc
if(matrices_->parallelInformation().coarsest().isRedistributed())
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
{
// We are still participating on this level
solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
}
else
solver_.reset();
}
else
{
solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
}
if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
}
else
{
if(matrices_->parallelInformation().coarsest().isRedistributed())
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
// We are still participating on this level
// we have to allocate these types using the rebound allocator
// in order to ensure that we fulfill the alignment requirements
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
*scalarProduct_,
*coarseSmoother_, 1E-2, 1000, 0));
else
solver_.reset();
}else
{
solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
*scalarProduct_,
*coarseSmoother_, 1E-2, 1000, 0));
// // we have to allocate these types using the rebound allocator
// // in order to ensure that we fulfill the alignment requirements
// using Alloc = typename std::allocator_traits<A>::template rebind_alloc<BiCGSTABSolver<X>>;
// Alloc alloc;
// auto p = alloc.allocate(1);
// std::allocator_traits<Alloc>::construct(alloc, p,
// const_cast<M&>(*matrices_->matrices().coarsest()),
// *scalarProduct_,
// *coarseSmoother_, 1E-2, 1000, 0);
// solver_.reset(p,[](BiCGSTABSolver<X>* p){
// Alloc alloc;
// std::allocator_traits<Alloc>::destroy(alloc, p);
// alloc.deallocate(p,1);
// });
}
}
}
if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
std::cout<<"Building hierarchy of "<<matrices_->maxlevels()<<" levels "
<<"(including coarse solver) took "<<watch.elapsed()<<" seconds."<<std::endl;
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::pre(Domain& x, Range& b)
{
// Detect Matrix rows where all offdiagonal entries are
// zero and set x such that A_dd*x_d=b_d
// Thus users can be more careless when setting up their linear
// systems.
typedef typename M::matrix_type Matrix;
typedef typename Matrix::ConstRowIterator RowIter;
typedef typename Matrix::ConstColIterator ColIter;
typedef typename Matrix::block_type Block;
Block zero;
zero=typename Matrix::field_type();
const Matrix& mat=matrices_->matrices().finest()->getmat();
for(RowIter row=mat.begin(); row!=mat.end(); ++row) {
bool isDirichlet = true;
bool hasDiagonal = false;
Block diagonal{};
for(ColIter col=row->begin(); col!=row->end(); ++col) {
if(row.index()==col.index()) {
diagonal = *col;
hasDiagonal = true;
}else{
if(*col!=zero)
isDirichlet = false;
}
}
if(isDirichlet && hasDiagonal)
{
auto&& xEntry = Impl::asVector(x[row.index()]);
auto&& bEntry = Impl::asVector(b[row.index()]);
Impl::asMatrix(diagonal).solve(xEntry, bEntry);
}
}
if(smoothers_->levels()>0)
smoothers_->finest()->pre(x,b);
else
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
rhs_ = std::make_shared<Hierarchy<Range,A>>(std::make_shared<Range>(b));
lhs_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
update_ = std::make_shared<Hierarchy<Domain,A>>(std::make_shared<Domain>(x));
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*update_);
// Preprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Range,A>::Iterator RIterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_->coarsest();
Iterator smoother = smoothers_->finest();
RIterator rhs = rhs_->finest();
DIterator lhs = lhs_->finest();
if(smoothers_->levels()>1) {
assert(lhs_->levels()==rhs_->levels());
assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
if(smoother!=coarsest)
for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
smoother->pre(*lhs,*rhs);
smoother->pre(*lhs,*rhs);
}
// The preconditioner might change x and b. So we have to
// copy the changes to the original vectors.
x = *lhs_->finest();
b = *rhs_->finest();
}
template<class M, class X, class S, class PI, class A>
std::size_t AMG<M,X,S,PI,A>::levels()
{
return matrices_->levels();
}
template<class M, class X, class S, class PI, class A>
std::size_t AMG<M,X,S,PI,A>::maxlevels()
{
return matrices_->maxlevels();
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::apply(Domain& v, const Range& d)
{
LevelContext levelContext;
if(additive) {
*(rhs_->finest())=d;
additiveMgc();
v=*lhs_->finest();
}else{
// Init all iterators for the current level
initIteratorsWithFineLevel(levelContext);
*levelContext.lhs = v;
*levelContext.rhs = d;
*levelContext.update=0;
levelContext.level=0;
mgc(levelContext);
if(postSteps_==0||matrices_->maxlevels()==1)
levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
v=*levelContext.update;
}
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel(LevelContext& levelContext)
{
levelContext.smoother = smoothers_->finest();
levelContext.matrix = matrices_->matrices().finest();
levelContext.pinfo = matrices_->parallelInformation().finest();
levelContext.redist =
matrices_->redistributeInformation().begin();
levelContext.aggregates = matrices_->aggregatesMaps().begin();
levelContext.lhs = lhs_->finest();
levelContext.update = update_->finest();
levelContext.rhs = rhs_->finest();
}
template<class M, class X, class S, class PI, class A>
bool AMG<M,X,S,PI,A>
::moveToCoarseLevel(LevelContext& levelContext)
{
bool processNextLevel=true;
if(levelContext.redist->isSetup()) {
levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.rhs),
levelContext.rhs.getRedistributed());
processNextLevel = levelContext.rhs.getRedistributed().size()>0;
if(processNextLevel) {
//restrict defect to coarse level right hand side.
typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
++levelContext.pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
static_cast<const Range&>(fineRhs.getRedistributed()),
*levelContext.pinfo);
}
}else{
//restrict defect to coarse level right hand side.
typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
++levelContext.pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*levelContext.aggregates),
*levelContext.rhs, static_cast<const Range&>(*fineRhs),
*levelContext.pinfo);
}
if(processNextLevel) {
// prepare coarse system
++levelContext.lhs;
++levelContext.update;
++levelContext.matrix;
++levelContext.level;
++levelContext.redist;
if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
// next level is not the globally coarsest one
++levelContext.smoother;
++levelContext.aggregates;
}
// prepare the update on the next level
*levelContext.update=0;
}
return processNextLevel;
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>
::moveToFineLevel(LevelContext& levelContext, bool processNextLevel)
{
if(processNextLevel) {
if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
// previous level is not the globally coarsest one
--levelContext.smoother;
--levelContext.aggregates;
}
--levelContext.redist;
--levelContext.level;
//prolongate and add the correction (update is in coarse left hand side)
--levelContext.matrix;
//typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
--levelContext.lhs;
--levelContext.pinfo;
}
if(levelContext.redist->isSetup()) {
// Need to redistribute during prolongateVector
levelContext.lhs.getRedistributed()=0;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
levelContext.lhs.getRedistributed(),
matrices_->getProlongationDampingFactor(),
*levelContext.pinfo, *levelContext.redist);
}else{
*levelContext.lhs=0;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
matrices_->getProlongationDampingFactor(),
*levelContext.pinfo);
}
if(processNextLevel) {
--levelContext.update;
--levelContext.rhs;
}
*levelContext.update += *levelContext.lhs;
}
template<class M, class X, class S, class PI, class A>
bool AMG<M,X,S,PI,A>::usesDirectCoarseLevelSolver() const
{
return IsDirectSolver< CoarseSolver>::value;
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::mgc(LevelContext& levelContext){
if(levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
// Solve directly
InverseOperatorResult res;
res.converged=true; // If we do not compute this flag will not get updated
if(levelContext.redist->isSetup()) {
levelContext.redist->redistribute(*levelContext.rhs, levelContext.rhs.getRedistributed());
if(levelContext.rhs.getRedistributed().size()>0) {
// We are still participating in the computation
levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
levelContext.rhs.getRedistributed());
solver_->apply(levelContext.update.getRedistributed(),
levelContext.rhs.getRedistributed(), res);
}
levelContext.redist->redistributeBackward(*levelContext.update, levelContext.update.getRedistributed());
levelContext.pinfo->copyOwnerToAll(*levelContext.update, *levelContext.update);
}else{
levelContext.pinfo->copyOwnerToAll(*levelContext.rhs, *levelContext.rhs);
solver_->apply(*levelContext.update, *levelContext.rhs, res);
}
if (!res.converged)
coarsesolverconverged = false;
}else{
// presmoothing
presmooth(levelContext, preSteps_);
#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
bool processNextLevel = moveToCoarseLevel(levelContext);
if(processNextLevel) {
// next level
for(std::size_t i=0; i<gamma_; i++){
mgc(levelContext);
if (levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels())
break;
if(i+1 < gamma_){
levelContext.matrix->applyscaleadd(-1., *levelContext.lhs, *levelContext.rhs);
}
}
}
moveToFineLevel(levelContext, processNextLevel);
#else
*lhs=0;
#endif
if(levelContext.matrix == matrices_->matrices().finest()) {
coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
if(!coarsesolverconverged)
DUNE_THROW(MathError, "Coarse solver did not converge");
}
// postsmoothing
postsmooth(levelContext, postSteps_);
}
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::additiveMgc(){
// restrict residual to all levels
typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest();
typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest();
typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates) {
++pinfo;
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
}
// pinfo is invalid, set to coarsest level
//pinfo = matrices_->parallelInformation().coarsest
// calculate correction for all levels
lhs = lhs_->finest();
typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
// presmoothing
*lhs=0;
smoother->apply(*lhs, *rhs);
}
// Coarse level solve
#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
InverseOperatorResult res;
pinfo->copyOwnerToAll(*rhs, *rhs);
solver_->apply(*lhs, *rhs, res);
if(!res.converged)
DUNE_THROW(MathError, "Coarse solver did not converge");
#else
*lhs=0;
#endif
// Prologate and add up corrections from all levels
--pinfo;
--aggregates;
for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo) {
Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*aggregates), *coarseLhs, *lhs, 1.0, *pinfo);
}
}
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::post([[maybe_unused]] Domain& x)
{
// Postprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
Iterator coarsest = smoothers_->coarsest();
Iterator smoother = smoothers_->finest();
DIterator lhs = lhs_->finest();
if(smoothers_->levels()>0) {
if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
smoother->post(*lhs);
if(smoother!=coarsest)
for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
smoother->post(*lhs);
smoother->post(*lhs);
}
lhs_ = nullptr;
update_ = nullptr;
rhs_ = nullptr;
}
template<class M, class X, class S, class PI, class A>
template<class A1>
void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont)
{
matrices_->getCoarsestAggregatesOnFinest(cont);
}
} // end namespace Amg
struct AMGCreator{
template<class> struct isValidMatrix : std::false_type{};
template<class T, int n, int m, class A> struct isValidMatrix<BCRSMatrix<FieldMatrix<T,n,m>, A>> : std::true_type{};
template<class OP>
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
{
DUNE_THROW(Dune::Exception, "Operator type not supported by AMG");
}
template<class M, class X, class Y>
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
{
using OP = MatrixAdapter<M,X,Y>;
if(smoother == "ssor")
return std::make_shared<Amg::AMG<OP, X, SeqSSOR<M,X,Y>>>(op, config);
if(smoother == "sor")
return std::make_shared<Amg::AMG<OP, X, SeqSOR<M,X,Y>>>(op, config);
if(smoother == "jac")
return std::make_shared<Amg::AMG<OP, X, SeqJac<M,X,Y>>>(op, config);
if(smoother == "gs")
return std::make_shared<Amg::AMG<OP, X, SeqGS<M,X,Y>>>(op, config);
if(smoother == "ilu")
return std::make_shared<Amg::AMG<OP, X, SeqILU<M,X,Y>>>(op, config);
else
DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
}
template<class M, class X, class Y, class C>
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
{
using OP = OverlappingSchwarzOperator<M,X,Y,C>;
auto cop = std::static_pointer_cast<const OP>(op);
if(smoother == "ssor")
return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
if(smoother == "sor")
return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y>>,C>>(cop, config, op->getCommunication());
if(smoother == "jac")
return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqJac<M,X,Y>>,C>>(cop, config, op->getCommunication());
if(smoother == "gs")
return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqGS<M,X,Y>>,C>>(cop, config, op->getCommunication());
if(smoother == "ilu")
return std::make_shared<Amg::AMG<OP, X, BlockPreconditioner<X,Y,C,SeqILU<M,X,Y>>,C>>(cop, config, op->getCommunication());
else
DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
}
template<class M, class X, class Y, class C>
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
{
using OP = NonoverlappingSchwarzOperator<M,X,Y,C>;
if(smoother == "ssor")
return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
if(smoother == "sor")
return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y>>,C>>(op, config, op->getCommunication());
if(smoother == "jac")
return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqJac<M,X,Y>>,C>>(op, config, op->getCommunication());
if(smoother == "gs")
return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqGS<M,X,Y>>,C>>(op, config, op->getCommunication());
if(smoother == "ilu")
return std::make_shared<Amg::AMG<OP, X, NonoverlappingBlockPreconditioner<C,SeqILU<M,X,Y>>,C>>(op, config, op->getCommunication());
else
DUNE_THROW(Dune::Exception, "Unknown smoother for AMG");
}
template<typename OpTraits, typename OP>
std::shared_ptr<Dune::Preconditioner<typename OpTraits::domain_type,
typename OpTraits::range_type>>
operator() (OpTraits opTraits, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
std::enable_if_t<isValidMatrix<typename OpTraits::matrix_type>::value,int> = 0) const
{
using field_type = typename OpTraits::matrix_type::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
if (!std::is_convertible<field_type, real_type>())
DUNE_THROW(UnsupportedType, "AMG needs field_type(" <<
className<field_type>() <<
") to be convertible to its real_type (" <<
className<real_type>() <<
").");
std::string smoother = config.get("smoother", "ssor");
// we can irgnore the OpTraits here. As the AMG can only work
// with actual matrices, the operator op must be of type
// MatrixAdapter or *SchwarzOperator. In any case these
// operators provide all necessary information about matrix,
// domain and range type
return makeAMG(op, smoother, config);
}
template<typename OpTraits, typename OP>
std::shared_ptr<Dune::Preconditioner<typename OpTraits::domain_type,
typename OpTraits::range_type>>
operator() (OpTraits opTraits, const std::shared_ptr<OP>& op, const Dune::ParameterTree& config,
std::enable_if_t<!isValidMatrix<typename OpTraits::matrix_type>::value,int> = 0) const
{
DUNE_THROW(UnsupportedType, "AMG needs a FieldMatrix as Matrix block_type");
}
};
DUNE_REGISTER_PRECONDITIONER("amg", AMGCreator());
} // end namespace Dune
#endif
solvers.hh
Implementations of the inverse operator interface.
matrixhierarchy.hh
Provides a classes representing the hierarchies in AMG.
transfer.hh
Prolongation and restriction for amg.
smoother.hh
Classes for the generic construction and application of the smoothers.
superlu.hh
Classes for using SuperLU with ISTL matrices.
solvertype.hh
Templates characterizing the type of a solver.
scalarproducts.hh
Define base class for scalar product and norm.
umfpack.hh
Classes for using UMFPack with ISTL matrices.
solverregistry.hh
DUNE_REGISTER_PRECONDITIONER
#define DUNE_REGISTER_PRECONDITIONER(name,...)
Definition solverregistry.hh:13
Dune::Amg::AMG::AMG
AMG(const AMG &amg)
Copy constructor.
Definition amg.hh:393
Dune::Amg::AMG::pre
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition amg.hh:802
Dune::Amg::DirectSolverSelector::create
static DirectSolver * create(const Matrix &mat, bool verbose, bool reusevector)
Definition amg.hh:682
Dune::AMGCreator::operator()
std::shared_ptr< Dune::Preconditioner< typename OpTraits::domain_type, typename OpTraits::range_type > > operator()(OpTraits opTraits, const std::shared_ptr< OP > &op, const Dune::ParameterTree &config, std::enable_if_t< isValidMatrix< typename OpTraits::matrix_type >::value, int >=0) const
Definition amg.hh:1249
Dune::Amg::DirectSolverSelector::name
static std::string name()
Definition amg.hh:681
Dune::Amg::SmootherTraits::Arguments
DefaultSmootherArgs< typename T::matrix_type::field_type > Arguments
Definition smoother.hh:67
Dune::Amg::AMG::LevelContext::update
Hierarchy< Domain, A >::Iterator update
The iterator over the updates.
Definition amg.hh:304
Dune::Amg::AMG::LevelContext::rhs
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition amg.hh:308
Dune::Amg::DirectSolverSelector::Solver< M, superlu >::name
static std::string name()
Definition amg.hh:673
Dune::Amg::AMG::usesDirectCoarseLevelSolver
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition amg.hh:1025
Dune::Amg::AMG::Domain
X Domain
The domain type.
Definition amg.hh:88
Dune::Amg::ConstructionTraits::construct
static std::shared_ptr< T > construct(Arguments &)
Construct an object with the specified arguments.
Definition construction.hh:52
Dune::Amg::DirectSolverSelector::Solver::create
static type * create(const M &mat, bool verbose, bool reusevector)
Definition amg.hh:645
Dune::Amg::AMG::AMG
AMG(OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const SmootherArgs &smootherArgs, const Parameters &parms)
Construct a new amg with a specific coarse solver.
Definition amg.hh:407
Dune::Amg::AMG::AMG
AMG(std::shared_ptr< const Operator > fineOperator, const ParameterTree &configuration, const ParallelInformation &pinfo=ParallelInformation())
Constructor an AMG via ParameterTree.
Definition amg.hh:453
Dune::Amg::AMG::LevelContext::pinfo
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition amg.hh:288
Dune::Amg::DirectSolverSelector::SolverType
SolverType
Definition amg.hh:628
Dune::Amg::AMG::LevelContext::aggregates
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition amg.hh:296
Dune::Amg::AMG::SmootherArgs
SmootherTraits< Smoother >::Arguments SmootherArgs
The argument type for the construction of the smoother.
Definition amg.hh:101
Dune::Amg::DirectSolverSelector::SelectedSolver
Solver< Matrix, solver > SelectedSolver
Definition amg.hh:678
Dune::AMGCreator::makeAMG
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:1183
Dune::Amg::AMG::ToLower::operator()
std::string operator()(const std::string &str)
Definition amg.hh:379
Dune::AMGCreator::makeAMG
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:1176
Dune::Amg::AMG::Smoother
S Smoother
The type of the smoother.
Definition amg.hh:98
Dune::Amg::DirectSolverSelector::Solver::name
static std::string name()
Definition amg.hh:650
Dune::Amg::AMG::LevelContext::smoother
Hierarchy< Smoother, A >::Iterator smoother
The iterator over the smoothers.
Definition amg.hh:280
Dune::Amg::AMG::Operator
M Operator
The matrix operator type.
Definition amg.hh:74
Dune::Amg::AMG::LevelContext::matrix
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition amg.hh:284
Dune::Amg::DirectSolverSelector::Solver< M, superlu >::type
SuperLU< M > type
Definition amg.hh:668
Dune::Amg::DirectSolverSelector::Solver< M, superlu >::create
static type * create(const M &mat, bool verbose, bool reusevector)
Definition amg.hh:669
Dune::Amg::AMG::LevelContext::redist
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition amg.hh:292
Dune::Amg::AMG::Range
X Range
The range type.
Definition amg.hh:90
Dune::Amg::presmooth
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition smoother.hh:406
Dune::Amg::AMG::getCoarsestAggregateNumbers
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition amg.hh:1163
Dune::Amg::AMG::levels
std::size_t levels()
Definition amg.hh:877
Dune::Amg::DirectSolverSelector::Solver::type
InverseOperator< Vector, Vector > type
Definition amg.hh:644
Dune::AMGCreator::makeAMG
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:1204
Dune::Amg::AMG::LevelContext::lhs
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition amg.hh:300
Dune::Amg::ConstructionTraits::Arguments
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition construction.hh:44
Dune::Amg::DirectSolverSelector::solver
static constexpr SolverType solver
Definition amg.hh:630
Dune::Amg::DirectSolverSelector::isDirectSolver
static constexpr bool isDirectSolver
Definition amg.hh:680
Dune::Amg::AMG::recalculateHierarchy
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition amg.hh:222
Dune::Amg::DefaultSmootherArgs::RelaxationFactor
FieldTraits< T >::real_type RelaxationFactor
The type of the relaxation factor.
Definition smoother.hh:42
Dune::Amg::DirectSolverSelector::field_type
Matrix::field_type field_type
Definition amg.hh:627
Dune::Amg::DirectSolverSelector::DirectSolver
SelectedSolver::type DirectSolver
Definition amg.hh:679
Dune::Amg::AMG::ParallelInformationHierarchy
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition amg.hh:85
Dune::Amg::AMG::CoarseSolver
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition amg.hh:92
Dune::Amg::AMG::post
void post(Domain &x)
Clean up.
Definition amg.hh:1140
Dune::Amg::AMG::maxlevels
std::size_t maxlevels()
Definition amg.hh:882
Dune::AMGCreator::makeAMG
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:1227
Dune::Amg::postsmooth
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition smoother.hh:428
Dune::Amg::AMG::LevelContext::level
std::size_t level
The level index.
Definition amg.hh:312
Dune::Amg::AMG::AMG
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
Dune::Amg::AMG::apply
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition amg.hh:889
Dune::Amg::AMG::LevelContext::SmootherType
Smoother SmootherType
Definition amg.hh:276
Dune::Amg::AMG::OperatorHierarchy
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition amg.hh:83
Dune::Amg::AMG::category
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category).
Definition amg.hh:195
Dune::Amg::AMG< M, X, S, PI, A >::KAMG
friend class KAMG
Definition amg.hh:68
Dune::Amg::AMG::ParallelInformation
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition amg.hh:81
Dune::Amg::DirectSolverSelector::none
@ none
Definition amg.hh:628
Dune::Amg::DirectSolverSelector::umfpack
@ umfpack
Definition amg.hh:628
Dune::Amg::DirectSolverSelector::superlu
@ superlu
Definition amg.hh:628
Dune::Amg::atOnceAccu
@ atOnceAccu
Accumulate data to one process at once.
Definition parameters.hh:243
Dune::Amg::noAccu
@ noAccu
No data accumulution.
Definition parameters.hh:237
Dune::Amg::successiveAccu
@ successiveAccu
Successively accumulate to fewer processes.
Definition parameters.hh:247
Dune
Definition allocator.hh:11
Dune::createScalarProduct
std::shared_ptr< ScalarProduct< X > > createScalarProduct(const Comm &comm, SolverCategory::Category category)
Definition scalarproducts.hh:242
Dune::col
MatrixSparsityPatternGatherScatter< M, I >::ColIter MatrixSparsityPatternGatherScatter< M, I >::col
Definition matrixredistribute.hh:591
Dune::type
const Dtype_t BaseGetSuperLUType< T >::type
Definition supermatrix.hh:135
Dune::Amg
Definition novlpschwarz.hh:256
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:467
Dune::Matrix
A generic dynamic dense matrix.
Definition matrix.hh:561
Dune::Matrix::ConstRowIterator
MatrixImp::DenseMatrixBase< T, A >::ConstIterator ConstRowIterator
Const iterator over the matrix rows.
Definition matrix.hh:586
Dune::Matrix::end
RowIterator end()
Get iterator to one beyond last row.
Definition matrix.hh:616
Dune::Matrix::begin
RowIterator begin()
Get iterator to first row.
Definition matrix.hh:610
Dune::Matrix::field_type
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition matrix.hh:565
Dune::Matrix::ConstColIterator
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition matrix.hh:589
Dune::Matrix::block_type
T block_type
Export the type representing the components.
Definition matrix.hh:568
Dune::FieldMatrix
Definition matrixutils.hh:27
Dune::NonoverlappingSchwarzOperator
A nonoverlapping operator with communication object.
Definition novlpschwarz.hh:61
Dune::MatrixAdapter
Adapter to turn a matrix into a linear operator.
Definition operators.hh:135
Dune::Amg::Diagonal
Norm that uses only the [N][N] entry of the block to determine couplings.
Definition aggregates.hh:381
Dune::Amg::RowSum
Functor using the row sum (infinity) norm to determine strong couplings.
Definition aggregates.hh:465
Dune::Amg::FrobeniusNorm
Definition aggregates.hh:486
Dune::Amg::AlwaysOneNorm
Definition aggregates.hh:502
Dune::Amg::SymmetricCriterion
Criterion taking advantage of symmetric matrices.
Definition aggregates.hh:525
Dune::Amg::UnSymmetricCriterion
Criterion suitable for unsymmetric matrices.
Definition aggregates.hh:545
Dune::Amg::KAMG
an algebraic multigrid method using a Krylov-cycle.
Definition kamg.hh:140
Dune::Amg::KAmgTwoGrid
Two grid operator for AMG with Krylov cycle.
Definition kamg.hh:33
Dune::Amg::AMG
Parallel algebraic multigrid based on agglomeration.
Definition amg.hh:66
Dune::Amg::DirectSolverSelector
Definition amg.hh:626
Dune::Amg::DirectSolverSelector::Solver
Definition amg.hh:643
Dune::AMGCreator
Definition amg.hh:1170
Dune::AMGCreator::isValidMatrix
Definition amg.hh:1171
Dune::OverlappingSchwarzOperator
An overlapping Schwarz operator.
Definition schwarz.hh:75
Dune::Amg::Hierarchy
A hierarchy of containers (e.g. matrices or vectors).
Definition hierarchy.hh:40
Dune::Amg::Hierarchy::Iterator
LevelIterator< Hierarchy< T, A >, T > Iterator
Type of the mutable iterator.
Definition hierarchy.hh:220
Dune::Amg::MatrixHierarchy
The hierarchies build by the coarsening process.
Definition matrixhierarchy.hh:61
Dune::Amg::MatrixHierarchy< M, ParallelInformation, A >::ParallelInformationHierarchy
Dune::Amg::Hierarchy< ParallelInformation, Allocator > ParallelInformationHierarchy
Definition matrixhierarchy.hh:82
Dune::Amg::CoarsenCriterion
The criterion describing the stop criteria for the coarsening process.
Definition matrixhierarchy.hh:283
Dune::Amg::Parameters
All parameters for AMG.
Definition parameters.hh:416
Dune::Amg::SequentialInformation
Definition pinfo.hh:28
Dune::Amg::Transfer
Definition transfer.hh:32
Dune::Preconditioner
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:33
Dune::Preconditioner< X, X >::field_type
X::field_type field_type
Definition preconditioner.hh:40
Dune::ScalarProduct
Base class for scalar product and norm computation.
Definition scalarproducts.hh:52
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition solver.hh:50
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition solver.hh:75
Dune::InverseOperator
Abstract base class for all solvers.
Definition solver.hh:101
Dune::SolverCategory
Categories for the solvers.
Definition solvercategory.hh:22
Dune::SolverCategory::Category
Category
Definition solvercategory.hh:23
Dune::SolverCategory::category
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
Dune::InvalidSolverCategory
Definition solvercategory.hh:54
Dune::UnsupportedType
Definition solverregistry.hh:97
Dune::IsDirectSolver::value
@ value
Whether this is a direct solver.
Definition solvertype.hh:24
Dune::SuperLU
SuperLu Solver.
Definition superlu.hh:271
Dune::UMFPackMethodChooser
Definition umfpack.hh:53
Dune::UMFPack
The UMFPack direct sparse solver.
Definition umfpack.hh:265