diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh index 688eca756368fd2407f39606c513189402ba670d..c307c3b550abe840c5dab277a0e98395dc8d7c07 100644 --- a/dune/istl/paamg/amg.hh +++ b/dune/istl/paamg/amg.hh @@ -175,10 +175,10 @@ namespace Dune /** * @brief Move the iterators to the finer level * @*/ - bool moveToFineLevel(); + void moveToFineLevel(bool processedFineLevel); /** @brief Move the iterators to the coarser level */ - void moveToCoarseLevel(bool processedFineLevel); + bool moveToCoarseLevel(); /** @brief Initialize iterators over levels with fine level */ void initIteratorsWithFineLevel(); @@ -429,7 +429,7 @@ namespace Dune template<class M, class X, class S, class P, class A> bool AMG<M,X,S,P,A> - ::moveToFineLevel() + ::moveToCoarseLevel() { bool processNextLevel=true; @@ -473,7 +473,7 @@ namespace Dune template<class M, class X, class S, class P, class A> void AMG<M,X,S,P,A> - ::moveToCoarseLevel(bool processNextLevel) + ::moveToFineLevel(bool processNextLevel) { if(processNextLevel) { if(matrix != matrices_->matrices().coarsest() || matrices_->levels()<matrices_->maxlevels()) { @@ -580,7 +580,7 @@ namespace Dune presmooth(); #ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION - bool processNextLevel = moveToFineLevel(); + bool processNextLevel = moveToCoarseLevel(); if(processNextLevel) { // next level @@ -588,7 +588,7 @@ namespace Dune mgc(); } - moveToCoarseLevel(processNextLevel); + moveToFineLevel(processNextLevel); #else *lhs=0; #endif diff --git a/dune/istl/paamg/kamg.hh b/dune/istl/paamg/kamg.hh index 5b61987cce049a7a5298670f0c0ac5a55861f92a..9ce4b9cf127b3044b70c52a147cf4377a9f946c6 100644 --- a/dune/istl/paamg/kamg.hh +++ b/dune/istl/paamg/kamg.hh @@ -64,22 +64,22 @@ namespace Dune { *amg_.rhs = d; *amg_.lhs = v; - *amg_.update=0; - amg_.level=0; // Copy data + *amg_.update=0; + amg_.presmooth(); bool processFineLevel = - amg_.moveToFineLevel(); + amg_.moveToCoarseLevel(); if(processFineLevel) { - typename AMG::Domain x=*amg_.update; typename AMG::Range b=*amg_.rhs; + typename AMG::Domain x=*amg_.update; InverseOperatorResult res; coarseSolver_->apply(x, b, res); *amg_.update=x; } - amg_.moveToCoarseLevel(processFineLevel); + amg_.moveToFineLevel(processFineLevel); amg_.postsmooth(); v=*amg_.update; @@ -162,8 +162,8 @@ namespace Dune */ KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, const SmootherArgs& smootherArgs, std::size_t gamma, - std::size_t preSmoothingSteps =1, - std::size_t postSmoothingSteps = 1); + std::size_t preSmoothingSteps =1, std::size_t postSmoothingSteps = 1, + std::size_t maxLevelKrylovSteps = 3 , double minDefectReduction =1e-1); /** * @brief Construct an AMG with an inexact coarse solver based on the smoother. @@ -177,12 +177,15 @@ namespace Dune * @param gamma 1 for V-cycle, 2 for W-cycle * @param preSmoothingSteps The number of smoothing steps for premoothing. * @param postSmoothingSteps The number of smoothing steps for postmoothing. + * @param maxLevelKrylovSteps The maximum number of Krylov steps allowed at each level. + * @param minDefectReduction The defect reduction to achieve on each krylov level. * @param pinfo The information about the parallel distribution of the data. */ template<class C> KAMG(const Operator& fineOperator, const C& criterion, const SmootherArgs& smootherArgs, std::size_t gamma=1, std::size_t preSmoothingSteps=1, std::size_t postSmoothingSteps=1, + std::size_t maxLevelKrylovSteps=3, double minDefectReduction=1e-1, const ParallelInformation& pinfo=ParallelInformation()); /** \copydoc Solver::pre(X&,Y&) */ @@ -191,10 +194,16 @@ namespace Dune void post(Domain& x); /** \copydoc Solver::apply(X&,Y&) */ void apply(Domain& x, const Range& b); - private: /** @brief The underlying amg. */ Amg amg; + + /** \brief The maximum number of Krylov steps allowed at each level. */ + std::size_t maxLevelKrylovSteps; + + /** \brief The defect reduction to achieve on each krylov level. */ + double levelDefectReduction; + /** @brief pointers to the allocated scalar products. */ std::vector<typename Amg::ScalarProduct*> scalarproducts; @@ -206,9 +215,10 @@ namespace Dune KAMG<M,X,S,P,K,A>::KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, const SmootherArgs& smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, - std::size_t postSmoothingSteps) + std::size_t postSmoothingSteps, + std::size_t ksteps, double reduction) : amg(matrices, coarseSolver, smootherArgs, gamma, preSmoothingSteps, - postSmoothingSteps) + postSmoothingSteps), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) {} template<class M, class X, class S, class P, class K, class A> @@ -216,22 +226,24 @@ namespace Dune KAMG<M,X,S,P,K,A>::KAMG(const Operator& fineOperator, const C& criterion, const SmootherArgs& smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, + std::size_t ksteps, double reduction, const ParallelInformation& pinfo) : amg(fineOperator, criterion, smootherArgs, gamma, preSmoothingSteps, - postSmoothingSteps, false, pinfo) + postSmoothingSteps, false, pinfo), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) {} + template<class M, class X, class S, class P, class K, class A> void KAMG<M,X,S,P,K,A>::pre(Domain& x, Range& b) { amg.pre(x,b); scalarproducts.reserve(amg.levels()); ksolvers.reserve(amg.levels()); + typename OperatorHierarchy::ParallelMatrixHierarchy::Iterator matrix = amg.matrices_->matrices().coarsest(); typename ParallelInformationHierarchy::Iterator pinfo = amg.matrices_->parallelInformation().coarsest(); - bool hasCoarsest=(amg.levels()==amg.maxlevels()); KrylovSolver* ks=0; @@ -244,13 +256,14 @@ namespace Dune }else ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); - double reduction = 1e-2; - int maxit = 3; + std::ostringstream s; + if(matrix!=amg.matrices_->matrices().finest()) while(true) { scalarproducts.push_back(Amg::ScalarProductChooser::construct(*pinfo)); ks = new KrylovSolver(*matrix, *(scalarproducts.back()), - *(ksolvers.back()), reduction, maxit, 0); + *(ksolvers.back()), levelDefectReduction, + maxLevelKrylovSteps, 0); ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); --matrix; --pinfo;