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> >
89 FastAMG(
const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
105 FastAMG(
const Operator& fineOperator,
const C& criterion,
118 void pre(Domain& x, Range& b);
121 void apply(Domain& v,
const Range& d);
130 void post(Domain& x);
153 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
170 void createHierarchies(C& criterion, Operator& matrix,
192 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
196 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
216 void mgc(LevelContext& levelContext, Domain& x,
const Range& b);
224 void presmooth(LevelContext& levelContext, Domain& x,
const Range& b);
232 void postsmooth(LevelContext& levelContext, Domain& x,
const Range& b);
240 void moveToFineLevel(LevelContext& levelContext,
bool processedFineLevel,
247 bool moveToCoarseLevel(LevelContext& levelContext);
253 void initIteratorsWithFineLevel(LevelContext& levelContext);
256 std::shared_ptr<OperatorHierarchy> matrices_;
258 std::shared_ptr<CoarseSolver> solver_;
269 std::shared_ptr<ScalarProduct> scalarProduct_;
273 std::size_t preSteps_;
275 std::size_t postSteps_;
277 bool buildHierarchy_;
279 bool coarsesolverconverged;
281 typedef std::shared_ptr<Smoother> SmootherPointer;
282 SmootherPointer coarseSmoother_;
284 std::size_t verbosity_;
287 template<
class M,
class X,
class PI,
class A>
289 : matrices_(amg.matrices_), solver_(amg.solver_),
290 rhs_(), lhs_(), residual_(), scalarProduct_(amg.scalarProduct_),
291 gamma_(amg.gamma_), preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
292 symmetric(amg.symmetric), coarsesolverconverged(amg.coarsesolverconverged),
293 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
303 template<
class M,
class X,
class PI,
class A>
306 : matrices_(&matrices), solver_(&coarseSolver),
307 rhs_(), lhs_(), residual_(), scalarProduct_(),
308 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
309 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
310 symmetric(symmetric_), coarsesolverconverged(true),
311 coarseSmoother_(), verbosity_(parms.debugLevel())
313 if(preSteps_>1||postSteps_>1)
315 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
316 preSteps_=postSteps_=0;
318 assert(matrices_->isBuilt());
319 static_assert(std::is_same<PI,SequentialInformation>::value,
320 "Currently only sequential runs are supported");
322 template<
class M,
class X,
class PI,
class A>
329 : solver_(), rhs_(), lhs_(), residual_(), scalarProduct_(), gamma_(parms.getGamma()),
330 preSteps_(parms.getNoPreSmoothSteps()), postSteps_(parms.getNoPostSmoothSteps()),
331 buildHierarchy_(true),
332 symmetric(symmetric_), coarsesolverconverged(true),
333 coarseSmoother_(), verbosity_(criterion.debugLevel())
335 if(preSteps_>1||postSteps_>1)
337 std::cerr<<
"WARNING only one step of smoothing is supported!"<<std::endl;
338 preSteps_=postSteps_=1;
340 static_assert(std::is_same<PI,SequentialInformation>::value,
341 "Currently only sequential runs are supported");
345 createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
348 template<
class M,
class X,
class PI,
class A>
351 if(buildHierarchy_) {
355 coarseSmoother_.reset();
368 template<
class M,
class X,
class PI,
class A>
376 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
378 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
379 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.elapsed()<<
" seconds."<<std::endl;
381 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
385 sargs.iterations = 1;
388 cargs.setArgs(sargs);
389 if(matrices_->redistributeInformation().back().isSetup()) {
391 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
392 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
394 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
395 cargs.setComm(*matrices_->parallelInformation().coarsest());
399 scalarProduct_ = createScalarProduct<X>(cargs.getComm(),
category());
401 #if HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK 402 #if HAVE_SUITESPARSE_UMFPACK 403 #define DIRECTSOLVER UMFPack 405 #define DIRECTSOLVER SuperLU 408 if(std::is_same<ParallelInformation,SequentialInformation>::value
409 || matrices_->parallelInformation().coarsest()->communicator().size()==1
410 || (matrices_->parallelInformation().coarsest().isRedistributed()
411 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
412 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
413 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
414 std::cout<<
"Using superlu"<<std::endl;
415 if(matrices_->parallelInformation().coarsest().isRedistributed())
417 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
419 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(),
false,
false));
423 solver_.reset(
new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(),
false,
false));
426 #endif // HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK 428 if(matrices_->parallelInformation().coarsest().isRedistributed())
430 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
432 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
434 *coarseSmoother_, 1E-2, 1000, 0));
438 solver_.reset(
new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
440 *coarseSmoother_, 1E-2, 1000, 0));
444 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
445 std::cout<<
"Building Hierarchy of "<<matrices_->maxlevels()<<
" levels took "<<watch.elapsed()<<
" seconds."<<std::endl;
449 template<
class M,
class X,
class PI,
class A>
457 typedef typename M::matrix_type
Matrix;
464 const Matrix&
mat=matrices_->matrices().finest()->getmat();
465 for(RowIter row=mat.
begin(); row!=mat.
end(); ++row) {
466 bool isDirichlet =
true;
467 bool hasDiagonal =
false;
469 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
470 if(row.index()==
col.index()) {
472 hasDiagonal = (*
col != zero);
478 if(isDirichlet && hasDiagonal)
479 diag->solve(x[row.index()], b[row.index()]);
482 std::cout<<
" Preprocessing Dirichlet took "<<watch1.elapsed()<<std::endl;
485 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
486 Range* copy =
new Range(b);
490 Domain* dcopy =
new Domain(x);
496 matrices_->coarsenVector(*rhs_);
497 matrices_->coarsenVector(*lhs_);
498 matrices_->coarsenVector(*residual_);
505 template<
class M,
class X,
class PI,
class A>
508 return matrices_->levels();
510 template<
class M,
class X,
class PI,
class A>
513 return matrices_->maxlevels();
517 template<
class M,
class X,
class PI,
class A>
520 LevelContext levelContext;
522 initIteratorsWithFineLevel(levelContext);
524 assert(v.two_norm()==0);
527 if(matrices_->maxlevels()==1){
530 mgc(levelContext, v, b);
532 mgc(levelContext, v, d);
533 if(postSteps_==0||matrices_->maxlevels()==1)
534 levelContext.pinfo->copyOwnerToAll(v, v);
537 template<
class M,
class X,
class PI,
class A>
540 levelContext.matrix = matrices_->matrices().finest();
541 levelContext.pinfo = matrices_->parallelInformation().finest();
542 levelContext.redist =
543 matrices_->redistributeInformation().begin();
544 levelContext.aggregates = matrices_->aggregatesMaps().begin();
545 levelContext.lhs = lhs_->
finest();
546 levelContext.residual = residual_->
finest();
547 levelContext.rhs = rhs_->
finest();
548 levelContext.level=0;
551 template<
class M,
class X,
class PI,
class A>
555 bool processNextLevel=
true;
557 if(levelContext.redist->isSetup()) {
559 levelContext.redist->redistribute(static_cast<const Range&>(*levelContext.residual),
560 levelContext.residual.getRedistributed());
561 processNextLevel = levelContext.residual.getRedistributed().size()>0;
562 if(processNextLevel) {
564 ++levelContext.pinfo;
567 static_cast<const Range&>(levelContext.residual.getRedistributed()),
568 *levelContext.pinfo);
573 ++levelContext.pinfo;
576 static_cast<const Range&>(*levelContext.residual), *levelContext.pinfo);
579 if(processNextLevel) {
581 ++levelContext.residual;
583 ++levelContext.matrix;
584 ++levelContext.level;
585 ++levelContext.redist;
587 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
589 ++levelContext.aggregates;
593 *levelContext.residual=0;
595 return processNextLevel;
598 template<
class M,
class X,
class PI,
class A>
602 if(processNextLevel) {
603 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
605 --levelContext.aggregates;
607 --levelContext.redist;
608 --levelContext.level;
610 --levelContext.matrix;
611 --levelContext.residual;
616 if(levelContext.redist->isSetup()) {
621 levelContext.lhs.getRedistributed(),
622 matrices_->getProlongationDampingFactor(),
623 *levelContext.pinfo, *levelContext.redist);
627 matrices_->getProlongationDampingFactor(), *levelContext.pinfo);
633 if(processNextLevel) {
640 template<
class M,
class X,
class PI,
class A>
642 ::presmooth(LevelContext& levelContext, Domain& x,
const Range& b)
646 *levelContext.residual,
650 template<
class M,
class X,
class PI,
class A>
652 ::postsmooth(LevelContext& levelContext, Domain& x,
const Range& b)
655 ::apply(levelContext.matrix->getmat(), x, *levelContext.residual, b);
659 template<
class M,
class X,
class PI,
class A>
665 template<
class M,
class X,
class PI,
class A>
668 if(levelContext.matrix == matrices_->matrices().coarsest() &&
levels()==
maxlevels()) {
672 if(levelContext.redist->isSetup()) {
673 levelContext.redist->redistribute(b, levelContext.rhs.getRedistributed());
674 if(levelContext.rhs.getRedistributed().size()>0) {
676 levelContext.pinfo.getRedistributed().copyOwnerToAll(levelContext.rhs.getRedistributed(),
677 levelContext.rhs.getRedistributed());
678 solver_->apply(levelContext.lhs.getRedistributed(), levelContext.rhs.getRedistributed(), res);
680 levelContext.redist->redistributeBackward(v, levelContext.lhs.getRedistributed());
681 levelContext.pinfo->copyOwnerToAll(v, v);
683 levelContext.pinfo->copyOwnerToAll(b, b);
684 solver_->apply(v, const_cast<Range&>(b), res);
690 coarsesolverconverged =
false;
693 presmooth(levelContext, v, b);
696 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 697 bool processNextLevel = moveToCoarseLevel(levelContext);
699 if(processNextLevel) {
701 for(std::size_t i=0; i<gamma_; i++)
702 mgc(levelContext, *levelContext.lhs, *levelContext.rhs);
705 moveToFineLevel(levelContext, processNextLevel, v);
710 if(levelContext.matrix == matrices_->matrices().finest()) {
711 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
712 if(!coarsesolverconverged)
713 DUNE_THROW(MathError,
"Coarse solver did not converge");
718 postsmooth(levelContext, v, b);
726 template<
class M,
class X,
class PI,
class A>
729 DUNE_UNUSED_PARAMETER(x);
738 template<
class M,
class X,
class PI,
class A>
742 matrices_->getCoarsestAggregatesOnFinest(cont);
Sequential SSOR preconditioner.
Definition: preconditioners.hh:135
void postsmooth(LevelContext &levelContext, size_t steps)
Apply post smoothing on the current level.
Definition: smoother.hh:476
MatrixHierarchy< M, ParallelInformation, A > OperatorHierarchy
The operator hierarchy type.
Definition: fastamg.hh:71
T block_type
Export the type representing the components.
Definition: matrix.hh:562
const void * Arguments
A type holding all the arguments needed to call the constructor.
Definition: construction.hh:44
void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: fastamg.hh:450
void getCoarsestAggregateNumbers(std::vector< std::size_t, A1 > &cont)
Get the aggregate number of each unknown on the coarsest level.
Definition: fastamg.hh:740
static void restrictVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, const Vector &fine, T &comm)
Classes for using SuperLU with ISTL matrices.
Definition: allocator.hh:7
Classes for using UMFPack with ISTL matrices.
Traits class for generically constructing non default constructable types.
Definition: novlpschwarz.hh:247
static void apply(const M &A, X &x, Y &d, const Y &b)
Definition: fastamgsmoother.hh:53
OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix
The iterator over the matrices.
Definition: fastamg.hh:184
Some generic functions for pretty printing vectors and matrices.
LevelIterator< Hierarchy< ParallelInformation, Allocator >, ParallelInformation > Iterator
Type of the mutable iterator.
Definition: hierarchy.hh:250
All parameters for AMG.
Definition: parameters.hh:390
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
X Domain
The domain type.
Definition: fastamg.hh:76
~FastAMG()
Definition: fastamg.hh:349
static void apply(const M &A, X &x, Y &d, const Y &b)
Definition: fastamgsmoother.hh:17
void recalculateHierarchy()
Recalculate the matrix hierarchy.
Definition: fastamg.hh:151
row_type::const_iterator ConstColIterator
Const iterator for the entries of each row.
Definition: matrix.hh:583
Define base class for scalar product and norm.
void post(Domain &x)
Clean up.
Definition: fastamg.hh:727
InverseOperator< X, X > CoarseSolver
the type of the coarse solver.
Definition: fastamg.hh:80
Bi-conjugate Gradient Stabilized (BiCG-STAB)
Definition: solvers.hh:528
Implementations of the inverse operator interface.
OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy
The parallal data distribution hierarchy type.
Definition: fastamg.hh:73
Matrix & mat
Definition: matrixmatrix.hh:345
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
Define general preconditioner interface.
OperatorHierarchy::RedistributeInfoList::const_iterator redist
The iterator over the redistribution information.
Definition: fastamg.hh:192
A generic dynamic dense matrix.
Definition: matrix.hh:554
void presmooth(LevelContext &levelContext, size_t steps)
Apply pre smoothing on the current level.
Definition: smoother.hh:454
T::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:559
Hierarchy< Range, A >::Iterator rhs
The iterator over the right hand sided.
Definition: fastamg.hh:208
static void prolongateVector(const AggregatesMap< Vertex > &aggregates, Vector &coarse, Vector &fine, Vector &fineRedist, T1 damp, R &redistributor=R())
Classes for the generic construction and application of the smoothers.
Hierarchy< Domain, A >::Iterator lhs
The iterator over the left hand side.
Definition: fastamg.hh:200
M Operator
The matrix operator type.
Definition: fastamg.hh:62
A fast (sequential) algebraic multigrid based on agglomeration that saves memory bandwidth.
Definition: fastamg.hh:58
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
Prolongation and restriction for amg.
LevelIterator< const Hierarchy< MatrixOperator, Allocator >, const MatrixOperator > ConstIterator
Type of the const iterator.
Definition: hierarchy.hh:253
Provides a classes representing the hierarchies in AMG.
The hierarchies build by the coarsening process.
Definition: hierarchy.hh:309
A hierarchy of coantainers (e.g. matrices or vectors)
Definition: hierarchy.hh:66
Category
Definition: solvercategory.hh:21
Definition: solvertype.hh:13
Traits class for getting the attribute class of a smoother.
Definition: smoother.hh:63
std::size_t maxlevels()
Definition: fastamg.hh:511
virtual SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: fastamg.hh:124
Col col
Definition: matrixmatrix.hh:349
ConstIterator class for sequential access.
Definition: matrix.hh:397
Base class for scalar product and norm computation.
Definition: scalarproducts.hh:46
void apply(Domain &v, const Range &d)
Apply one step of the preconditioner to the system A(v)=d.
Definition: fastamg.hh:518
OperatorHierarchy::AggregatesMapList::const_iterator aggregates
The iterator over the aggregates maps.
Definition: fastamg.hh:196
PI ParallelInformation
The type of the parallel information. Either OwnerOverlapCommunication or another type describing the...
Definition: fastamg.hh:69
Statistics about the application of an inverse operator.
Definition: solver.hh:40
Category for sequential solvers.
Definition: solvercategory.hh:23
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:30
bool usesDirectCoarseLevelSolver() const
Check whether the coarse solver used is a direct solver.
Definition: fastamg.hh:660
ParallelInformationHierarchy::Iterator pinfo
The iterator over the parallel information.
Definition: fastamg.hh:188
FastAMG(const OperatorHierarchy &matrices, CoarseSolver &coarseSolver, const Parameters &parms, bool symmetric=true)
Construct a new amg with a specific coarse solver.
Definition: fastamg.hh:304
Hierarchy< Domain, A >::Iterator residual
The iterator over the residuals.
Definition: fastamg.hh:204
X Range
The range type.
Definition: fastamg.hh:78
std::size_t level
The level index.
Definition: fastamg.hh:212
Templates characterizing the type of a solver.
Iterator finest()
Get an iterator positioned at the finest level.
Definition: hierarchy.hh:1363
std::size_t levels()
Definition: fastamg.hh:506