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>
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>
393 typedef typename M::matrix_type
Matrix;
400 const Matrix&
mat=matrices_->matrices().finest()->getmat();
402 bool isDirichlet =
true;
403 bool hasDiagonal =
false;
406 if(
row.index()==
col.index()){
414 if(isDirichlet && hasDiagonal)
415 diagonal.solve(x[
row.index()], b[
row.index()]);
418 if(smoothers_.levels()>0)
419 smoothers_.finest()->pre(x,b);
422 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
429 matrices_->coarsenVector(*rhs_);
430 matrices_->coarsenVector(*lhs_);
431 matrices_->coarsenVector(*update_);
437 Iterator coarsest = smoothers_.
coarsest();
438 Iterator smoother = smoothers_.finest();
439 RIterator rhs = rhs_->finest();
440 DIterator lhs = lhs_->finest();
441 if(smoothers_.levels()>0){
443 assert(lhs_->levels()==rhs_->levels());
444 assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
445 assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
447 if(smoother!=coarsest)
448 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
449 smoother->pre(*lhs,*rhs);
450 smoother->pre(*lhs,*rhs);
459 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){
462 sargs.iterations = 1;
465 cargs.setArgs(sargs);
466 if(matrices_->redistributeInformation().back().isSetup()){
468 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
469 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
472 scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
474 cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
475 cargs.setComm(*matrices_->parallelInformation().coarsest());
478 scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
482 if(is_same<ParallelInformation,SequentialInformation>::value
483 || matrices_->parallelInformation().coarsest()->communicator().size()==1
484 || (matrices_->parallelInformation().coarsest().isRedistributed()
485 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
486 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){
487 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
488 std::cout<<
"Using superlu"<<std::endl;
489 if(matrices_->parallelInformation().coarsest().isRedistributed())
491 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
501 if(matrices_->parallelInformation().coarsest().isRedistributed())
503 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
505 solver_ =
new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
507 *coarseSmoother_, 1E-2, 10000, 0);
511 solver_ =
new BiCGSTABSolver<X>(
const_cast<M&
>(*matrices_->matrices().coarsest()),
513 *coarseSmoother_, 1E-2, 1000, 0);
517 template<
class M,
class X,
class S,
class PI,
class A>
520 return matrices_->levels();
522 template<
class M,
class X,
class S,
class PI,
class A>
525 return matrices_->maxlevels();
529 template<
class M,
class X,
class S,
class PI,
class A>
538 initIteratorsWithFineLevel();
548 if(postSteps_==0||matrices_->maxlevels()==1)
549 pinfo->copyOwnerToAll(*update, *update);
556 template<
class M,
class X,
class S,
class PI,
class A>
559 smoother = smoothers_.finest();
560 matrix = matrices_->matrices().finest();
561 pinfo = matrices_->parallelInformation().finest();
563 matrices_->redistributeInformation().begin();
564 aggregates = matrices_->aggregatesMaps().begin();
565 lhs = lhs_->finest();
566 update = update_->finest();
567 rhs = rhs_->finest();
570 template<
class M,
class X,
class S,
class PI,
class A>
572 ::moveToCoarseLevel()
575 bool processNextLevel=
true;
577 if(redist->isSetup()){
578 redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed());
579 processNextLevel =rhs.getRedistributed().size()>0;
580 if(processNextLevel){
585 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo);
592 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
595 if(processNextLevel){
603 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
611 return processNextLevel;
614 template<
class M,
class X,
class S,
class PI,
class A>
616 ::moveToFineLevel(
bool processNextLevel)
618 if(processNextLevel){
619 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){
633 if(redist->isSetup()){
635 lhs.getRedistributed()=0;
637 ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(),
643 matrices_->getProlongationDampingFactor(), *pinfo);
647 if(processNextLevel){
656 template<
class M,
class X,
class S,
class PI,
class A>
661 for(std::size_t i=0; i < preSteps_; ++i){
668 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
669 pinfo->project(*rhs);
673 template<
class M,
class X,
class S,
class PI,
class A>
678 for(std::size_t i=0; i < postSteps_; ++i){
680 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs);
682 pinfo->project(*rhs);
690 template<
class M,
class X,
class S,
class PI,
class A>
696 template<
class M,
class X,
class S,
class PI,
class A>
698 if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){
702 if(redist->isSetup()){
703 redist->redistribute(*rhs, rhs.getRedistributed());
704 if(rhs.getRedistributed().size()>0){
706 pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed());
707 solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res);
709 redist->redistributeBackward(*update, update.getRedistributed());
710 pinfo->copyOwnerToAll(*update, *update);
712 pinfo->copyOwnerToAll(*rhs, *rhs);
713 solver_->apply(*update, *rhs, res);
717 coarsesolverconverged =
false;
722 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
723 bool processNextLevel = moveToCoarseLevel();
725 if(processNextLevel){
727 for(std::size_t i=0; i<gamma_; i++)
731 moveToFineLevel(processNextLevel);
736 if(matrix == matrices_->matrices().finest()){
737 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
738 if(!coarsesolverconverged)
739 DUNE_THROW(MathError,
"Coarse solver did not converge");
747 template<
class M,
class X,
class S,
class PI,
class A>
748 void AMG<M,X,S,PI,A>::additiveMgc(){
751 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest();
754 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin();
759 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo);
765 lhs = lhs_->finest();
768 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){
771 smoother->apply(*lhs, *rhs);
775 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION
776 InverseOperatorResult res;
777 pinfo->copyOwnerToAll(*rhs, *rhs);
778 solver_->apply(*lhs, *rhs, res);
781 DUNE_THROW(MathError,
"Coarse solver did not converge");
791 ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo);
797 template<
class M,
class X,
class S,
class PI,
class A>
811 Iterator coarsest = smoothers_.
coarsest();
812 Iterator smoother = smoothers_.finest();
813 DIterator lhs = lhs_->finest();
814 if(smoothers_.levels()>0){
815 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
816 smoother->post(*lhs);
817 if(smoother!=coarsest)
818 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs)
819 smoother->post(*lhs);
820 smoother->post(*lhs);
823 delete &(*lhs_->finest());
825 delete &(*update_->finest());
827 delete &(*rhs_->finest());
831 template<
class M,
class X,
class S,
class PI,
class A>
835 matrices_->getCoarsestAggregatesOnFinest(cont);