dune-istl
2.2.0
|
00001 // $Id: amg.hh 1547 2012-04-19 17:19:40Z mblatt $ 00002 #ifndef DUNE_AMG_AMG_HH 00003 #define DUNE_AMG_AMG_HH 00004 00005 #include<memory> 00006 #include<dune/common/exceptions.hh> 00007 #include<dune/istl/paamg/smoother.hh> 00008 #include<dune/istl/paamg/transfer.hh> 00009 #include<dune/istl/paamg/hierarchy.hh> 00010 #include<dune/istl/solvers.hh> 00011 #include<dune/istl/scalarproducts.hh> 00012 #include<dune/istl/superlu.hh> 00013 #include<dune/istl/solvertype.hh> 00014 #include<dune/common/typetraits.hh> 00015 #include<dune/common/exceptions.hh> 00016 00017 namespace Dune 00018 { 00019 namespace Amg 00020 { 00038 template<class M, class X, class S, class P, class K, class A> 00039 class KAMG; 00040 00041 template<class T> 00042 class KAmgTwoGrid; 00043 00051 template<class M, class X, class S, class PI=SequentialInformation, 00052 class A=std::allocator<X> > 00053 class AMG : public Preconditioner<X,X> 00054 { 00055 template<class M1, class X1, class S1, class P1, class K1, class A1> 00056 friend class KAMG; 00057 00058 friend class KAmgTwoGrid<AMG>; 00059 00060 public: 00062 typedef M Operator; 00069 typedef PI ParallelInformation; 00071 typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy; 00073 typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy; 00074 00076 typedef X Domain; 00078 typedef X Range; 00080 typedef InverseOperator<X,X> CoarseSolver; 00086 typedef S Smoother; 00087 00089 typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs; 00090 00091 enum { 00093 category = S::category 00094 }; 00095 00111 AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00112 const SmootherArgs& smootherArgs, std::size_t gamma, 00113 std::size_t preSmoothingSteps, 00114 std::size_t postSmoothingSteps, 00115 bool additive=false) DUNE_DEPRECATED; 00116 00126 AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00127 const SmootherArgs& smootherArgs, const Parameters& parms); 00128 00147 template<class C> 00148 AMG(const Operator& fineOperator, const C& criterion, 00149 const SmootherArgs& smootherArgs, std::size_t gamma, 00150 std::size_t preSmoothingSteps, 00151 std::size_t postSmoothingSteps, 00152 bool additive=false, 00153 const ParallelInformation& pinfo=ParallelInformation()) DUNE_DEPRECATED; 00154 00166 template<class C> 00167 AMG(const Operator& fineOperator, const C& criterion, 00168 const SmootherArgs& smootherArgs, 00169 const ParallelInformation& pinfo=ParallelInformation()); 00170 00171 ~AMG(); 00172 00174 void pre(Domain& x, Range& b); 00175 00177 void apply(Domain& v, const Range& d); 00178 00180 void post(Domain& x); 00181 00186 template<class A1> 00187 void getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont); 00188 00189 std::size_t levels(); 00190 00191 std::size_t maxlevels(); 00192 00201 void recalculateHierarchy() 00202 { 00203 matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>()); 00204 } 00205 00210 bool usesDirectCoarseLevelSolver() const; 00211 00212 private: 00214 void mgc(); 00215 00216 typename Hierarchy<Smoother,A>::Iterator smoother; 00217 typename OperatorHierarchy::ParallelMatrixHierarchy::ConstIterator matrix; 00218 typename ParallelInformationHierarchy::Iterator pinfo; 00219 typename OperatorHierarchy::RedistributeInfoList::const_iterator redist; 00220 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates; 00221 typename Hierarchy<Domain,A>::Iterator lhs; 00222 typename Hierarchy<Domain,A>::Iterator update; 00223 typename Hierarchy<Range,A>::Iterator rhs; 00224 00225 void additiveMgc(); 00226 00228 void presmooth(); 00229 00231 void postsmooth(); 00232 00236 void moveToFineLevel(bool processedFineLevel); 00237 00239 bool moveToCoarseLevel(); 00240 00242 void initIteratorsWithFineLevel(); 00243 00245 OperatorHierarchy* matrices_; 00247 SmootherArgs smootherArgs_; 00249 Hierarchy<Smoother,A> smoothers_; 00251 CoarseSolver* solver_; 00253 Hierarchy<Range,A>* rhs_; 00255 Hierarchy<Domain,A>* lhs_; 00257 Hierarchy<Domain,A>* update_; 00259 typedef Dune::ScalarProductChooser<X,PI,M::category> ScalarProductChooser; 00261 typedef typename ScalarProductChooser::ScalarProduct ScalarProduct; 00263 ScalarProduct* scalarProduct_; 00265 std::size_t gamma_; 00267 std::size_t preSteps_; 00269 std::size_t postSteps_; 00270 std::size_t level; 00271 bool buildHierarchy_; 00272 bool additive; 00273 bool coarsesolverconverged; 00274 Smoother *coarseSmoother_; 00276 std::size_t verbosity_; 00277 }; 00278 00279 template<class M, class X, class S, class PI, class A> 00280 AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00281 const SmootherArgs& smootherArgs, 00282 std::size_t gamma, std::size_t preSmoothingSteps, 00283 std::size_t postSmoothingSteps, bool additive_) 00284 : matrices_(&matrices), smootherArgs_(smootherArgs), 00285 smoothers_(), solver_(&coarseSolver), scalarProduct_(0), 00286 gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false), 00287 additive(additive_), coarsesolverconverged(true), 00288 coarseSmoother_(), verbosity_(2) 00289 { 00290 assert(matrices_->isBuilt()); 00291 00292 // build the necessary smoother hierarchies 00293 matrices_->coarsenSmoother(smoothers_, smootherArgs_); 00294 } 00295 00296 template<class M, class X, class S, class PI, class A> 00297 AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, 00298 const SmootherArgs& smootherArgs, 00299 const Parameters& parms) 00300 : matrices_(&matrices), smootherArgs_(smootherArgs), 00301 smoothers_(), solver_(&coarseSolver), scalarProduct_(0), 00302 gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()), 00303 postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false), 00304 additive(parms.getAdditive()), coarsesolverconverged(true), 00305 coarseSmoother_(), verbosity_(parms.debugLevel()) 00306 { 00307 assert(matrices_->isBuilt()); 00308 00309 // build the necessary smoother hierarchies 00310 matrices_->coarsenSmoother(smoothers_, smootherArgs_); 00311 } 00312 00313 template<class M, class X, class S, class PI, class A> 00314 template<class C> 00315 AMG<M,X,S,PI,A>::AMG(const Operator& matrix, 00316 const C& criterion, 00317 const SmootherArgs& smootherArgs, 00318 std::size_t gamma, std::size_t preSmoothingSteps, 00319 std::size_t postSmoothingSteps, 00320 bool additive_, 00321 const PI& pinfo) 00322 : smootherArgs_(smootherArgs), 00323 smoothers_(), solver_(), scalarProduct_(0), gamma_(gamma), 00324 preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true), 00325 additive(additive_), coarsesolverconverged(true), 00326 coarseSmoother_(), verbosity_(criterion.debugLevel()) 00327 { 00328 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category), 00329 "Matrix and Solver must match in terms of category!"); 00330 // TODO: reestablish compile time checks. 00331 //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category), 00332 // "Matrix and Solver must match in terms of category!"); 00333 Timer watch; 00334 matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo); 00335 00336 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion); 00337 00338 // build the necessary smoother hierarchies 00339 matrices_->coarsenSmoother(smoothers_, smootherArgs_); 00340 00341 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) 00342 std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl; 00343 } 00344 00345 template<class M, class X, class S, class PI, class A> 00346 template<class C> 00347 AMG<M,X,S,PI,A>::AMG(const Operator& matrix, 00348 const C& criterion, 00349 const SmootherArgs& smootherArgs, 00350 const PI& pinfo) 00351 : smootherArgs_(smootherArgs), 00352 smoothers_(), solver_(), scalarProduct_(0), 00353 gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()), 00354 postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true), 00355 additive(criterion.getAdditive()), coarsesolverconverged(true), 00356 coarseSmoother_(), verbosity_(criterion.debugLevel()) 00357 { 00358 dune_static_assert(static_cast<int>(M::category)==static_cast<int>(S::category), 00359 "Matrix and Solver must match in terms of category!"); 00360 // TODO: reestablish compile time checks. 00361 //dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category), 00362 // "Matrix and Solver must match in terms of category!"); 00363 Timer watch; 00364 matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo); 00365 00366 matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion); 00367 00368 // build the necessary smoother hierarchies 00369 matrices_->coarsenSmoother(smoothers_, smootherArgs_); 00370 00371 if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0) 00372 std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl; 00373 } 00374 00375 template<class M, class X, class S, class PI, class A> 00376 AMG<M,X,S,PI,A>::~AMG() 00377 { 00378 if(buildHierarchy_){ 00379 delete matrices_; 00380 } 00381 if(scalarProduct_) 00382 delete scalarProduct_; 00383 } 00384 00385 00386 template<class M, class X, class S, class PI, class A> 00387 void AMG<M,X,S,PI,A>::pre(Domain& x, Range& b) 00388 { 00389 if(smoothers_.levels()>0) 00390 smoothers_.finest()->pre(x,b); 00391 else 00392 // No smoother to make x consistent! Do it by hand 00393 matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x); 00394 Range* copy = new Range(b); 00395 rhs_ = new Hierarchy<Range,A>(*copy); 00396 Domain* dcopy = new Domain(x); 00397 lhs_ = new Hierarchy<Domain,A>(*dcopy); 00398 dcopy = new Domain(x); 00399 update_ = new Hierarchy<Domain,A>(*dcopy); 00400 matrices_->coarsenVector(*rhs_); 00401 matrices_->coarsenVector(*lhs_); 00402 matrices_->coarsenVector(*update_); 00403 00404 // Preprocess all smoothers 00405 typedef typename Hierarchy<Smoother,A>::Iterator Iterator; 00406 typedef typename Hierarchy<Range,A>::Iterator RIterator; 00407 typedef typename Hierarchy<Domain,A>::Iterator DIterator; 00408 Iterator coarsest = smoothers_.coarsest(); 00409 Iterator smoother = smoothers_.finest(); 00410 RIterator rhs = rhs_->finest(); 00411 DIterator lhs = lhs_->finest(); 00412 if(smoothers_.levels()>0){ 00413 00414 assert(lhs_->levels()==rhs_->levels()); 00415 assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels()); 00416 assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels()); 00417 00418 if(smoother!=coarsest) 00419 for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs) 00420 smoother->pre(*lhs,*rhs); 00421 smoother->pre(*lhs,*rhs); 00422 } 00423 00424 00425 // The preconditioner might change x and b. So we have to 00426 // copy the changes to the original vectors. 00427 x = *lhs_->finest(); 00428 b = *rhs_->finest(); 00429 00430 if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()){ 00431 // We have the carsest level. Create the coarse Solver 00432 SmootherArgs sargs(smootherArgs_); 00433 sargs.iterations = 1; 00434 00435 typename ConstructionTraits<Smoother>::Arguments cargs; 00436 cargs.setArgs(sargs); 00437 if(matrices_->redistributeInformation().back().isSetup()){ 00438 // Solve on the redistributed partitioning 00439 cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat()); 00440 cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed()); 00441 00442 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs); 00443 scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed()); 00444 }else{ 00445 cargs.setMatrix(matrices_->matrices().coarsest()->getmat()); 00446 cargs.setComm(*matrices_->parallelInformation().coarsest()); 00447 00448 coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs); 00449 scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest()); 00450 } 00451 #if HAVE_SUPERLU 00452 // Use superlu if we are purely sequential or with only one processor on the coarsest level. 00453 if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode 00454 || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor 00455 || (matrices_->parallelInformation().coarsest().isRedistributed() 00456 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1 00457 && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)){ // redistribute and 1 proc 00458 if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0) 00459 std::cout<<"Using superlu"<<std::endl; 00460 if(matrices_->parallelInformation().coarsest().isRedistributed()) 00461 { 00462 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0) 00463 // We are still participating on this level 00464 solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat()); 00465 else 00466 solver_ = 0; 00467 }else 00468 solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat()); 00469 }else 00470 #endif 00471 { 00472 if(matrices_->parallelInformation().coarsest().isRedistributed()) 00473 { 00474 if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0) 00475 // We are still participating on this level 00476 solver_ = new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()), 00477 *scalarProduct_, 00478 *coarseSmoother_, 1E-2, 10000, 0); 00479 else 00480 solver_ = 0; 00481 }else 00482 solver_ = new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()), 00483 *scalarProduct_, 00484 *coarseSmoother_, 1E-2, 1000, 0); 00485 } 00486 } 00487 } 00488 template<class M, class X, class S, class PI, class A> 00489 std::size_t AMG<M,X,S,PI,A>::levels() 00490 { 00491 return matrices_->levels(); 00492 } 00493 template<class M, class X, class S, class PI, class A> 00494 std::size_t AMG<M,X,S,PI,A>::maxlevels() 00495 { 00496 return matrices_->maxlevels(); 00497 } 00498 00500 template<class M, class X, class S, class PI, class A> 00501 void AMG<M,X,S,PI,A>::apply(Domain& v, const Range& d) 00502 { 00503 if(additive){ 00504 *(rhs_->finest())=d; 00505 additiveMgc(); 00506 v=*lhs_->finest(); 00507 }else{ 00508 // Init all iterators for the current level 00509 initIteratorsWithFineLevel(); 00510 00511 00512 *lhs = v; 00513 *rhs = d; 00514 *update=0; 00515 level=0; 00516 00517 mgc(); 00518 00519 if(postSteps_==0||matrices_->maxlevels()==1) 00520 pinfo->copyOwnerToAll(*update, *update); 00521 00522 v=*update; 00523 } 00524 00525 } 00526 00527 template<class M, class X, class S, class PI, class A> 00528 void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel() 00529 { 00530 smoother = smoothers_.finest(); 00531 matrix = matrices_->matrices().finest(); 00532 pinfo = matrices_->parallelInformation().finest(); 00533 redist = 00534 matrices_->redistributeInformation().begin(); 00535 aggregates = matrices_->aggregatesMaps().begin(); 00536 lhs = lhs_->finest(); 00537 update = update_->finest(); 00538 rhs = rhs_->finest(); 00539 } 00540 00541 template<class M, class X, class S, class PI, class A> 00542 bool AMG<M,X,S,PI,A> 00543 ::moveToCoarseLevel() 00544 { 00545 00546 bool processNextLevel=true; 00547 00548 if(redist->isSetup()){ 00549 redist->redistribute(static_cast<const Range&>(*rhs), rhs.getRedistributed()); 00550 processNextLevel =rhs.getRedistributed().size()>0; 00551 if(processNextLevel){ 00552 //restrict defect to coarse level right hand side. 00553 typename Hierarchy<Range,A>::Iterator fineRhs = rhs++; 00554 ++pinfo; 00555 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00556 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(fineRhs.getRedistributed()), *pinfo); 00557 } 00558 }else{ 00559 //restrict defect to coarse level right hand side. 00560 typename Hierarchy<Range,A>::Iterator fineRhs = rhs++; 00561 ++pinfo; 00562 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00563 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo); 00564 } 00565 00566 if(processNextLevel){ 00567 // prepare coarse system 00568 ++lhs; 00569 ++update; 00570 ++matrix; 00571 ++level; 00572 ++redist; 00573 00574 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){ 00575 // next level is not the globally coarsest one 00576 ++smoother; 00577 ++aggregates; 00578 } 00579 // prepare the update on the next level 00580 *update=0; 00581 } 00582 return processNextLevel; 00583 } 00584 00585 template<class M, class X, class S, class PI, class A> 00586 void AMG<M,X,S,PI,A> 00587 ::moveToFineLevel(bool processNextLevel) 00588 { 00589 if(processNextLevel){ 00590 if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()){ 00591 // previous level is not the globally coarsest one 00592 --smoother; 00593 --aggregates; 00594 } 00595 --redist; 00596 --level; 00597 //prolongate and add the correction (update is in coarse left hand side) 00598 --matrix; 00599 00600 //typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; 00601 --lhs; 00602 --pinfo; 00603 } 00604 if(redist->isSetup()){ 00605 // Need to redistribute during prolongate 00606 lhs.getRedistributed()=0; 00607 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00608 ::prolongate(*(*aggregates), *update, *lhs, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(), 00609 *pinfo, *redist); 00610 }else{ 00611 *lhs=0; 00612 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00613 ::prolongate(*(*aggregates), *update, *lhs, 00614 matrices_->getProlongationDampingFactor(), *pinfo); 00615 } 00616 00617 00618 if(processNextLevel){ 00619 --update; 00620 --rhs; 00621 } 00622 00623 *update += *lhs; 00624 } 00625 00626 00627 template<class M, class X, class S, class PI, class A> 00628 void AMG<M,X,S,PI,A> 00629 ::presmooth() 00630 { 00631 00632 for(std::size_t i=0; i < preSteps_; ++i){ 00633 *lhs=0; 00634 SmootherApplier<S>::preSmooth(*smoother, *lhs, *rhs); 00635 // Accumulate update 00636 *update += *lhs; 00637 00638 // update defect 00639 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs); 00640 pinfo->project(*rhs); 00641 } 00642 } 00643 00644 template<class M, class X, class S, class PI, class A> 00645 void AMG<M,X,S,PI,A> 00646 ::postsmooth() 00647 { 00648 00649 for(std::size_t i=0; i < postSteps_; ++i){ 00650 // update defect 00651 matrix->applyscaleadd(-1,static_cast<const Domain&>(*lhs), *rhs); 00652 *lhs=0; 00653 pinfo->project(*rhs); 00654 SmootherApplier<S>::postSmooth(*smoother, *lhs, *rhs); 00655 // Accumulate update 00656 *update += *lhs; 00657 } 00658 } 00659 00660 00661 template<class M, class X, class S, class PI, class A> 00662 bool AMG<M,X,S,PI,A>::usesDirectCoarseLevelSolver() const 00663 { 00664 return IsDirectSolver< CoarseSolver>::value; 00665 } 00666 00667 template<class M, class X, class S, class PI, class A> 00668 void AMG<M,X,S,PI,A>::mgc(){ 00669 if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()){ 00670 // Solve directly 00671 InverseOperatorResult res; 00672 res.converged=true; // If we do not compute this flag will not get updated 00673 if(redist->isSetup()){ 00674 redist->redistribute(*rhs, rhs.getRedistributed()); 00675 if(rhs.getRedistributed().size()>0){ 00676 // We are still participating in the computation 00677 pinfo.getRedistributed().copyOwnerToAll(rhs.getRedistributed(), rhs.getRedistributed()); 00678 solver_->apply(update.getRedistributed(), rhs.getRedistributed(), res); 00679 } 00680 redist->redistributeBackward(*update, update.getRedistributed()); 00681 pinfo->copyOwnerToAll(*update, *update); 00682 }else{ 00683 pinfo->copyOwnerToAll(*rhs, *rhs); 00684 solver_->apply(*update, *rhs, res); 00685 } 00686 00687 if (!res.converged) 00688 coarsesolverconverged = false; 00689 }else{ 00690 // presmoothing 00691 presmooth(); 00692 00693 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 00694 bool processNextLevel = moveToCoarseLevel(); 00695 00696 if(processNextLevel){ 00697 // next level 00698 for(std::size_t i=0; i<gamma_; i++) 00699 mgc(); 00700 } 00701 00702 moveToFineLevel(processNextLevel); 00703 #else 00704 *lhs=0; 00705 #endif 00706 00707 if(matrix == matrices_->matrices().finest()){ 00708 coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged); 00709 if(!coarsesolverconverged) 00710 DUNE_THROW(MathError, "Coarse solver did not converge"); 00711 } 00712 // postsmoothing 00713 postsmooth(); 00714 00715 } 00716 } 00717 00718 template<class M, class X, class S, class PI, class A> 00719 void AMG<M,X,S,PI,A>::additiveMgc(){ 00720 00721 // restrict residual to all levels 00722 typename ParallelInformationHierarchy::Iterator pinfo=matrices_->parallelInformation().finest(); 00723 typename Hierarchy<Range,A>::Iterator rhs=rhs_->finest(); 00724 typename Hierarchy<Domain,A>::Iterator lhs = lhs_->finest(); 00725 typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates=matrices_->aggregatesMaps().begin(); 00726 00727 for(typename Hierarchy<Range,A>::Iterator fineRhs=rhs++; fineRhs != rhs_->coarsest(); fineRhs=rhs++, ++aggregates){ 00728 ++pinfo; 00729 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00730 ::restrict(*(*aggregates), *rhs, static_cast<const Range&>(*fineRhs), *pinfo); 00731 } 00732 00733 // pinfo is invalid, set to coarsest level 00734 //pinfo = matrices_->parallelInformation().coarsest 00735 // calculate correction for all levels 00736 lhs = lhs_->finest(); 00737 typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest(); 00738 00739 for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother){ 00740 // presmoothing 00741 *lhs=0; 00742 smoother->apply(*lhs, *rhs); 00743 } 00744 00745 // Coarse level solve 00746 #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION 00747 InverseOperatorResult res; 00748 pinfo->copyOwnerToAll(*rhs, *rhs); 00749 solver_->apply(*lhs, *rhs, res); 00750 00751 if(!res.converged) 00752 DUNE_THROW(MathError, "Coarse solver did not converge"); 00753 #else 00754 *lhs=0; 00755 #endif 00756 // Prologate and add up corrections from all levels 00757 --pinfo; 00758 --aggregates; 00759 00760 for(typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--; coarseLhs != lhs_->finest(); coarseLhs = lhs--, --aggregates, --pinfo){ 00761 Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation> 00762 ::prolongate(*(*aggregates), *coarseLhs, *lhs, 1, *pinfo); 00763 } 00764 } 00765 00766 00768 template<class M, class X, class S, class PI, class A> 00769 void AMG<M,X,S,PI,A>::post(Domain& x) 00770 { 00771 if(buildHierarchy_){ 00772 if(solver_) 00773 delete solver_; 00774 if(coarseSmoother_) 00775 ConstructionTraits<Smoother>::deconstruct(coarseSmoother_); 00776 } 00777 00778 // Postprocess all smoothers 00779 typedef typename Hierarchy<Smoother,A>::Iterator Iterator; 00780 typedef typename Hierarchy<Range,A>::Iterator RIterator; 00781 typedef typename Hierarchy<Domain,A>::Iterator DIterator; 00782 Iterator coarsest = smoothers_.coarsest(); 00783 Iterator smoother = smoothers_.finest(); 00784 DIterator lhs = lhs_->finest(); 00785 if(smoothers_.levels()>0){ 00786 if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels()) 00787 smoother->post(*lhs); 00788 if(smoother!=coarsest) 00789 for(++smoother, ++lhs; smoother != coarsest; ++smoother, ++lhs) 00790 smoother->post(*lhs); 00791 smoother->post(*lhs); 00792 } 00793 00794 delete &(*lhs_->finest()); 00795 delete lhs_; 00796 delete &(*update_->finest()); 00797 delete update_; 00798 delete &(*rhs_->finest()); 00799 delete rhs_; 00800 } 00801 00802 template<class M, class X, class S, class PI, class A> 00803 template<class A1> 00804 void AMG<M,X,S,PI,A>::getCoarsestAggregateNumbers(std::vector<std::size_t,A1>& cont) 00805 { 00806 matrices_->getCoarsestAggregatesOnFinest(cont); 00807 } 00808 00809 } // end namespace Amg 00810 }// end namespace Dune 00811 00812 #endif