portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

PVdagM MG

Browse Source
This commit is contained in:
Peter Boyle 2025-02-01 11:04:13 -05:00
parent c4fc972fec
commit c74d11e3d7
3 changed files with 112 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h
View File

@ -74,7 +74,7 @@ public:
void operator() (const Field &src, Field &psi){ void operator() (const Field &src, Field &psi){
psi=Zero(); // psi=Zero();
RealD cp, ssq,rsq; RealD cp, ssq,rsq;
ssq=norm2(src); ssq=norm2(src);
rsq=Tolerance*Tolerance*ssq; rsq=Tolerance*Tolerance*ssq;

74
Grid/algorithms/multigrid/Aggregates.h
View File

@ -30,6 +30,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#pragma once #pragma once
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
inline RealD AggregatePowerLaw(RealD x) inline RealD AggregatePowerLaw(RealD x)
@ -124,6 +126,53 @@ public:
} }
} }
virtual void CreateSubspaceGCR(GridParallelRNG &RNG,LinearOperatorBase<FineField> &DiracOp,int nn=nbasis)
{
RealD scale;
TrivialPrecon<FineField> simple_fine;
PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
FineField noise(FineGrid);
FineField src(FineGrid);
FineField guess(FineGrid);
FineField Mn(FineGrid);
for(int b=0;b<nn;b++){
subspace[b] = Zero();
gaussian(RNG,noise);
scale = std::pow(norm2(noise),-0.5);
noise=noise*scale;
DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|Op|n> "<<innerProduct(noise,Mn)<<std::endl;
for(int i=0;i<3;i++){
// void operator() (const Field &src, Field &psi){
#if 1
std::cout << GridLogMessage << " inverting on noise "<<std::endl;
src = noise;
guess=Zero();
GCR(src,guess);
subspace[b] = guess;
#else
std::cout << GridLogMessage << " inverting on zero "<<std::endl;
src=Zero();
guess = noise;
GCR(src,guess);
subspace[b] = guess;
#endif
noise = subspace[b];
scale = std::pow(norm2(noise),-0.5);
noise=noise*scale;
}
DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<std::endl;
subspace[b] = noise;
}
}
//////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////
// World of possibilities here. But have tried quite a lot of experiments (250+ jobs run on Summit) // World of possibilities here. But have tried quite a lot of experiments (250+ jobs run on Summit)
// and this is the best I found // and this is the best I found
@ -160,14 +209,21 @@ public:
int b =0; int b =0;
{ {
ComplexD ip;
// Filter // Filter
Chebyshev<FineField> Cheb(lo,hi,orderfilter); Chebyshev<FineField> Cheb(lo,hi,orderfilter);
Cheb(hermop,noise,Mn); Cheb(hermop,noise,Mn);
// normalise // normalise
scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale; scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
subspace[b] = Mn; subspace[b] = Mn;
hermop.Op(Mn,tmp);
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl; hermop.Op(Mn,tmp);
ip= innerProduct(Mn,tmp);
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|Op|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
hermop.AdjOp(Mn,tmp);
ip = innerProduct(Mn,tmp);
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|AdjOp|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
b++; b++;
} }
@ -213,8 +269,18 @@ public:
Mn=*Tnp; Mn=*Tnp;
scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale; scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
subspace[b] = Mn; subspace[b] = Mn;
hermop.Op(Mn,tmp);
std::cout<<GridLogMessage << n<<" filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
ComplexD ip;
hermop.Op(Mn,tmp);
ip= innerProduct(Mn,tmp);
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|Op|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
hermop.AdjOp(Mn,tmp);
ip = innerProduct(Mn,tmp);
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|AdjOp|n> "<<norm2(tmp)<<" "<<ip<<std::endl;
b++; b++;
} }

51
tests/debug/Test_general_coarse_pvdagm.cc
View File

@ -154,6 +154,8 @@ public:
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl; // std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t; double t;
// Fine Smoother // Fine Smoother
// out = in;
out = Zero();
t=-usecond(); t=-usecond();
_PreSmoother(in,out); _PreSmoother(in,out);
t+=usecond(); t+=usecond();
@ -172,6 +174,7 @@ public:
// Coarse correction // Coarse correction
t=-usecond(); t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol); _CoarseSolve(Csrc,Csol);
//Csol=Zero(); //Csol=Zero();
t+=usecond(); t+=usecond();
@ -191,6 +194,8 @@ public:
// Fine Smoother // Fine Smoother
t=-usecond(); t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2); _PostSmoother(vec1,vec2);
t+=usecond(); t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl; std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
@ -215,7 +220,8 @@ int main (int argc, char ** argv)
// Construct a coarsened grid // Construct a coarsened grid
Coordinate clatt = GridDefaultLatt(); Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){ for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/4; clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
} }
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());; GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d); GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
@ -244,7 +250,7 @@ int main (int argc, char ** argv)
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5); DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
const int nbasis = 8; const int nbasis = 20;
const int cb = 0 ; const int cb = 0 ;
LatticeFermion prom(FGrid); LatticeFermion prom(FGrid);
@ -260,7 +266,25 @@ int main (int argc, char ** argv)
typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t; typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t; typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv); PVdagM_t PVdagM(Ddwf,Dpv);
ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
// ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
// ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ...
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ... with 24 vecs
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ... with 24 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ... with 16 vecs and 2^4 blocking, sloppier
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking, looser coarse
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64 ... with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); //
// Run power method on HOA?? // Run power method on HOA??
@ -269,6 +293,7 @@ int main (int argc, char ** argv)
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace; typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace AggregatesPD(Coarse5d,FGrid,cb); Subspace AggregatesPD(Coarse5d,FGrid,cb);
/*
AggregatesPD.CreateSubspaceChebyshev(RNG5, AggregatesPD.CreateSubspaceChebyshev(RNG5,
PVdagM, PVdagM,
nbasis, nbasis,
@ -278,6 +303,10 @@ int main (int argc, char ** argv)
200, 200,
200, 200,
0.0); 0.0);
*/
AggregatesPD.CreateSubspaceGCR(RNG5,
PVdagM,
nbasis);
LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD); LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
@ -334,12 +363,13 @@ int main (int argc, char ** argv)
/////////////////////////////////////// ///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl; std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve "<<std::endl; std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl; std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple; TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV); NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-8, 100, LinOpCoarse,simple,10,10); // PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
L2PGCR.Level(2); PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(3.0e-2, 100, LinOpCoarse,simple,10,10);
L2PGCR.Level(3);
c_res=Zero(); c_res=Zero();
L2PGCR(c_src,c_res); L2PGCR(c_src,c_res);
@ -347,11 +377,12 @@ int main (int argc, char ** argv)
// Fine grid smoother // Fine grid smoother
//////////////////////////////////////// ////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl; std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother "<<std::endl; std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl; std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine; TrivialPrecon<LatticeFermionD> simple_fine;
// NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM); // NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,10,ShiftedPVdagM,simple_fine,4,4); PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid); LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid); LatticeFermionD f_res(FGrid);
@ -364,12 +395,12 @@ int main (int argc, char ** argv)
TwoLevelMG TwoLevelPrecon(AggregatesPD, TwoLevelMG TwoLevelPrecon(AggregatesPD,
PVdagM, PVdagM,
SmootherGCR, simple_fine,
SmootherGCR, SmootherGCR,
LinOpCoarse, LinOpCoarse,
L2PGCR); L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,8,8); PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,16,16);
L1PGCR.Level(1); L1PGCR.Level(1);
f_res=Zero(); f_res=Zero();