portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2026-08-24 03:19:35 +01:00
Code Releases Activity

Verbose trasing

Browse Source
This commit is contained in:
Peter Boyle
2026-08-20 11:46:33 -04:00
parent 97c9178785
commit 3ed05dbbb1
1 changed files with 90 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files
examples/Example_pvdagm_mrhs_3level_DenseCoarseMatrix.cc
+90 -36
View File
@@ -198,6 +198,7 @@ public:
_Dense.ApplyBatch(split_in, split_out); _Dense.ApplyBatch(split_in, split_out);
for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0); for(int r=0;r<_nrhs;r++) InsertSliceFast(split_out[r], out, r, 0);
} else { } else {
GRID_TRACE("MrhsDensCCSolve::ApplyBatch6D");
_Dense.ApplyBatch6D(in, out, _nrhs); _Dense.ApplyBatch6D(in, out, _nrhs);
} }
} }
@@ -226,11 +227,12 @@ public:
void SetZeroGuess(int z){ ZeroGuess=z; } void SetZeroGuess(int z){ ZeroGuess=z; }
MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,MrhsLinearFunction<Field> &Prec,int _mmax,int _nstep) MrhsPGCRNonHermitian(RealD tol,Integer maxit,LinearOperatorBase<Field> &_Linop,MrhsLinearFunction<Field> &Prec,int _mmax,int _nstep)
: Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; } : Tolerance(tol),MaxIterations(maxit),Linop(_Linop),Preconditioner(Prec),mmax(_mmax),nstep(_nstep){ level=1; }
static RealD vnorm2(std::vector<Field> &x){ RealD s=0; for(auto &f:x) s+=norm2(f); return s; } static RealD vnorm2(std::vector<Field> &x){ GRID_TRACE("GCR-vnorm2"); RealD s=0; for(auto &f:x) s+=norm2(f); return s; }
static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &y){ ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; } static ComplexD vinnerProduct(std::vector<Field> &x,std::vector<Field> &y){ GRID_TRACE("GCR-vinnerProduct"); ComplexD s(0); for(int r=0;r<(int)x.size();r++) s+=innerProduct(x[r],y[r]); return s; }
static void vaxpy(std::vector<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); } static void vaxpy(std::vector<Field> &z,ComplexD a,std::vector<Field> &x,std::vector<Field> &y){ GRID_TRACE("GCR-vaxpy"); for(int r=0;r<(int)z.size();r++) axpy(z[r],a,x[r],y[r]); }
void vOp(std::vector<Field> &in,std::vector<Field> &out){ for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); } void vOp(std::vector<Field> &in,std::vector<Field> &out){ GRID_TRACE("GCR-vOp"); for(int r=0;r<(int)in.size();r++) Linop.Op(in[r],out[r]); }
void operator()(std::vector<Field> &src,std::vector<Field> &psi){ void operator()(std::vector<Field> &src,std::vector<Field> &psi){
GRID_TRACE((name+"MrhsPGCRNonHermitian").c_str());
RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid(); RealD cp,ssq,rsq; int nrhs=src.size(); GridBase *grid=src[0].Grid();
ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq; ssq=vnorm2(src); rsq=Tolerance*Tolerance*ssq;
std::vector<Field> r(nrhs,grid); std::vector<Field> r(nrhs,grid);
@@ -264,17 +266,32 @@ public:
p[0]=z; q[0]=Az; qq[0]=zAAz; cp=vnorm2(r); p[0]=z; q[0]=Az; qq[0]=zAAz; cp=vnorm2(r);
for(int k=0;k<nstep;k++){ for(int k=0;k<nstep;k++){
steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax; steps++; int kp=k+1, peri_k=k%mmax, peri_kp=kp%mmax;
rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k]; {
vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r); GRID_TRACE("GCR-project");
rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
}
// rq=vinnerProduct(q[peri_k],r); a=rq/qq[peri_k];
// vaxpy(psi,a,p[peri_k],psi); vaxpy(r,-a,q[peri_k],r); cp=vnorm2(r);
std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"] resid "<<cp<<" target "<<rsq<<std::endl; std::cout<<GridLogMessage<<std::string(level,'\t')<<" "<<name<<" MrhsPGCR step["<<steps<<"] resid "<<cp<<" target "<<rsq<<std::endl;
if((k==nstep-1)||(cp<rsq)) return cp; if((k==nstep-1)||(cp<rsq)) return cp;
Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az); Preconditioner(r,z); vOp(z,Az); zAAz=vnorm2(Az);
q[peri_kp]=Az; p[peri_kp]=z;
int northog=((kp)>(mmax-1))?(mmax-1):(kp); {
for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0); GRID_TRACE("GCR-orthog");
b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back]; q[peri_kp]=Az; p[peri_kp]=z;
vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); } int northog=((kp)>(mmax-1))?(mmax-1):(kp);
qq[peri_kp]=vnorm2(q[peri_kp]); for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0);
b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back];
vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); }
qq[peri_kp]=vnorm2(q[peri_kp]);
}
// q[peri_kp]=Az; p[peri_kp]=z;
// int northog=((kp)>(mmax-1))?(mmax-1):(kp);
// for(int back=0;back<northog;back++){ int peri_back=(k-back)%mmax; GRID_ASSERT((k-back)>=0);
// b=-real(vinnerProduct(q[peri_back],Az))/qq[peri_back];
// vaxpy(p[peri_kp],b,p[peri_back],p[peri_kp]); vaxpy(q[peri_kp],b,q[peri_back],q[peri_kp]); }
// qq[peri_kp]=vnorm2(q[peri_kp]);
} }
GRID_ASSERT(0); return cp; GRID_ASSERT(0); return cp;
} }
@@ -323,33 +340,45 @@ public:
CoarseCoarseField CCsol(_CoarseCoarseMrhs); CoarseCoarseField CCsol(_CoarseCoarseMrhs);
t=-usecond(); t=-usecond();
for(int r=0;r<nrhs;r++) ExtractSliceFast(csplit[r],vec1,r,0); {
_Projector.blockProject(csplit,ccsplit); GRID_TRACE("L2L3-Vcycle - Extract/blockProject/Insert");
for(int r=0;r<nrhs;r++) InsertSliceFast(ccsplit[r],CCsrc,r,0); for(int r=0;r<nrhs;r++) ExtractSliceFast(csplit[r],vec1,r,0);
_Projector.blockProject(csplit,ccsplit);
for(int r=0;r<nrhs;r++) InsertSliceFast(ccsplit[r],CCsrc,r,0);
}
t+=usecond(); t+=usecond();
std::cout<<GridLogMessage<<"L2->L3 restrict took "<<t/1000.0<<"ms"<<std::endl; std::cout<<GridLogMessage<<"L2->L3 restrict took "<<t/1000.0<<"ms"<<std::endl;
// L3 solve (dense mrhs GEMM, or PGCR) // L3 solve (dense mrhs GEMM, or PGCR)
t=-usecond(); t=-usecond();
CCsol=Zero(); {
_CoarseCoarseSolve(CCsrc,CCsol); GRID_TRACE("L2L3-Vcycle - CoarseCoarseSolve");
CCsol=Zero();
_CoarseCoarseSolve(CCsrc,CCsol);
}
t+=usecond(); t+=usecond();
std::cout<<GridLogMessage<<"L3 coarse-coarse solve took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl; std::cout<<GridLogMessage<<"L3 coarse-coarse solve took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
// prolong: unpack 6D cc -> blockPromote -> pack 6D coarse; add correction // prolong: unpack 6D cc -> blockPromote -> pack 6D coarse; add correction
t=-usecond(); t=-usecond();
for(int r=0;r<nrhs;r++) ExtractSliceFast(ccsplit[r],CCsol,r,0); {
_Projector.blockPromote(csplit,ccsplit); GRID_TRACE("L2L3-Vcycle - Ext/blockPromote/Ins");
for(int r=0;r<nrhs;r++) InsertSliceFast(csplit[r],vec1,r,0); for(int r=0;r<nrhs;r++) ExtractSliceFast(ccsplit[r],CCsol,r,0);
add(out,out,vec1); _Projector.blockPromote(csplit,ccsplit);
for(int r=0;r<nrhs;r++) InsertSliceFast(csplit[r],vec1,r,0);
add(out,out,vec1);
}
t+=usecond(); t+=usecond();
std::cout<<GridLogMessage<<"L2->L3 prolong took "<<t/1000.0<<"ms"<<std::endl; std::cout<<GridLogMessage<<"L2->L3 prolong took "<<t/1000.0<<"ms"<<std::endl;
// residual + coarse smoother (6D coarse) // residual + coarse smoother (6D coarse)
_CoarseOp.Op(out,vec1); sub(vec1,in,vec1); {
vec2=Zero(); GRID_TRACE("L2L3-Vcycle - Resid+CoarseSmooth");
_CoarseSmoother(vec1,vec2); _CoarseOp.Op(out,vec1); sub(vec1,in,vec1);
add(out,out,vec2); vec2=Zero();
_CoarseSmoother(vec1,vec2);
add(out,out,vec2);
}
} }
}; };
@@ -372,25 +401,50 @@ public:
_CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){} _CoarseGrid(CoarseGrid),_CoarseGridMrhs(CoarseGridMrhs){}
virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){ virtual void operator()(std::vector<FineField> &in, std::vector<FineField> &out){
int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t; int nrhs=in.size(); GridBase *fgrid=in[0].Grid(); double t;
std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid); std::vector<FineField> vec1(nrhs,fgrid),vec2(nrhs,fgrid);
for(int r=0;r<nrhs;r++) out[r]=in[r]; {
for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); } GRID_TRACE("L1L2-Vcycle - Resid");
for(int r=0;r<nrhs;r++) out[r]=in[r];
for(int r=0;r<nrhs;r++){
_FineOperator.Op(out[r],vec1[r]);
sub(vec1[r],in[r],vec1[r]);
}
}
std::vector<CoarseVector> Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid); std::vector<CoarseVector> Csrc_split(nrhs,_CoarseGrid), Csol_split(nrhs,_CoarseGrid);
CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs); CoarseVector CsrcMrhs(_CoarseGridMrhs), CsolMrhs(_CoarseGridMrhs);
t=-usecond(); t=-usecond();
_Projector.blockProject(vec1,Csrc_split); {
for(int r=0;r<nrhs;r++) InsertSliceFast(Csrc_split[r],CsrcMrhs,r,0); GRID_TRACE("L1L2-Vcycle - blockProject");
_Projector.blockProject(vec1,Csrc_split);
for(int r=0;r<nrhs;r++) InsertSliceFast(Csrc_split[r],CsrcMrhs,r,0);
}
t+=usecond(); std::cout<<GridLogMessage<<"Mrhs project+pack took "<<t/1000.0<<"ms"<<std::endl; t+=usecond(); std::cout<<GridLogMessage<<"Mrhs project+pack took "<<t/1000.0<<"ms"<<std::endl;
t=-usecond(); CsolMrhs=Zero(); _CoarseSolve(CsrcMrhs,CsolMrhs); t+=usecond(); t=-usecond();
{
GRID_TRACE("L1L2-Vcycle - CoarseSolve");
CsolMrhs=Zero();
_CoarseSolve(CsrcMrhs,CsolMrhs);
}
t+=usecond();
std::cout<<GridLogMessage<<"Mrhs coarse solve took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl; std::cout<<GridLogMessage<<"Mrhs coarse solve took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
t=-usecond(); t=-usecond();
for(int r=0;r<nrhs;r++) ExtractSliceFast(Csol_split[r],CsolMrhs,r,0); {
_Projector.blockPromote(vec1,Csol_split); GRID_TRACE("L1L2-Vcycle - ExtractSlice/blockPromote/add delta");
for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]); for(int r=0;r<nrhs;r++) ExtractSliceFast(Csol_split[r],CsolMrhs,r,0);
_Projector.blockPromote(vec1,Csol_split);
for(int r=0;r<nrhs;r++) add(out[r],out[r],vec1[r]);
}
t+=usecond(); std::cout<<GridLogMessage<<"Mrhs unpack+promote took "<<t/1000.0<<"ms"<<std::endl; t+=usecond(); std::cout<<GridLogMessage<<"Mrhs unpack+promote took "<<t/1000.0<<"ms"<<std::endl;
for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); } {
GRID_TRACE("L1L2-Vcycle - Residuals");
for(int r=0;r<nrhs;r++){ _FineOperator.Op(out[r],vec1[r]); sub(vec1[r],in[r],vec1[r]); }
}
t=-usecond(); t=-usecond();
for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); } {
GRID_TRACE("L1L2-Vcycle - Smoothers");
for(int r=0;r<nrhs;r++){ vec2[r]=Zero(); _PostSmoother(vec1[r],vec2[r]); add(out[r],out[r],vec2[r]); }
}
t+=usecond(); std::cout<<GridLogMessage<<"Mrhs post-smooth took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl; t+=usecond(); std::cout<<GridLogMessage<<"Mrhs post-smooth took "<<t/1000.0<<"ms ("<<t/1000.0/nrhs<<"ms/rhs)"<<std::endl;
} }
}; };
@@ -419,7 +473,7 @@ int main (int argc, char ** argv)
// Level 2 blocking: SUPERCOARSE default 8,4,3,6 -> CC [3,6,8,8], the dense floor. // Level 2 blocking: SUPERCOARSE default 8,4,3,6 -> CC [3,6,8,8], the dense floor.
Coordinate cclatt = clatt; Coordinate cclatt = clatt;
Coordinate Block2({8,4,3,6}); Coordinate Block2({4,4,3,6});
if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); } if ( getenv("BLOCK2") ){ GridCmdOptionIntVector(std::string(getenv("BLOCK2")),Block2); GRID_ASSERT(Block2.size()==4); }
for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; } for(int d=0;d<4;d++){ GRID_ASSERT(clatt[d]%Block2[d]==0); cclatt[d]=clatt[d]/Block2[d]; }
std::cout << GridLogMessage << "Block2 " << Block2 << " coarse-coarse lattice " << cclatt << std::endl; std::cout << GridLogMessage << "Block2 " << Block2 << " coarse-coarse lattice " << cclatt << std::endl;