198 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator
matrix;
206 typename OperatorHierarchy::RedistributeInfoList::const_iterator
redist;
210 typename OperatorHierarchy::AggregatesMapList::const_iterator
aggregates;
256 void initIteratorsWithFineLevel(LevelContext&
levelContext);
259 std::shared_ptr<OperatorHierarchy> matrices_;
263 std::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
265 std::shared_ptr<CoarseSolver> solver_;
275 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
276 typedef std::shared_ptr<ScalarProduct> ScalarProductPointer;
278 ScalarProductPointer scalarProduct_;
282 std::size_t preSteps_;
284 std::size_t postSteps_;
285 bool buildHierarchy_;
287 bool coarsesolverconverged;
288 std::shared_ptr<Smoother> coarseSmoother_;
290 std::size_t verbosity_;
293 template<
class M,
class X,
class S,
class PI,
class A>
295 : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
296 smoothers_(amg.smoothers_), solver_(amg.solver_),
297 rhs_(), lhs_(), update_(),
298 scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
299 preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
300 buildHierarchy_(amg.buildHierarchy_),
301 additive(amg.additive), coarsesolverconverged(amg.coarsesolverconverged),
302 coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
312 template<
class M,
class X,
class S,
class PI,
class A>
318 rhs_(), lhs_(), update_(), scalarProduct_(0),
319 gamma_(
parms.getGamma()), preSteps_(
parms.getNoPreSmoothSteps()),
320 postSteps_(
parms.getNoPostSmoothSteps()), buildHierarchy_(
false),
321 additive(
parms.getAdditive()), coarsesolverconverged(
true),
322 coarseSmoother_(), verbosity_(
parms.debugLevel())
324 assert(matrices_->isBuilt());
327 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
330 template<
class M,
class X,
class S,
class PI,
class A>
338 rhs_(), lhs_(), update_(), scalarProduct_(),
340 postSteps_(
criterion.getNoPostSmoothSteps()), buildHierarchy_(
true),
341 additive(
criterion.getAdditive()), coarsesolverconverged(
true),
342 coarseSmoother_(), verbosity_(
criterion.debugLevel())
344 static_assert(
static_cast<int>(M::category)==
static_cast<int>(S::category),
345 "Matrix and Solver must match in terms of category!");
353 template<
class M,
class X,
class S,
class PI,
class A>
356 if(buildHierarchy_) {
360 coarseSmoother_.reset();
373 template<
class M,
class X,
class S,
class PI,
class A>
379 matrices_.reset(
new OperatorHierarchy(matrix, pinfo));
384 matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
386 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
388 SmootherArgs
sargs(smootherArgs_);
389 sargs.iterations = 1;
393 if(matrices_->redistributeInformation().back().isSetup()) {
395 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
396 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
398 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
399 cargs.setComm(*matrices_->parallelInformation().coarsest());
403 scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
405#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
406#if HAVE_SUITESPARSE_UMFPACK
407#define DIRECTSOLVER UMFPack
409#define DIRECTSOLVER SuperLU
412 if(std::is_same<ParallelInformation,SequentialInformation>::value
413 || matrices_->parallelInformation().coarsest()->communicator().size()==1
414 || (matrices_->parallelInformation().coarsest().isRedistributed()
415 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
416 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) {
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));
426 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
427 std::cout<<
"Using a direct coarse solver (" <<
static_cast< DIRECTSOLVER<typename M::matrix_type>*
>(solver_.get())->name() <<
")" << std::endl;
432 if(matrices_->parallelInformation().coarsest().isRedistributed())
434 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
436 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(matrices_->matrices().coarsest().getRedistributed()),
438 *coarseSmoother_, 1E-2, 1000, 0));
442 solver_.reset(
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
444 *coarseSmoother_, 1E-2, 1000, 0));
448 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
449 std::cout<<
"Building hierarchy of "<<matrices_->maxlevels()<<
" levels "
450 <<
"(inclusive coarse solver) took "<<watch.elapsed()<<
" seconds."<<std::endl;
454 template<
class M,
class X,
class S,
class PI,
class A>
461 typedef typename M::matrix_type
Matrix;
468 const Matrix&
mat=matrices_->matrices().finest()->getmat();
473 for(ColIter
col=row->begin();
col!=row->end(); ++
col) {
474 if(row.index()==
col.index()) {
483 diagonal.solve(x[row.index()], b[row.index()]);
486 if(smoothers_->levels()>0)
487 smoothers_->finest()->pre(x,b);
490 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
503 matrices_->coarsenVector(*rhs_);
504 matrices_->coarsenVector(*lhs_);
505 matrices_->coarsenVector(*update_);
511 Iterator coarsest = smoothers_->coarsest();
512 Iterator smoother = smoothers_->finest();
515 if(smoothers_->levels()>0) {
517 assert(lhs_->levels()==rhs_->levels());
518 assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
519 assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
521 if(smoother!=coarsest)
522 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
523 smoother->pre(*lhs,*rhs);
524 smoother->pre(*lhs,*rhs);
534 template<
class M,
class X,
class S,
class PI,
class A>
537 return matrices_->levels();
539 template<
class M,
class X,
class S,
class PI,
class A>
542 return matrices_->maxlevels();
546 template<
class M,
class X,
class S,
class PI,
class A>
567 if(postSteps_==0||matrices_->maxlevels()==1)
575 template<
class M,
class X,
class S,
class PI,
class A>
580 levelContext.pinfo = matrices_->parallelInformation().finest();
582 matrices_->redistributeInformation().begin();
583 levelContext.aggregates = matrices_->aggregatesMaps().begin();
589 template<
class M,
class X,
class S,
class PI,
class A>
604 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
606 static_cast<const Range&
>(
fineRhs.getRedistributed()),
611 typename Hierarchy<Range,A>::Iterator fineRhs = levelContext.rhs++;
612 ++levelContext.pinfo;
613 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
614 ::restrictVector(*(*levelContext.aggregates),
615 *levelContext.rhs,
static_cast<const Range&
>(*fineRhs),
616 *levelContext.pinfo);
619 if(processNextLevel) {
622 ++levelContext.update;
623 ++levelContext.matrix;
624 ++levelContext.level;
625 ++levelContext.redist;
627 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
629 ++levelContext.smoother;
630 ++levelContext.aggregates;
633 *levelContext.update=0;
635 return processNextLevel;
638 template<
class M,
class X,
class S,
class PI,
class A>
640 ::moveToFineLevel(LevelContext& levelContext,
bool processNextLevel)
642 if(processNextLevel) {
643 if(levelContext.matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) {
645 --levelContext.smoother;
646 --levelContext.aggregates;
648 --levelContext.redist;
649 --levelContext.level;
651 --levelContext.matrix;
655 --levelContext.pinfo;
657 if(levelContext.redist->isSetup()) {
659 levelContext.lhs.getRedistributed()=0;
660 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
661 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
662 levelContext.lhs.getRedistributed(),
663 matrices_->getProlongationDampingFactor(),
664 *levelContext.pinfo, *levelContext.redist);
667 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
668 ::prolongateVector(*(*levelContext.aggregates), *levelContext.update, *levelContext.lhs,
669 matrices_->getProlongationDampingFactor(),
670 *levelContext.pinfo);
674 if(processNextLevel) {
675 --levelContext.update;
679 *levelContext.update += *levelContext.lhs;
682 template<
class M,
class X,
class S,
class PI,
class A>