3#ifndef DUNE_ISTL_FASTAMG_HH
4#define DUNE_ISTL_FASTAMG_HH
7#include <dune/common/exceptions.hh>
8#include <dune/common/typetraits.hh>
9#include <dune/common/unused.hh>
57 template<
class M,
class X,
class PI=SequentialInformation,
class A=std::allocator<X> >
183 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
191 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
195 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
252 void initIteratorsWithFineLevel(LevelContext&
levelContext);
255 std::shared_ptr<OperatorHierarchy> matrices_;
257 std::shared_ptr<CoarseSolver> solver_;
268 typedef typename ScalarProductChooserType::ScalarProduct ScalarProduct;
269 typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
271 ScalarProductPointer scalarProduct_;
275 std::size_t preSteps_;
277 std::size_t postSteps_;
279 bool buildHierarchy_;
281 bool coarsesolverconverged;
283 typedef std::shared_ptr<Smoother> SmootherPointer;
284 SmootherPointer coarseSmoother_;
286 std::size_t verbosity_;
289 template<
class M,
class X,
class PI,
class A>
291 : matrices_(amg.matrices_), solver_(amg.solver_),
292 rhs_(), lhs_(), residual_(), scalarProduct_(amg.scalarProduct_),
293 gamma_(amg.gamma_), preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
294 symmetric(amg.symmetric), coarsesolverconverged(amg.coarsesolverconverged),
295 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
305 template<
class M,
class X,
class PI,
class A>
308 : matrices_(&matrices), solver_(&coarseSolver),
309 rhs_(), lhs_(), residual_(), scalarProduct_(),
310 gamma_(
parms.getGamma()), preSteps_(
parms.getNoPreSmoothSteps()),
311 postSteps_(
parms.getNoPostSmoothSteps()), buildHierarchy_(
false),
313 coarseSmoother_(), verbosity_(
parms.debugLevel())
315 if(preSteps_>1||postSteps_>1)
317 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
318 preSteps_=postSteps_=0;
320 assert(matrices_->isBuilt());
321 static_assert(std::is_same<PI,SequentialInformation>::value,
322 "Currently only sequential runs are supported");
324 template<
class M,
class X,
class PI,
class A>
331 : solver_(), rhs_(), lhs_(), residual_(), scalarProduct_(), gamma_(
parms.getGamma()),
332 preSteps_(
parms.getNoPreSmoothSteps()), postSteps_(
parms.getNoPostSmoothSteps()),
333 buildHierarchy_(
true),
335 coarseSmoother_(), verbosity_(
criterion.debugLevel())
337 if(preSteps_>1||postSteps_>1)
339 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
340 preSteps_=postSteps_=1;
342 static_assert(std::is_same<PI,SequentialInformation>::value,
343 "Currently only sequential runs are supported");
350 template<
class M,
class X,
class PI,
class A>
353 if(buildHierarchy_) {
357 coarseSmoother_.reset();
370 template<
class M,
class X,
class PI,
class A>
376 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
380 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
381 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<
watch.elapsed()<<
" seconds."<<std::endl;
383 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
387 sargs.iterations = 1;
391 if(matrices_->redistributeInformation().back().isSetup()) {
393 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
394 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
396 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
397 cargs.setComm(*matrices_->parallelInformation().coarsest());
401 scalarProduct_.reset(ScalarProductChooserType::construct(cargs.getComm()));
403#if HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK
404#if HAVE_SUITESPARSE_UMFPACK
405#define DIRECTSOLVER UMFPack
407#define DIRECTSOLVER SuperLU
410 if(std::is_same<ParallelInformation,SequentialInformation>::value
411 || matrices_->parallelInformation().coarsest()->communicator().size()==1
412 || (matrices_->parallelInformation().coarsest().isRedistributed()
413 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
414 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
415 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
416 std::cout<<
"Using superlu"<<std::endl;
417 if(matrices_->parallelInformation().coarsest().isRedistributed())
419 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
421 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
425 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
430 if(matrices_->parallelInformation().coarsest().isRedistributed())
432 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
434 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
436 *coarseSmoother_, 1E-2, 1000, 0));
440 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
442 *coarseSmoother_, 1E-2, 1000, 0));
446 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
447 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.elapsed()<<
" seconds."<<std::endl;
451 template<
class M,
class X,
class PI,
class A>
459 typedef typename M::matrix_type
Matrix;
466 const Matrix&
mat=matrices_->matrices().finest()->getmat();
471 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
472 if(row.index()==
col.index()) {
481 diag->solve(x[row.index()], b[row.index()]);
484 std::cout<<
" Preprocessing Dirichlet took "<<
watch1.elapsed()<<std::endl;
487 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
498 matrices_->coarsenVector(*rhs_);
499 matrices_->coarsenVector(*lhs_);
500 matrices_->coarsenVector(*residual_);
507 template<
class M,
class X,
class PI,
class A>
510 return matrices_->levels();
512 template<
class M,
class X,
class PI,
class A>
515 return matrices_->maxlevels();
519 template<
class M,
class X,
class PI,
class A>
529 if(matrices_->maxlevels()==1){
535 if(postSteps_==0||matrices_->maxlevels()==1)
539 template<
class M,
class X,
class PI,
class A>
543 levelContext.pinfo = matrices_->parallelInformation().finest();
545 matrices_->redistributeInformation().begin();
546 levelContext.aggregates = matrices_->aggregatesMaps().begin();
553 template<
class M,
class X,
class PI,
class A>
554 bool FastAMG<M,X,PI,A>
567 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
569 static_cast<const Range&
>(
levelContext.residual.getRedistributed()),
575 ++levelContext.pinfo;
576 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
577 ::restrictVector(*(*levelContext.aggregates), *levelContext.rhs,
578 static_cast<const Range&
>(*levelContext.residual), *levelContext.pinfo);
581 if(processNextLevel) {
583 ++levelContext.residual;
585 ++levelContext.matrix;
586 ++levelContext.level;
587 ++levelContext.redist;
589 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
591 ++levelContext.aggregates;
595 *levelContext.residual=0;
597 return processNextLevel;
600 template<
class M,
class X,
class PI,
class A>
601 void FastAMG<M,X,PI,A>
602 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel, Domain& x)
604 if(processNextLevel) {
605 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
607 --levelContext.aggregates;
609 --levelContext.redist;
610 --levelContext.level;
612 --levelContext.matrix;
613 --levelContext.residual;
617 typename Hierarchy<Domain,A>::Iterator coarseLhs = levelContext.lhs--;
618 if(levelContext.redist->isSetup()) {
621 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
622 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
623 levelContext.lhs.getRedistributed(),
624 matrices_->getProlongationDampingFactor(),
625 *levelContext.pinfo, *levelContext.redist);
627 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
628 ::prolongateVector(*(*levelContext.aggregates), *coarseLhs, x,
629 matrices_->getProlongationDampingFactor(), *levelContext.pinfo);
635 if(processNextLevel) {
642 template<
class M,
class X,
class PI,
class A>
643 void FastAMG<M,X,PI,A>
644 ::presmooth(LevelContext& levelContext, Domain& x,
const Range& b)
648 *levelContext.residual,
652 template<
class M,
class X,
class PI,
class A>
653 void FastAMG<M,X,PI,A>
654 ::postsmooth(LevelContext& levelContext, Domain& x,
const Range& b)
656 GaussSeidelPostsmoothDefect<M::matrix_type::blocklevel>
657 ::apply(levelContext.matrix->getmat(), x, *levelContext.residual, b);
661 template<
class M,
class X,
class PI,
class A>
667 template<
class M,
class X,
class PI,
class A>
670 if(
levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
685 levelContext.pinfo->copyOwnerToAll(b, b);
686 solver_->apply(v,
const_cast<Range&
>(b), res);
692 coarsesolverconverged =
false;
698#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
699 bool processNextLevel = moveToCoarseLevel(levelContext);
701 if(processNextLevel) {
703 for(std::size_t i=0; i<gamma_; i++)
704 mgc(levelContext, *levelContext.lhs, *levelContext.rhs);
707 moveToFineLevel(levelContext, processNextLevel, v);
712 if(levelContext.matrix == matrices_->matrices().finest()) {
713 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
714 if(!coarsesolverconverged)
715 DUNE_THROW(MathError,
"Coarse solver did not converge");
728 template<
class M,
class X,
class PI,
class A>
740 template<
class M,
class X,
class PI,
class A>
744 matrices_->getCoarsestAggregatesOnFinest(
cont);
Classes for using UMFPack with ISTL matrices.
Classes for using SuperLU with ISTL matrices.
Define base class for scalar product and norm.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Define general preconditioner interface.
Some generic functions for pretty printing vectors and matrices.
Prolongation and restriction for amg.
Provides a classes representing the hierarchies in AMG.
Classes for the generic construction and application of the smoothers.
Col col
Definition matrixmatrix.hh:347
Matrix & mat
Definition matrixmatrix.hh:343
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition smoother.hh:408
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition smoother.hh:430
static T * construct(Arguments &args)
Construct an object with the specified arguments.
Definition construction.hh:52
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition fastamg.hh:183
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition fastamg.hh:742
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition fastamg.hh:187
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition fastamg.hh:150
void post(Domain &x)
Clean up.
Definition fastamg.hh:729
std::size_t maxlevels()
Definition fastamg.hh:513
X Domain
The domain type.
Definition fastamg.hh:76
Hierarchy< Domain, A >::Iterator residual
The iterator over the residuals.
Definition fastamg.hh:203
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition fastamg.hh:71
~FastAMG()
Definition fastamg.hh:351
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition fastamg.hh:191
X Range
The range type.
Definition fastamg.hh:78
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition fastamg.hh:69
M Operator
The matrix operator type.
Definition fastamg.hh:62
std::size_t levels()
Definition fastamg.hh:508
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition fastamg.hh:80
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition fastamg.hh:662
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition fastamg.hh:199
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition fastamg.hh:207
std::size_t level
The level index.
Definition fastamg.hh:211
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition fastamg.hh:520
FastAMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const Parameters &parms, bool symmetric=true)
Construct a new amg with a specific coarse solver.
Definition fastamg.hh:306
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition fastamg.hh:73
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition fastamg.hh:452
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition fastamg.hh:195
@ category
The solver category.
Definition fastamg.hh:84
Definition basearray.hh:19
Iterator implementation class
Definition basearray.hh:84
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
ConstIterator class for sequential access.
Definition matrix.hh:398
A generic dynamic dense matrix.
Definition matrix.hh:555
T::field_type field_type
Export the type representing the underlying field.
Definition matrix.hh:559
T block_type
Export the type representing the components.
Definition matrix.hh:562
A fast (sequential) algebraic multigrid based on agglomeration that saves memory bandwidth.
Definition fastamg.hh:59
static void apply(const M &A, X &x, Y &d, const Y &b)
Definition fastamgsmoother.hh:17
The hierarchies build by the coarsening process.
Definition hierarchy.hh:317
All parameters for AMG.
Definition parameters.hh:391
Base class for matrix free definition of preconditioners.
Definition preconditioner.hh:26
Statistics about the application of an inverse operator.
Definition solver.hh:32
bool converged
True if convergence criterion has been met.
Definition solver.hh:56
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:21
Definition solvertype.hh:14