From b4a6dbfa65c5ce85ce358fff3a085030b7e8f992 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 20 Jun 2015 22:22:56 +0100 Subject: [PATCH] Patches for beginnings of an overlap multigrid --- lib/Algorithms.h | 4 +- lib/algorithms/CoarsenedMatrix.h | 126 ++++++++++++++++--- lib/algorithms/iterative/ConjugateGradient.h | 25 ++-- lib/algorithms/iterative/ConjugateResidual.h | 15 +-- lib/lattice/Lattice_ET.h | 2 +- lib/lattice/Lattice_transfer.h | 2 +- lib/qcd/utils/SUn.h | 54 +++++++- tests/Test_cayley_coarsen_support.cc | 19 ++- tests/Test_cayley_ldop_cg.cc | 49 +++++--- tests/Test_cayley_ldop_cr.cc | 63 +++++++--- tests/Test_cayley_ldop_cr_chebyshev.cc | 124 ------------------ 11 files changed, 285 insertions(+), 198 deletions(-) delete mode 100644 tests/Test_cayley_ldop_cr_chebyshev.cc diff --git a/lib/Algorithms.h b/lib/Algorithms.h index 2007cf60..751cf6b2 100644 --- a/lib/Algorithms.h +++ b/lib/Algorithms.h @@ -3,7 +3,6 @@ #include #include -#include #include #include @@ -17,6 +16,9 @@ #include + +#include + // Eigen/lanczos // EigCg // MCR diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h index fe801106..2576a696 100644 --- a/lib/algorithms/CoarsenedMatrix.h +++ b/lib/algorithms/CoarsenedMatrix.h @@ -12,9 +12,6 @@ namespace Grid { std::vector directions ; std::vector displacements; - // FIXME -- don't like xposing the operator directions - // as different to the geometrical dirs - // Also don't like special casing five dim.. should pass an object in template Geometry(int _d) { int base = (_d==5) ? 1:0; @@ -64,6 +61,99 @@ namespace Grid { }; + template + class Aggregation { + public: + typedef iVector siteVector; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; + + typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field + typedef Lattice FineField; + + GridBase *CoarseGrid; + GridBase *FineGrid; + std::vector > subspace; + + Aggregation(GridBase *_CoarseGrid,GridBase *_FineGrid) : + CoarseGrid(_CoarseGrid), + FineGrid(_FineGrid), + subspace(nbasis,_FineGrid) + { + }; + + void Orthogonalise(void){ + CoarseScalar InnerProd(CoarseGrid); + blockOrthogonalise(InnerProd,subspace); +#if 1 + // CheckOrthogonal(); +#endif + + } + void CheckOrthogonal(void){ + CoarseVector iProj(CoarseGrid); + CoarseVector eProj(CoarseGrid); + Lattice pokey(CoarseGrid); + + + for(int i=0;ioSites();ss++){ + eProj._odata[ss](i)=CComplex(1.0); + } + eProj=eProj - iProj; + std::cout<<"Orthog check error "< &hermop) { + + RealD scale; + + ConjugateGradient CG(1.0e-4,10000); + FineField noise(FineGrid); + FineField Mn(FineGrid); + + for(int b=0;b "< @@ -145,7 +235,8 @@ namespace Grid { comm_buf.resize(Stencil._unified_buffer_size); }; - void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase > &linop,std::vector > & subspace){ + void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase > &linop, + Aggregation & Subspace){ FineField iblock(FineGrid); // contributions from within this block FineField oblock(FineGrid); // contributions from outwith this block @@ -162,8 +253,11 @@ namespace Grid { CoarseScalar InnerProd(Grid()); // Orthogonalise the subblocks over the basis - blockOrthogonalise(InnerProd,subspace); - blockProject(iProj,subspace[0],subspace); + blockOrthogonalise(InnerProd,Subspace.subspace); + //Subspace.Orthogonalise(); + // Subspace.CheckOrthogonal(); + //Subspace.Orthogonalise(); + // Subspace.CheckOrthogonal(); // Compute the matrix elements of linop between this orthonormal // set of vectors. @@ -177,7 +271,7 @@ namespace Grid { assert(self_stencil!=-1); for(int i=0;ioSites();ss++){ for(int j=0;j bc(FineGrid->_ndimension,0); blockPick(Grid(),phi,tmp,bc); // Pick out a block linop.Op(tmp,Mphi); // Apply big dop - blockProject(iProj,Mphi,subspace); // project it and print it + blockProject(iProj,Mphi,Subspace.subspace); // project it and print it std::cout<< " Computed matrix elements from block zero only "< ident; ident=one; + iMatrix ident; ident=one; val = val*adj(val); val = val + 1.0; @@ -279,7 +375,7 @@ namespace Grid { int dd=d+1; A[2*d] = adj(Cshift(A[2*d+1],dd,1)); } - A[8] = 0.5*(A[8] + adj(A[8])); + // A[8] = 0.5*(A[8] + adj(A[8])); } void AssertHermitian(void) { CoarseMatrix AA (Grid()); diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h index e6915d83..0f366ebf 100644 --- a/lib/algorithms/iterative/ConjugateGradient.h +++ b/lib/algorithms/iterative/ConjugateGradient.h @@ -15,7 +15,7 @@ public: Integer MaxIterations; int verbose; ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { - verbose=1; + verbose=0; }; @@ -58,7 +58,7 @@ public: return; } - std::cout << std::setprecision(4)<< "ConjugateGradient: k=0 residual "< &Linop,const Field &src, Field &psi){ @@ -37,8 +37,8 @@ namespace Grid { Linop.HermOpAndNorm(p,Ap,pAp,pAAp); Linop.HermOpAndNorm(r,Ar,rAr,rAAr); - std::cout << "pAp, pAAp"<< pAp<<" "< &ip,Lattice &fineX) GridBase *coarse = ip._grid; Lattice zz(fineX._grid); zz=zero; blockInnerProduct(ip,fineX,fineX); - ip = rsqrt(ip); + ip = pow(ip,-0.5); blockZAXPY(fineX,ip,fineX,zz); } // useful in multigrid project; diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index 0ede7965..647699c8 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -522,22 +522,66 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g } // reunitarise?? + static void LieRandomize(GridParallelRNG &pRNG,LatticeMatrix &out,double scale=1.0){ + GridBase *grid = out._grid; + LatticeComplex ca (grid); + LatticeMatrix lie(grid); + LatticeMatrix la (grid); + Complex ci(0.0,scale); + Matrix ta; + + lie=zero; + for(int a=0;a subspace(nbasis,FGrid); LatticeFermion prom(FGrid); - for(int b=0;b subspace(nbasis,FGrid); + + std::cout<<"Calling Aggregation class" < HermDefOp(Ddwf); + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FGrid); + Aggregates.CreateSubspaceRandom(RNG5); + + subspace=Aggregates.subspace; + + std::cout << "Called aggregation class"<< std::endl; typedef CoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; LittleDiracOperator LittleDiracOp(*Coarse5d); - LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace); + LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); CoarseVector c_src (Coarse5d); CoarseVector c_res (Coarse5d); diff --git a/tests/Test_cayley_ldop_cg.cc b/tests/Test_cayley_ldop_cg.cc index b146002d..f960ab7b 100644 --- a/tests/Test_cayley_ldop_cg.cc +++ b/tests/Test_cayley_ldop_cg.cc @@ -1,4 +1,6 @@ #include +#include +#include using namespace std; using namespace Grid; @@ -59,67 +61,84 @@ int main (int argc, char ** argv) LatticeFermion ref(FGrid); ref=zero; LatticeFermion tmp(FGrid); LatticeFermion err(FGrid); - LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu); + LatticeGaugeField Umu(UGrid); + + //gaussian(RNG4,Umu); //random(RNG4,Umu); + NerscField header; std::string file("./ckpoint_lat.400"); readNerscConfiguration(Umu,header,file); + // SU3::ColdConfiguration(RNG4,Umu); + // SU3::TepidConfiguration(RNG4,Umu); + // SU3::HotConfiguration(RNG4,Umu); + // Umu=zero; #if 0 LatticeColourMatrix U(UGrid); - U=zero; for(int nn=0;nn2) { - pokeIndex(Umu,U,nn); - } + U=peekIndex(Umu,nn); + U=U*adj(U)-1.0; + std::cout<<"SU3 test "< HermIndefOp(Ddwf); const int nbasis = 8; + +#if 0 std::vector subspace(nbasis,FGrid); LatticeFermion noise(FGrid); LatticeFermion ms(FGrid); for(int b=0;b HermDefOp(Ddwf); - ConjugateGradient CG(1.0e-5,10000); + ConjugateGradient CG(1.0e-4,10000); - for(int i=0;i<4;i++){ + for(int i=0;i<1;i++){ CG(HermDefOp,noise,subspace[b]); noise = subspace[b]; scale = pow(norm2(noise),-0.5); noise=noise*scale; - HermDefOp.HermOp(noise,ms); std::cout << "filt "< "< "< HermDefOp(Ddwf); + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FGrid); + Aggregates.CreateSubspace(RNG5,HermDefOp); + std::cout << "Called aggregation class"<< std::endl; +#endif + typedef CoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; LittleDiracOperator LittleDiracOp(*Coarse5d); - LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace); + LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); CoarseVector c_src (Coarse5d); CoarseVector c_res (Coarse5d); @@ -128,7 +147,7 @@ int main (int argc, char ** argv) std::cout << "Solving CG on coarse space "<< std::endl; MdagMLinearOperator PosdefLdop(LittleDiracOp); - ConjugateGradient CG(1.0e-8,10000); + ConjugateGradient CG(1.0e-6,10000); CG(PosdefLdop,c_src,c_res); std::cout << "Done "<< std::endl; diff --git a/tests/Test_cayley_ldop_cr.cc b/tests/Test_cayley_ldop_cr.cc index 698bbb28..c4001b55 100644 --- a/tests/Test_cayley_ldop_cr.cc +++ b/tests/Test_cayley_ldop_cr.cc @@ -1,4 +1,6 @@ #include +#include +#include using namespace std; using namespace Grid; @@ -59,67 +61,100 @@ int main (int argc, char ** argv) LatticeFermion ref(FGrid); ref=zero; LatticeFermion tmp(FGrid); LatticeFermion err(FGrid); - LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu); + LatticeGaugeField Umu(UGrid); + + //gaussian(RNG4,Umu); //random(RNG4,Umu); - NerscField header; - std::string file("./ckpoint_lat.4000"); - readNerscConfiguration(Umu,header,file); - RealD mass=0.01; - RealD M5=1.8; + NerscField header; + std::string file("./ckpoint_lat.400"); + readNerscConfiguration(Umu,header,file); + // SU3::ColdConfiguration(RNG4,Umu); + // SU3::TepidConfiguration(RNG4,Umu); + // SU3::HotConfiguration(RNG4,Umu); + // Umu=zero; + +#if 0 + LatticeColourMatrix U(UGrid); + for(int nn=0;nn(Umu,nn); + U=U*adj(U)-1.0; + std::cout<<"SU3 test "< HermIndefOp(Ddwf); const int nbasis = 8; + +#if 0 std::vector subspace(nbasis,FGrid); LatticeFermion noise(FGrid); LatticeFermion ms(FGrid); for(int b=0;b HermDefOp(Ddwf); ConjugateGradient CG(1.0e-4,10000); - for(int i=0;i<4;i++){ + for(int i=0;i<1;i++){ CG(HermDefOp,noise,subspace[b]); noise = subspace[b]; scale = pow(norm2(noise),-0.5); noise=noise*scale; - HermDefOp.HermOp(noise,ms); std::cout << "filt "< "< "< HermDefOp(Ddwf); + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FGrid); + Aggregates.CreateSubspace(RNG5,HermDefOp); + std::cout << "Called aggregation class"<< std::endl; +#endif + typedef CoarsenedMatrix LittleDiracOperator; typedef LittleDiracOperator::CoarseVector CoarseVector; LittleDiracOperator LittleDiracOp(*Coarse5d); - LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace); + LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); CoarseVector c_src (Coarse5d); CoarseVector c_res (Coarse5d); gaussian(CRNG,c_src); c_res=zero; + + std::cout << "Solving CG on coarse space "<< std::endl; + MdagMLinearOperator PosdefLdop(LittleDiracOp); + ConjugateGradient CG(1.0e-6,10000); + CG(PosdefLdop,c_src,c_res); + + std::cout << "Solving MCR on coarse space "<< std::endl; HermitianLinearOperator HermIndefLdop(LittleDiracOp); - ConjugateResidual MCR(1.0e-8,10000); + ConjugateResidual MCR(1.0e-6,10000); MCR(HermIndefLdop,c_src,c_res); std::cout << "Done "<< std::endl; diff --git a/tests/Test_cayley_ldop_cr_chebyshev.cc b/tests/Test_cayley_ldop_cr_chebyshev.cc deleted file mode 100644 index 091b8eff..00000000 --- a/tests/Test_cayley_ldop_cr_chebyshev.cc +++ /dev/null @@ -1,124 +0,0 @@ -#include - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - - -template -struct scal { - d internal; -}; - -Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT -}; - -double lowpass(double x) -{ - return pow(x*x,-2); -} - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - Chebyshev filter(-120.0, 120.0,256, lowpass); - ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc); - filter.csv(csv); - csv.close(); - - const int Ls=8; - - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi()); - GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - - // Construct a coarsened grid - std::vector clatt = GridDefaultLatt(); - for(int d=0;d seeds4({1,2,3,4}); - std::vector seeds5({5,6,7,8}); - std::vector cseeds({5,6,7,8}); - GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); - GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); - - LatticeFermion src(FGrid); gaussian(RNG5,src); - LatticeFermion result(FGrid); result=zero; - LatticeFermion ref(FGrid); ref=zero; - LatticeFermion tmp(FGrid); - LatticeFermion err(FGrid); - LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu); - - - //random(RNG4,Umu); - NerscField header; - std::string file("./ckpoint_lat.400"); - readNerscConfiguration(Umu,header,file); - -#if 0 - LatticeColourMatrix U(UGrid); - Complex cone(1.0,0.0); - for(int nn=0;nn(Umu,U,nn); - } - } -#endif - RealD mass=0.5; - RealD M5=1.8; - - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); - - const int nbasis = 8; - std::vector subspace(nbasis,FGrid); - LatticeFermion noise(FGrid); - LatticeFermion ms(FGrid); - for(int b=0;b LittleDiracOperator; - typedef LittleDiracOperator::CoarseVector CoarseVector; - - LittleDiracOperator LittleDiracOp(*Coarse5d); - LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace); - - CoarseVector c_src (Coarse5d); - CoarseVector c_res (Coarse5d); - gaussian(CRNG,c_src); - c_res=zero; - - std::cout << "Solving MCR on coarse space "<< std::endl; - HermitianLinearOperator HermIndefLdop(LittleDiracOp); - ConjugateResidual MCR(1.0e-8,10000); - MCR(HermIndefLdop,c_src,c_res); - - std::cout << "Done "<< std::endl; - Grid_finalize(); -}