2 #ifndef DUNE_AMG_AMG_HH
3 #define DUNE_AMG_AMG_HH
6 #include<dune/common/exceptions.hh>
14 #include<dune/common/typetraits.hh>
15 #include<dune/common/exceptions.hh>
38 template<
class M,
class X,
class S,
class P,
class K,
class A>
51 template<
class M,
class X,
class S,
class PI=SequentialInformation,
52 class A=std::allocator<X> >
55 template<
class M1,
class X1,
class S1,
class P1,
class K1,
class A1>
113 std::size_t preSmoothingSteps,
114 std::size_t postSmoothingSteps,
115 bool additive=
false) DUNE_DEPRECATED;
148 AMG(const
Operator& fineOperator, const C& criterion,
150 std::
size_t preSmoothingSteps,
151 std::
size_t postSmoothingSteps,
167 AMG(const
Operator& fineOperator, const C& criterion,
219 typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
220 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
236 void moveToFineLevel(
bool processedFineLevel);
239 bool moveToCoarseLevel();
242 void initIteratorsWithFineLevel();
261 typedef typename ScalarProductChooser::ScalarProduct
ScalarProduct;
267 std::size_t preSteps_;
269 std::size_t postSteps_;
271 bool buildHierarchy_;
273 bool coarsesolverconverged;
276 std::size_t verbosity_;
279 template<
class M,
class X,
class S,
class PI,
class A>
282 std::size_t gamma, std::size_t preSmoothingSteps,
283 std::size_t postSmoothingSteps,
bool additive_)
284 : matrices_(&matrices), smootherArgs_(smootherArgs),
285 smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
286 gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
287 additive(additive_), coarsesolverconverged(true),
288 coarseSmoother_(), verbosity_(2)
296 template<
class M,
class X,
class S,
class PI,
class A>
300 : matrices_(&matrices), smootherArgs_(smootherArgs),
301 smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
302 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
303 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
304 additive(parms.getAdditive()), coarsesolverconverged(true),
305 coarseSmoother_(), verbosity_(parms.debugLevel())
313 template<
class M,
class X,
class S,
class PI,
class A>
318 std::size_t gamma, std::size_t preSmoothingSteps,
319 std::size_t postSmoothingSteps,
322 : smootherArgs_(smootherArgs),
323 smoothers_(), solver_(), scalarProduct_(0), gamma_(gamma),
324 preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
325 additive(additive_), coarsesolverconverged(true),
326 coarseSmoother_(), verbosity_(criterion.debugLevel())
328 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
329 "Matrix and Solver must match in terms of category!");
336 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
342 std::cout<<
"Building Hierarchy of "<<matrices_->
maxlevels()<<
" levels took "<<watch.elapsed()<<
" seconds."<<std::endl;
345 template<
class M,
class X,
class S,
class PI,
class A>
351 : smootherArgs_(smootherArgs),
352 smoothers_(), solver_(), scalarProduct_(0),
353 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
354 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
355 additive(criterion.getAdditive()), coarsesolverconverged(true),
356 coarseSmoother_(), verbosity_(criterion.debugLevel())
358 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category),
359 "Matrix and Solver must match in terms of category!");
366 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
372 std::cout<<
"Building Hierarchy of "<<matrices_->
maxlevels()<<
" levels took "<<watch.elapsed()<<
" seconds."<<std::endl;
375 template<
class M,
class X,
class S,
class PI,
class A>
382 delete scalarProduct_;
386 template<
class M,
class X,
class S,
class PI,
class A>
389 if(smoothers_.levels()>0)
390 smoothers_.finest()->pre(x,b);
393 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
400 matrices_->coarsenVector(*rhs_);
401 matrices_->coarsenVector(*lhs_);
402 matrices_->coarsenVector(*update_);
408 Iterator coarsest = smoothers_.
coarsest();
409 Iterator smoother = smoothers_.finest();
410 RIterator rhs = rhs_->finest();
411 DIterator lhs = lhs_->finest();
412 if(smoothers_.levels()>0){
414 assert(lhs_->levels()==rhs_->levels());
415 assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
416 assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
418 if(smoother!=coarsest)
419 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
420 smoother->pre(*lhs,*rhs);
421 smoother->pre(*lhs,*rhs);
430 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){
433 sargs.iterations = 1;
436 cargs.setArgs(sargs);
437 if(matrices_->redistributeInformation().back().isSetup()){
439 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
440 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
443 scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
445 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
446 cargs.setComm(*matrices_->parallelInformation().coarsest());
449 scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
453 if(is_same<ParallelInformation,SequentialInformation>::value
454 || matrices_->parallelInformation().coarsest()->communicator().size()==1
455 || (matrices_->parallelInformation().coarsest().isRedistributed()
456 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
457 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){
458 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
459 std::cout<<
"Using superlu"<<std::endl;
460 if(matrices_->parallelInformation().coarsest().isRedistributed())
462 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
472 if(matrices_->parallelInformation().coarsest().isRedistributed())
474 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
476 solver_ =
new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
478 *coarseSmoother_, 1E-2, 10000, 0);
482 solver_ =
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
484 *coarseSmoother_, 1E-2, 1000, 0);
488 template<
class M,
class X,
class S,
class PI,
class A>
491 return matrices_->levels();
493 template<
class M,
class X,
class S,
class PI,
class A>
496 return matrices_->maxlevels();
500 template<
class M,
class X,
class S,
class PI,
class A>
509 initIteratorsWithFineLevel();
519 if(postSteps_==0||matrices_->maxlevels()==1)
520 pinfo->copyOwnerToAll(*update, *update);
527 template<
class M,
class X,
class S,
class PI,
class A>
530 smoother = smoothers_.finest();
531 matrix = matrices_->matrices().finest();
532 pinfo = matrices_->parallelInformation().finest();
534 matrices_->redistributeInformation().begin();
535 aggregates = matrices_->aggregatesMaps().begin();
536 lhs = lhs_->finest();
537 update = update_->finest();
538 rhs = rhs_->finest();
541 template<
class M,
class X,
class S,
class PI,
class A>
543 ::moveToCoarseLevel()
546 bool processNextLevel=
true;
548 if(redist->isSetup()){
549 redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed());
550 processNextLevel =rhs.getRedistributed().size()>0;
551 if(processNextLevel){
556 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo);
563 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
566 if(processNextLevel){
574 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
582 return processNextLevel;
585 template<
class M,
class X,
class S,
class PI,
class A>
587 ::moveToFineLevel(
bool processNextLevel)
589 if(processNextLevel){
590 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
604 if(redist->isSetup()){
606 lhs.getRedistributed()=0;
608 ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(),
614 matrices_->getProlongationDampingFactor(), *pinfo);
618 if(processNextLevel){
627 template<
class M,
class X,
class S,
class PI,
class A>
632 for(std::size_t i=0; i < preSteps_; ++i){
639 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
640 pinfo->project(*rhs);
644 template<
class M,
class X,
class S,
class PI,
class A>
649 for(std::size_t i=0; i < postSteps_; ++i){
651 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
653 pinfo->project(*rhs);
661 template<
class M,
class X,
class S,
class PI,
class A>
667 template<
class M,
class X,
class S,
class PI,
class A>
669 if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){
673 if(redist->isSetup()){
674 redist->redistribute(*rhs, rhs.getRedistributed());
675 if(rhs.getRedistributed().size()>0){
677 pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed());
678 solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res);
680 redist->redistributeBackward(*update, update.getRedistributed());
681 pinfo->copyOwnerToAll(*update, *update);
683 pinfo->copyOwnerToAll(*rhs, *rhs);
684 solver_->apply(*update, *rhs, res);
688 coarsesolverconverged =
false;
693 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
694 bool processNextLevel = moveToCoarseLevel();
696 if(processNextLevel){
698 for(std::size_t i=0; i<gamma_; i++)
702 moveToFineLevel(processNextLevel);
707 if(matrix == matrices_->matrices().finest()){
708 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
709 if(!coarsesolverconverged)
710 DUNE_THROW(MathError,
"Coarse solver did not converge");
718 template<
class M,
class X,
class S,
class PI,
class A>
719 void AMG<M,X,S,PI,A>::additiveMgc(){
722 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
725 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
730 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
736 lhs = lhs_->finest();
739 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){
742 smoother->apply(*lhs, *rhs);
746 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
747 InverseOperatorResult res;
748 pinfo->copyOwnerToAll(*rhs, *rhs);
749 solver_->apply(*lhs, *rhs, res);
752 DUNE_THROW(MathError,
"Coarse solver did not converge");
762 ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo);
768 template<
class M,
class X,
class S,
class PI,
class A>
782 Iterator coarsest = smoothers_.
coarsest();
783 Iterator smoother = smoothers_.finest();
784 DIterator lhs = lhs_->finest();
785 if(smoothers_.levels()>0){
786 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
787 smoother->post(*lhs);
788 if(smoother!=coarsest)
789 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
790 smoother->post(*lhs);
791 smoother->post(*lhs);
794 delete &(*lhs_->finest());
796 delete &(*update_->finest());
798 delete &(*rhs_->finest());
802 template<
class M,
class X,
class S,
class PI,
class A>
806 matrices_->getCoarsestAggregatesOnFinest(cont);