diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh index 1e05f663ef04c369946628b3ba43d003eb48121f..3825ce851a4bcea81f4fe98bf5113aebd7c82bf2 100644 --- a/dune/istl/paamg/amg.hh +++ b/dune/istl/paamg/amg.hh @@ -1057,8 +1057,14 @@ namespace Dune if(processNextLevel) { // next level - for(std::size_t i=0; i<gamma_; i++) + for(std::size_t i=0; i<gamma_; i++){ mgc(levelContext); + if (levelContext.matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) + break; + if(i+1 < gamma_){ + levelContext.matrix->applyscaleadd(-1., *levelContext.lhs, *levelContext.rhs); + } + } } moveToFineLevel(levelContext, processNextLevel); diff --git a/dune/istl/paamg/test/amgtest.cc b/dune/istl/paamg/test/amgtest.cc index b8dbd0e9a52f8ad16b931e04602a6fe2a8d6ac81..c7e1843a0d08c97a10cf75e16d61e7b2fe471372 100644 --- a/dune/istl/paamg/test/amgtest.cc +++ b/dune/istl/paamg/test/amgtest.cc @@ -77,7 +77,7 @@ void randomize(const M& mat, V& b) template <class Matrix, class Vector> -void testAMG(int N, int coarsenTarget, int ml) +Dune::InverseOperatorResult testAMG(int N, int coarsenTarget, int ml, int gamma = 1) { std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl; @@ -139,6 +139,7 @@ void testAMG(int N, int coarsenTarget, int ml) criterion.setDefaultValuesIsotropic(2); criterion.setAlpha(.67); criterion.setBeta(1.0e-4); + criterion.setGamma(gamma); criterion.setMaxLevel(ml); criterion.setSkipIsolated(false); // specify pre/post smoother steps @@ -176,6 +177,7 @@ void testAMG(int N, int coarsenTarget, int ml) std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl; */ + return r; } @@ -195,11 +197,19 @@ try if(argc>3) ml = atoi(argv[3]); - { - using Matrix = Dune::BCRSMatrix<XREAL>; - using Vector = Dune::BlockVector<XREAL>; - - testAMG<Matrix,Vector>(N, coarsenTarget, ml); + Dune::InverseOperatorResult gamma1_res; + for(int gamma = 1; gamma<=2;++gamma){ + { + using Matrix = Dune::BCRSMatrix<XREAL>; + using Vector = Dune::BlockVector<XREAL>; + + Dune::InverseOperatorResult res = testAMG<Matrix,Vector>(N, coarsenTarget, ml, gamma); + if(gamma==1){ + gamma1_res = res; + }else{ + assert(res.conv_rate < gamma1_res.conv_rate); + } + } } {