1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00

Mixed precision for Laplace. Main program with Metric

This commit is contained in:
Chulwoo Jung 2024-02-08 17:13:10 -05:00
parent 026eb8a695
commit 9ad6836b0f
4 changed files with 148 additions and 145 deletions

View File

@ -148,116 +148,117 @@ public:
GaugeField& der , std::vector< std::vector<GaugeLinkField> >& prev_solns ) {
// get rid of this please
std::cout<<GridLogMessage << "LaplaceStart " <<std::endl;
RealD fac = - 1. / (double(4 * Nd)) ;
RealD coef=0.5;
LapStencil.GaugeImport(Usav);
LapStencilF.GaugeImport(UsavF);
for (int nu=0;nu<Nd;nu++){
GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
GaugeLinkField LMinvMom(left.Grid());
GaugeLinkField GMom(left.Grid());
GaugeLinkField LMinvGMom(left.Grid());
GaugeLinkField AGMom(left.Grid());
GaugeLinkField MinvAGMom(left.Grid());
GaugeLinkField LMinvAGMom(left.Grid());
GaugeLinkField AMinvMom(left.Grid());
GaugeLinkField LMinvAMom(left.Grid());
GaugeLinkField temp(left.Grid());
GaugeLinkField temp2(left.Grid());
std::vector<GaugeLinkField> MinvMom(par.order,left.Grid());
GaugeLinkField MinvGMom(left.Grid());
GaugeLinkField Gtemp(left.Grid());
GaugeLinkField Gtemp2(left.Grid());
ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000,false);
// ConjugateGradient<GaugeFieldF> CG_f(par.tolerance,10000,false);
LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,GaugeLinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
GMom = par.offset * right_nu;
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
#if USE_CHRONO
MinvMom[i] = Forecast(QuadOp, right_nu, prev_solns[nu]);
#endif
#ifndef MIXED_CG
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(right_nu,MinvMom[i]);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(MinvMom[i]);
#endif
for (int nu=0;nu<Nd;nu++){
GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
GaugeLinkField LMinvMom(left.Grid());
GMom += par.a0[i]*MinvMom[i];
LapStencil.M(MinvMom[i],Gtemp2);
GMom += par.a1[i]*fac*Gtemp2;
GaugeLinkField GMom(left.Grid());
GaugeLinkField LMinvGMom(left.Grid());
GaugeLinkField AGMom(left.Grid());
GaugeLinkField MinvAGMom(left.Grid());
GaugeLinkField LMinvAGMom(left.Grid());
GaugeLinkField AMinvMom(left.Grid());
GaugeLinkField LMinvAMom(left.Grid());
GaugeLinkField temp(left.Grid());
GaugeLinkField temp2(left.Grid());
std::vector<GaugeLinkField> MinvMom(par.order,left.Grid());
GaugeLinkField MinvGMom(left.Grid());
GaugeLinkField Gtemp(left.Grid());
GaugeLinkField Gtemp2(left.Grid());
ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000,false);
// ConjugateGradient<GaugeFieldF> CG_f(par.tolerance,10000,false);
LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,GaugeLinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
GMom = par.offset * right_nu;
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
#if USE_CHRONO
MinvMom[i] = Forecast(QuadOp, right_nu, prev_solns[nu]);
#endif
#ifndef MIXED_CG
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(right_nu,MinvMom[i]);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(MinvMom[i]);
#endif
GMom += par.a0[i]*MinvMom[i];
LapStencil.M(MinvMom[i],Gtemp2);
GMom += par.a1[i]*fac*Gtemp2;
}
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
MinvGMom = Forecast(QuadOp, GMom, prev_solns[nu]);
#ifndef MIXED_CG
CG(QuadOp,GMom,MinvGMom);
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(GMom,MinvGMom);
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
// Laplacian.M(MinvGMom, LMinvGMom);
MixedCG(right_nu,MinvMom[i]);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(MinvGMom);
#endif
LapStencil.M(MinvMom[i], Gtemp2); LMinvMom=fac*Gtemp2;
AMinvMom = par.a1[i]*LMinvMom;
AMinvMom += par.a0[i]*MinvMom[i];
LapStencil.M(AMinvMom, Gtemp2); LMinvAMom=fac*Gtemp2;
LapStencil.M(MinvGMom, Gtemp2); temp=fac*Gtemp2;
MinvAGMom = par.a1[i]*temp;
MinvAGMom += par.a0[i]*MinvGMom;
LapStencil.M(MinvAGMom, Gtemp2); LMinvAGMom=fac*Gtemp2;
GaugeField tempDer(left.Grid());
std::cout<<GridLogMessage << "force contraction "<< i <<std::endl;
// roctxRangePushA("RMHMC force contraction");
MDerivLink(GMom,MinvMom[i],tempDer); der += coef*2*par.a1[i]*tempDer;
MDerivLink(left_nu,MinvGMom,tempDer); der += coef*2*par.a1[i]*tempDer;
MDerivLink(LMinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(LMinvAMom,MinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(MinvAGMom,LMinvMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(AMinvMom,LMinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(MinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b1[i]*tempDer;
MDerivLink(AMinvMom,MinvGMom,tempDer); der += coef*-2.*par.b1[i]*tempDer;
std::cout<<GridLogMessage << "coef = force contraction "<< i << "done "<< coef <<std::endl;
// roctxRangePop();
}
}
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
MinvGMom = Forecast(QuadOp, GMom, prev_solns[nu]);
#ifndef MIXED_CG
CG(QuadOp,GMom,MinvGMom);
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
CG(QuadOp,right_nu,MinvMom[i]);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(GMom,MinvGMom);
LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
// Laplacian.M(MinvGMom, LMinvGMom);
MixedCG(right_nu,MinvMom[i]);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(MinvGMom);
#endif
LapStencil.M(MinvMom[i], Gtemp2); LMinvMom=fac*Gtemp2;
AMinvMom = par.a1[i]*LMinvMom;
AMinvMom += par.a0[i]*MinvMom[i];
LapStencil.M(AMinvMom, Gtemp2); LMinvAMom=fac*Gtemp2;
LapStencil.M(MinvGMom, Gtemp2); temp=fac*Gtemp2;
MinvAGMom = par.a1[i]*temp;
MinvAGMom += par.a0[i]*MinvGMom;
LapStencil.M(MinvAGMom, Gtemp2); LMinvAGMom=fac*Gtemp2;
GaugeField tempDer(left.Grid());
std::cout<<GridLogMessage << "force contraction "<< i <<std::endl;
// roctxRangePushA("RMHMC force contraction");
MDerivLink(GMom,MinvMom[i],tempDer); der += coef*2*par.a1[i]*tempDer;
MDerivLink(left_nu,MinvGMom,tempDer); der += coef*2*par.a1[i]*tempDer;
MDerivLink(LMinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(LMinvAMom,MinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(MinvAGMom,LMinvMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(AMinvMom,LMinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
MDerivLink(MinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b1[i]*tempDer;
MDerivLink(AMinvMom,MinvGMom,tempDer); der += coef*-2.*par.b1[i]*tempDer;
std::cout<<GridLogMessage << "coef = force contraction "<< i << "done "<< coef <<std::endl;
// roctxRangePop();
}
}
std::cout<<GridLogMessage << "LaplaceEnd " <<std::endl;
// exit(-42);
}
@ -283,45 +284,47 @@ for (int nu=0;nu<Nd;nu++){
void MSquareRootInt(LaplacianRatParams &par, GaugeField& P, std::vector< std::vector<GaugeLinkField> > & prev_solns ){
std::cout<<GridLogMessage << "LaplaceStart " <<std::endl;
RealD fac = -1. / (double(4 * Nd));
LapStencil.GaugeImport(Usav);
LapStencilF.GaugeImport(UsavF);
for(int nu=0; nu<Nd;nu++){
GaugeLinkField P_nu = PeekIndex<LorentzIndex>(P, nu);
GaugeLinkField Gp(P.Grid());
Gp = par.offset * P_nu;
ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000);
// ConjugateGradient<GaugeLinkFieldF> CG_f(1.0e-8,10000);
ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
GaugeLinkField Gtemp(P.Grid());
GaugeLinkField Gtemp2(P.Grid());
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
Gtemp = Forecast(QuadOp, P_nu, prev_solns[nu]);
#ifndef MIXED_CG
CG(QuadOp,P_nu,Gtemp);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(P_nu,Gtemp);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(Gtemp);
#endif
Gp += par.a0[i]*Gtemp;
LapStencil.M(Gtemp,Gtemp2);
Gp += par.a1[i]*fac*Gtemp2;
for(int nu=0; nu<Nd;nu++){
GaugeLinkField P_nu = PeekIndex<LorentzIndex>(P, nu);
GaugeLinkField Gp(P.Grid());
Gp = par.offset * P_nu;
ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000);
// ConjugateGradient<GaugeLinkFieldF> CG_f(1.0e-8,10000);
ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
GaugeLinkField Gtemp(P.Grid());
GaugeLinkField Gtemp2(P.Grid());
for(int i =0;i<par.order;i++){
QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
Gtemp = Forecast(QuadOp, P_nu, prev_solns[nu]);
#ifndef MIXED_CG
CG(QuadOp,P_nu,Gtemp);
#else
QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
// QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
MixedCG.InnerTolerance=par.tolerance;
MixedCG(P_nu,Gtemp);
#endif
#if USE_CHRONO
prev_solns[nu].push_back(Gtemp);
#endif
Gp += par.a0[i]*Gtemp;
LapStencil.M(Gtemp,Gtemp2);
Gp += par.a1[i]*fac*Gtemp2;
}
PokeIndex<LorentzIndex>(P, Gp, nu);
}
PokeIndex<LorentzIndex>(P, Gp, nu);
}
std::cout<<GridLogMessage << "LaplaceEnd " <<std::endl;
}
void MSquareRoot(GaugeField& P){

View File

@ -36,7 +36,7 @@ directory
#endif
// second level EOFA
#undef EOFA_H
#define USE_OBC
#undef USE_OBC
#define DO_IMPLICIT
NAMESPACE_BEGIN(Grid);
@ -203,9 +203,9 @@ int main(int argc, char **argv) {
HMCparams.MD.name =std::string("ImplicitMinimumNorm2");
#else
// typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
HMCparams.MD.name =std::string("MinimumNorm2");
typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
HMCparams.MD.name =std::string("ForceGradient");
#endif
std::cout << GridLogMessage<< HMCparams <<std::endl;
@ -344,7 +344,7 @@ int main(int argc, char **argv) {
ConjugateGradient<FermionField> ActionCG(ActionStoppingCondition,MaxCGIterations);
ConjugateGradient<FermionField> DerivativeCG(DerivativeStoppingCondition,MaxCGIterations);
#ifdef MIXED_PRECISION
const int MX_inner = 5000;
const int MX_inner = 50000;
// Mixed precision EOFA
LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L);

View File

@ -83,7 +83,7 @@ int main(int argc, char **argv)
// need wrappers of the fermionic classes
// that have a complex construction
// standard
RealD beta = 6.4 ;
RealD beta = 6.6 ;
#if 0
WilsonGaugeActionR Waction(beta);

View File

@ -30,7 +30,7 @@ directory
/* END LEGAL */
#include <Grid/Grid.h>
#define USE_OBC
#undef USE_OBC
#define DO_IMPLICIT
@ -138,7 +138,7 @@ int main(int argc, char **argv)
// that have a complex construction
// standard
RealD beta = 6.4;
RealD beta = 6.6;
std::cout << "Wilson Gauge beta= " <<beta <<std::endl;
#ifndef USE_OBC
WilsonGaugeActionR Waction(beta);