diff --git a/dune/istl/gsetc.hh b/dune/istl/gsetc.hh index 2770cf4ca88cff1d9336ebde635e23541b9e2892..821426aa1bb89e1886d292247c31a94e6bfe5349 100644 --- a/dune/istl/gsetc.hh +++ b/dune/istl/gsetc.hh @@ -373,7 +373,7 @@ namespace Dune { template<int I> struct algmeta_itsteps { -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION template<typename T11, typename T12, typename T13, typename T14, @@ -421,7 +421,7 @@ namespace Dune { x.axpy(K(1)-w,xold); } -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION template<typename T11, typename T12, typename T13, typename T14, @@ -471,7 +471,7 @@ namespace Dune { } } -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION template<typename T11, typename T12, typename T13, typename T14, @@ -521,7 +521,7 @@ namespace Dune { } } -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION template<typename T11, typename T12, typename T13, typename T14, diff --git a/dune/istl/multitypeblockmatrix.hh b/dune/istl/multitypeblockmatrix.hh index da0012ed9e7ea04fe622a2ca893278d4a843dfc2..eb7ab3e2d32ad04caedc588c3f7328958b8dd405 100644 --- a/dune/istl/multitypeblockmatrix.hh +++ b/dune/istl/multitypeblockmatrix.hh @@ -8,7 +8,7 @@ #include "istlexception.hh" -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION #include <boost/fusion/sequence.hpp> @@ -463,5 +463,5 @@ namespace Dune { } // end namespace #endif // HAVE_BOOST_FUSION -#endif // HAVE_DUNE_BOOST +#endif // HAVE_BOOST #endif diff --git a/dune/istl/multitypeblockvector.hh b/dune/istl/multitypeblockvector.hh index 39495b9dee762e6988aade0c730dd1b7e8fd32c6..63da325dc25c821f2cc12c0f7d3bc07102ec60f0 100644 --- a/dune/istl/multitypeblockvector.hh +++ b/dune/istl/multitypeblockvector.hh @@ -3,7 +3,7 @@ #ifndef DUNE_MULTITYPEVECTOR_HH #define DUNE_MULTITYPEVECTOR_HH -#if HAVE_DUNE_BOOST +#if HAVE_BOOST #ifdef HAVE_BOOST_FUSION #include <cmath> @@ -317,6 +317,6 @@ namespace Dune { } // end namespace #endif // end HAVE_BOOST_FUSION -#endif // end HAVE_DUNE_BOOST +#endif // end HAVE_BOOST #endif diff --git a/dune/istl/paamg/kamg.hh b/dune/istl/paamg/kamg.hh index 7ea418399a36e598936b2f530e13e24da9d1f1fc..eacf0c5a024d58269b6046cab1199e3d611c60fe 100644 --- a/dune/istl/paamg/kamg.hh +++ b/dune/istl/paamg/kamg.hh @@ -47,9 +47,8 @@ namespace Dune * @param coarseSolver The solver used for the coarse grid correction. */ - KAmgTwoGrid(AMG& amg, InverseOperator<Domain,Range>* coarseSolver, bool - isCoarsest) - : amg_(amg), coarseSolver_(coarseSolver), isCoarsest_(isCoarsest) + KAmgTwoGrid(AMG& amg, InverseOperator<Domain,Range>* coarseSolver) + : amg_(amg), coarseSolver_(coarseSolver) {} /** \copydoc Preconditioner::pre(X&,Y&) */ @@ -63,29 +62,27 @@ namespace Dune /** \copydoc Preconditioner::apply(X&,const Y&) */ void apply(typename AMG::Domain& v, const typename AMG::Range& d) { - Range rhs(d); + *amg_.rhs = d; *amg_.lhs = v; // Copy data + *amg_.update=0; + amg_.presmooth(); bool processFineLevel = - amg_.moveToCoarseLevel(rhs); + amg_.moveToCoarseLevel(); if(processFineLevel) { - typename AMG::Domain x(*amg_.update=0); // update is zero - if(!isCoarsest_) - amg_.presmooth(); - typename AMG::Range b(*amg_.rhs); // b might get overridden + typename AMG::Range b=*amg_.rhs; + typename AMG::Domain x=*amg_.update; InverseOperatorResult res; coarseSolver_->apply(x, b, res); - *amg_.update+=x; // coarse level correction has to be in lhs - - if(!isCoarsest_) { - *amg_.lhs=x; - amg_.postsmooth(); - } + *amg_.update=x; } - amg_.moveToFineLevel(v, processFineLevel); + amg_.moveToFineLevel(processFineLevel); + + amg_.postsmooth(); + v=*amg_.update; } /** @@ -106,7 +103,6 @@ namespace Dune AMG& amg_; /** @brief The coarse grid solver.*/ InverseOperator<Domain,Range>* coarseSolver_; - bool isCoarsest_; }; @@ -167,9 +163,9 @@ namespace Dune * @param minDefectReduction The minimal defect reduction to achieve on each Krylov level. */ KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, - const SmootherArgs& smootherArgs, std::size_t maxLevelKrylovSteps = 3 , + const SmootherArgs& smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps =1, std::size_t postSmoothingSteps = 1, - double minDefectReduction =1e-1); + std::size_t maxLevelKrylovSteps = 3 , double minDefectReduction =1e-1); /** * @brief Construct an AMG with an inexact coarse solver based on the smoother. @@ -189,9 +185,9 @@ namespace Dune */ template<class C> KAMG(const Operator& fineOperator, const C& criterion, - const SmootherArgs& smootherArgs, std::size_t maxLevelKrylovSteps=3, + const SmootherArgs& smootherArgs, std::size_t gamma=1, std::size_t preSmoothingSteps=1, std::size_t postSmoothingSteps=1, - double minDefectReduction=1e-1, + std::size_t maxLevelKrylovSteps=3, double minDefectReduction=1e-1, const ParallelInformation& pinfo=ParallelInformation()); /** \copydoc Preconditioner::pre(X&,Y&) */ @@ -217,30 +213,27 @@ namespace Dune std::vector<typename Amg::ScalarProduct*> scalarproducts; /** @brief pointers to the allocated krylov solvers. */ - std::vector<KrylovSolver*> ksolvers; - - /** @brief pointers to the allocated krylov solvers. */ - std::vector<KAmgTwoGrid<Amg>*> twogridMethods; + std::vector<KAmgTwoGrid<Amg>*> ksolvers; }; template<class M, class X, class S, class P, class K, class A> KAMG<M,X,S,P,K,A>::KAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver, const SmootherArgs& smootherArgs, - std::size_t ksteps, std::size_t preSmoothingSteps, + std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, - double reduction) - : amg(matrices, coarseSolver, smootherArgs, 1, preSmoothingSteps, + std::size_t ksteps, double reduction) + : amg(matrices, coarseSolver, smootherArgs, gamma, preSmoothingSteps, postSmoothingSteps), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) {} template<class M, class X, class S, class P, class K, class A> template<class C> KAMG<M,X,S,P,K,A>::KAMG(const Operator& fineOperator, const C& criterion, - const SmootherArgs& smootherArgs, std::size_t ksteps, + const SmootherArgs& smootherArgs, std::size_t gamma, std::size_t preSmoothingSteps, std::size_t postSmoothingSteps, - double reduction, + std::size_t ksteps, double reduction, const ParallelInformation& pinfo) - : amg(fineOperator, criterion, smootherArgs, 1, preSmoothingSteps, + : amg(fineOperator, criterion, smootherArgs, gamma, preSmoothingSteps, postSmoothingSteps, false, pinfo), maxLevelKrylovSteps(ksteps), levelDefectReduction(reduction) {} @@ -257,43 +250,31 @@ namespace Dune typename ParallelInformationHierarchy::Iterator pinfo = amg.matrices_->parallelInformation().coarsest(); bool hasCoarsest=(amg.levels()==amg.maxlevels()); + KrylovSolver* ks=0; - if(hasCoarsest && matrix==amg.matrices_->matrices().finest()) - return; - - --matrix; - --pinfo; - twogridMethods.push_back(new KAmgTwoGrid<Amg>(amg, amg.solver_,true)); - ksolvers.push_back(new KrylovSolver(*matrix, *twogridMethods.back(), - levelDefectReduction,maxLevelKrylovSteps,1)); - std::cout<<"push back true"<<std::endl; - if(matrix!=amg.matrices_->matrices().finest()) - for(--matrix, --pinfo; true; --matrix, --pinfo) { - twogridMethods.push_back(new KAmgTwoGrid<Amg>(amg, ksolvers.back(),false)); - ksolvers.push_back(new KrylovSolver(*matrix, *twogridMethods.back(), - levelDefectReduction,maxLevelKrylovSteps,1)); - std::cout<<"push back false"<<std::endl; - if(matrix==amg.matrices_->matrices().finest()) - break; - } + if(hasCoarsest) { + if(matrix==amg.matrices_->matrices().finest()) + return; + --matrix; + --pinfo; + ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, amg.solver_)); + }else + ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); - /* - std::ostringstream s; + std::ostringstream s; - while(true){ + if(matrix!=amg.matrices_->matrices().finest()) + while(true) { scalarproducts.push_back(Amg::ScalarProductChooser::construct(*pinfo)); ks = new KrylovSolver(*matrix, *(scalarproducts.back()), - *(ksolvers.back()), levelDefectReduction, + *(ksolvers.back()), levelDefectReduction, maxLevelKrylovSteps, 0); + ksolvers.push_back(new KAmgTwoGrid<Amg>(amg, ks)); + --matrix; + --pinfo; if(matrix==amg.matrices_->matrices().finest()) break; - else{ - --matrix; - --pinfo; - ksolvers.push_front(new KAmgTwoGrid<Amg>(amg, ks)); - std::cout<<"push back false 2"<<std::endl; - } - }*/ + } } @@ -302,7 +283,7 @@ namespace Dune { typedef typename std::vector<KAmgTwoGrid<Amg>*>::reverse_iterator KIterator; typedef typename std::vector<ScalarProduct*>::iterator SIterator; - for(KIterator kiter = twogridMethods.rbegin(); kiter != twogridMethods.rend(); + for(KIterator kiter = ksolvers.rbegin(); kiter != ksolvers.rend(); ++kiter) { if((*kiter)->coarseSolver()!=amg.solver_) @@ -320,26 +301,13 @@ namespace Dune void KAMG<M,X,S,P,K,A>::apply(Domain& v, const Range& d) { amg.initIteratorsWithFineLevel(); - *amg.lhs =v; - *amg.rhs = d; - *amg.update = 0; if(ksolvers.size()==0) { Range td=d; InverseOperatorResult res; amg.solver_->apply(v,td,res); - }else{ - amg.presmooth(); - *amg.lhs=0; -#ifndef DUNE_AMG_NO_COARSEGRIDCORRECTION - Range td=*amg.rhs; - InverseOperatorResult res; - ksolvers.back()->apply(*amg.lhs,td, res); - *amg.update+=*amg.lhs; -#endif - amg.postsmooth(); - v=*amg.update; - } + }else + ksolvers.back()->apply(v,d); } template<class M, class X, class S, class P, class K, class A> diff --git a/dune/istl/paamg/parameters.hh b/dune/istl/paamg/parameters.hh index f6a1ba3e6e86cd307bcddd766a6c9abc8d7d94e3..edff1ef89b69bdf37eef814f08134dd449a3bce4 100644 --- a/dune/istl/paamg/parameters.hh +++ b/dune/istl/paamg/parameters.hh @@ -106,7 +106,7 @@ namespace Dune */ void setDefaultValuesIsotropic(std::size_t dim, std::size_t diameter=2) { - maxDistance_=0; + maxDistance_=diameter-1; std::size_t csize=1; for(; dim>0; dim--) { diff --git a/dune/istl/paamg/test/amgtest.cc b/dune/istl/paamg/test/amgtest.cc index b8c5f869b43d5c2f5d7c7acc6fe54369b46ff6f8..e69405cbddd26353988dfc577fab19afa14999fa 100644 --- a/dune/istl/paamg/test/amgtest.cc +++ b/dune/istl/paamg/test/amgtest.cc @@ -77,10 +77,15 @@ void testAMG(int N, int coarsenTarget, int ml) int n; Comm c; - BCRSMat mat = setupAnisotropic2d<BS,XREAL>(N, indices, c, &n, 1000); + BCRSMat mat = setupAnisotropic2d<BS,XREAL>(N, indices, c, &n, 1); Vector b(mat.N()), x(mat.M()); + b=0; + x=100; + + setBoundary(x, b, N); + x=0; randomize(mat, b); @@ -94,12 +99,10 @@ void testAMG(int N, int coarsenTarget, int ml) watch.reset(); Operator fop(mat); - //typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> > - // Criterion; - typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<BCRSMat,Dune::Amg::FirstDiagonal> > CriterionBase; - typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion; - //typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother; - typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother; + typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> > + Criterion; + typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother; + //typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother; //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother; //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother; //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother; @@ -116,12 +119,11 @@ void testAMG(int N, int coarsenTarget, int ml) smootherArgs.relaxationFactor = 1; - Criterion criterion(Dune::Amg::Parameters(15,coarsenTarget)); - criterion.setDefaultValuesIsotropic(2, 2); + Criterion criterion(15,coarsenTarget); + criterion.setDefaultValuesIsotropic(2); criterion.setAlpha(.67); criterion.setBeta(1.0e-4); criterion.setMaxLevel(ml); - criterion.setProlongationDampingFactor(1.6); criterion.setSkipIsolated(false); Dune::SeqScalarProduct<Vector> sp; diff --git a/dune/istl/paamg/test/anisotropic.hh b/dune/istl/paamg/test/anisotropic.hh index 065f7bb40290710d0f3a842d9d546c5fac404743..678b97febb44baedd5b86fc657c0722f8f7acc87 100644 --- a/dune/istl/paamg/test/anisotropic.hh +++ b/dune/istl/paamg/test/anisotropic.hh @@ -105,7 +105,7 @@ void fillValues(int N, M& mat, int overlapStart, int overlapEnd, int start, int ColIterator endi = (*i).end(); - if(x<start || x >= end || x==0 || x==N-1 || y==0 || y==N-1) { + if(x<start || x >= end) { // || x==0 || x==N-1 || y==0 || y==N-1){ // overlap node is dirichlet border ColIterator j = (*i).begin(); @@ -140,7 +140,7 @@ void setBoundary(Dune::BlockVector<Dune::FieldVector<T,BS> >& lhs, if(x==0 || y ==0 || x==n-1 || y==n-1) { //double h = 1.0 / ((double) (n-1)); - lhs[i->local()]=rhs[i->local()]; //((double)x)*((double)y)*h*h; + lhs[i->local()]=rhs[i->local()]=0; //((double)x)*((double)y)*h*h; } } } @@ -164,7 +164,7 @@ void setBoundary(Dune::BlockVector<Dune::FieldVector<T,BS> >& lhs, else y = ((T) j)*h; - lhs[j*N+i]=rhs[j*N+i]; //=0;//x*y; + lhs[j*N+i]=rhs[j*N+i]=0; //x*y; } } diff --git a/dune/istl/paamg/test/kamgtest.cc b/dune/istl/paamg/test/kamgtest.cc index 9fd12f4dc55726a025cc234dd41bec729d5a2d4d..c50ccc3cdc5aa847b33b39f8a7a14a7a59f0aaa5 100644 --- a/dune/istl/paamg/test/kamgtest.cc +++ b/dune/istl/paamg/test/kamgtest.cc @@ -7,7 +7,6 @@ #include <dune/common/parallel/collectivecommunication.hh> #include <dune/istl/paamg/kamg.hh> #include <dune/istl/paamg/pinfo.hh> -#include <dune/istl/matrixmarket.hh> #include <cstdlib> #include <ctime> @@ -45,7 +44,7 @@ void testAMG(int N, int coarsenTarget, int ml) int n; Comm c; - BCRSMat mat = setupAnisotropic2d<BS,double>(N, indices, c, &n, 1); + BCRSMat mat = setupAnisotropic2d<BS>(N, indices, c, &n, .001); Vector b(mat.N()), x(mat.M()); @@ -56,16 +55,10 @@ void testAMG(int N, int coarsenTarget, int ml) x=0; randomize(mat, b); - x=1; - mat.mv(x, b); - x=0; if(N<6) { Dune::printmatrix(std::cout, mat, "A", "row"); Dune::printvector(std::cout, x, "x", "row"); - }else{ - Dune::storeMatrixMarket(mat, "Mat.mm"); - Dune::storeMatrixMarket(b, "b.mm"); } Dune::Timer watch; @@ -94,7 +87,7 @@ void testAMG(int N, int coarsenTarget, int ml) smootherArgs.relaxationFactor = 1; Criterion criterion(15,coarsenTarget); - criterion.setDefaultValuesIsotropic(2,2); + criterion.setDefaultValuesIsotropic(2); criterion.setAlpha(.67); criterion.setBeta(1.0e-4); criterion.setMaxLevel(ml); @@ -150,6 +143,6 @@ int main(int argc, char** argv) ml = atoi(argv[3]); testAMG<1>(N, coarsenTarget, ml); - //testAMG<2>(N, coarsenTarget, ml); + testAMG<2>(N, coarsenTarget, ml); } diff --git a/dune/istl/test/Makefile.am b/dune/istl/test/Makefile.am index 2608ced727b0864d2685d99d7d3f24c3ed9f0908..bad932b37ae77056579dbd8f8b47be374f67c1c2 100644 --- a/dune/istl/test/Makefile.am +++ b/dune/istl/test/Makefile.am @@ -30,7 +30,6 @@ NORMALTESTS = basearraytest \ mv \ scaledidmatrixtest \ seqmatrixmarkettest \ - solvertest \ vbvectortest # list of tests to run (indicestest is special case) @@ -106,8 +105,6 @@ iotest_SOURCES = iotest.cc scaledidmatrixtest_SOURCES = scaledidmatrixtest.cc -solvertest_SOURCES = solvertest.cc - if MPI vectorcommtest_SOURCES = vectorcommtest.cc vectorcommtest_CPPFLAGS = $(AM_CPPFLAGS) \