Integrators for the HMC. More...

#include <stdbool.h>
#include <random.h>
#include <sizes.h>

Include dependency graph for integrate.h:

This graph shows which files directly or indirectly include this file:

Functions
int	Gauge_Update (const double d, double pp, Complex u11t, Complex u12t, Complex_f u11t_f, Complex_f *u12t_f)
	Gauge update for the integration step of the HMC.

int	Momentum_Update (const double d, const double dSdpi, double pp)
	Wrapper for the momentum update during the integration step of the HMC.

int	Leapfrog (Complex u11t, Complex u12t, Complex_f u11t_f, Complex_f u12t_f, Complex X0, Complex X1, Complex Phi, double dk4m, double dk4p, float dk4m_f, float dk4p_f, double dSdpi, double pp, int iu, int id, Complex gamval, Complex_f gamval_f, int gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double ancg, int itot, float proby)
	Leapfrog integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.

int	OMF2 (Complex u11t, Complex u12t, Complex_f u11t_f, Complex_f u12t_f, Complex X0, Complex X1, Complex Phi, double dk4m, double dk4p, float dk4m_f, float dk4p_f, double dSdpi, double pp, int iu, int id, Complex gamval, Complex_f gamval_f, int gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double ancg, int itot, float proby)
	OMF integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.

int	OMF4 (Complex u11t, Complex u12t, Complex_f u11t_f, Complex_f u12t_f, Complex X0, Complex X1, Complex Phi, double dk4m, double dk4p, float dk4m_f, float dk4p_f, double dSdpi, double pp, int iu, int id, Complex gamval, Complex_f gamval_f, int gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double ancg, int itot, float proby)

Detailed Description

Integrators for the HMC.

Definition in file integrate.h.

Function Documentation

◆ Gauge_Update()

int Gauge_Update	(	const double	d,
		double *	pp,
		Complex *	u11t,
		Complex *	u12t,
		Complex_f *	u11t_f,
		Complex_f *	u12t_f )

Gauge update for the integration step of the HMC.

Parameters

dt	Gauge step size
pp	Momentum field
u11t,u12t	Double precision gauge fields
u11t_f,u12t_f	Single precision gauge fields

Returns: Zero on success, integer error code otherwise.

Definition at line 3 of file integrate.c.

                                                                                                              {
   /*
    * @brief Generates new trial fields
    *
    * @see cuGauge_Update (CUDA Wrapper)
    * 
    * @param   d:    Half lattice spacing
    * @param   pp:      Momentum field
    * @param   u11t:    First colour field
    * @param   u12t:    Second colour field
    *
    * @returns Zero on success, integer error code otherwise
    */
   char *funcname = "Gauge_Update"; 
#ifdef __NVCC__
   cuGauge_Update(d,pp,u11t,u12t,dimGrid,dimBlock);
#else
#pragma omp parallel for simd collapse(2) aligned(pp,u11t,u12t:AVX) 
   for(int i=0;i<kvol;i++)
      for(int mu = 0; mu<ndim; mu++){
         /*
          * Sticking to what was in the FORTRAN for variable names.
          * CCC for cosine SSS for sine AAA for...
          * Re-exponentiating the force field. Can be done analytically in SU(2)
          * using sine and cosine which is nice
          */
 
         double AAA = d*sqrt(pp[i*nadj*ndim+mu]*pp[i*nadj*ndim+mu]\
               +pp[(i*nadj+1)*ndim+mu]*pp[(i*nadj+1)*ndim+mu]\
               +pp[(i*nadj+2)*ndim+mu]*pp[(i*nadj+2)*ndim+mu]);
         double CCC = cos(AAA);
         double SSS = d*sin(AAA)/AAA;
         Complex a11 = CCC+I*SSS*pp[(i*nadj+2)*ndim+mu];
         Complex a12 = pp[(i*nadj+1)*ndim+mu]*SSS + I*SSS*pp[i*nadj*ndim+mu];
         //b11 and b12 are u11t and u12t terms, so we'll use u12t directly
         //but use b11 for u11t to prevent RAW dependency
         complex b11 = u11t[i*ndim+mu];
         u11t[i*ndim+mu] = a11*b11-a12*conj(u12t[i*ndim+mu]);
         u12t[i*ndim+mu] = a11*u12t[i*ndim+mu]+a12*conj(b11);
      }
#endif
   Reunitarise(u11t,u12t);
   //Get trial fields from accelerator for halo exchange
   Trial_Exchange(u11t,u12t,u11t_f,u12t_f);
   return 0;
}

References Complex, Gauge_Update(), kvol, nadj, ndim, Reunitarise(), and Trial_Exchange().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Leapfrog()

int Leapfrog	(	Complex *	u11t,
		Complex *	u12t,
		Complex_f *	u11t_f,
		Complex_f *	u12t_f,
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk4m,
		double *	dk4p,
		float *	dk4m_f,
		float *	dk4p_f,
		double *	dSdpi,
		double *	pp,
		int *	iu,
		int *	id,
		Complex *	gamval,
		Complex_f *	gamval_f,
		int *	gamin,
		Complex	jqq,
		float	beta,
		float	akappa,
		int	stepl,
		float	dt,
		double *	ancg,
		int *	itot,
		float	proby )

Leapfrog integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.

Parameters

u11t,u12t	Double precision colour fields
u11t_f,u12t_f	Single precision colour fields
X0	Up/down partitioned pseudofermion field
X1	Holder for the partitioned fermion field, then the conjugate gradient output
Phi	Pseudofermion field
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
dk4m_f	\(\left(1+\gamma_0\right)e^{-\mu}\) float
dk4p_f	\(\left(1-\gamma_0\right)e^\mu\) float
dSdpi	The force
pp	Momentum field
iu,id	Lattice indices
gamin	Gamma indices
gamval	Double precision gamma matrices rescaled by kappa
gamval_f	Single precision gamma matrices rescaled by kappa
jqq	Diquark source
akappa	Hopping parameter
beta	Inverse gauge coupling
stepl	Steps per trajectory
dt	Step size
ancg	Counter for average conjugate gradient iterations
itot	Total average conjugate gradient iterations
proby	Termination probability for random trajectory length

Returns: Zero on success, integer error code otherwise

Definition at line 60 of file integrate.c.

{
   //This was originally in the half-step of the FORTRAN code, but it makes more sense to declare
   //it outside the loop. Since it's always being subtracted we'll define it as negative
   const double d =-dt*0.5;
   //Half step forward for p
   //=======================
#ifdef _DEBUG
   printf("Evaluating force on rank %i\n", rank);
#endif
   Force(dSdpi, 1, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
         dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
   Momentum_Update(d,dSdpi,pp);
   //Main loop for classical time evolution
   //======================================
   bool end_traj=false; int step =1;
   // for(int step = 1; step<=stepmax; step++){
   do{
#ifdef _DEBUG
      if(!rank)
         printf("Step: %d\n",  step);
#endif
      //The FORTRAN redefines d=dt here, which makes sense if you have a limited line length.
      //I'll stick to using dt though.
      //step (i) st(t+dt)=st(t)+p(t+dt/2)*dt;
      //Note that we are moving from kernel to kernel within the default streams so don't need a Device_Sync here
      Gauge_Update(dt,pp,u11t,u12t,u11t_f,u12t_f);
 
      //p(t+3et/2)=p(t+dt/2)-dSds(t+dt)*dt
      // Force(dSdpi, 0, rescgg);
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
 
   // if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
      if(step==stepl){
         //Final trajectory has a half momentum step
         Momentum_Update(d,dSdpi,pp);
         *itot+=step;
         *ancg/=step;
         end_traj=true;
         break;
      }
      else{
         //Otherwise, there's a half step at the end and start of each trajectory, so we combine them into one full step.
         Momentum_Update(-dt,dSdpi,pp);
         step++;
      }
   }while(!end_traj);
   return 0;
}

References Force(), Gauge_Update(), Leapfrog(), Momentum_Update(), rank, and rescgg.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Momentum_Update()

int Momentum_Update	(	const double	d,
		const double *	dSdpi,
		double *	pp )

inline

Wrapper for the momentum update during the integration step of the HMC.

Parameters

dt	Gauge step size
pp	Momentum field
dSdpi	Force field

Returns: Zero on success, integer error code otherwise.

Definition at line 49 of file integrate.c.

{
#ifdef __NVCC__
   cublasDaxpy(cublas_handle,kmom, &d, dSdpi, 1, pp, 1);
#elif defined USE_BLAS
   cblas_daxpy(kmom, d, dSdpi, 1, pp, 1);
#else
   for(int i=0;i<kmom;i++)
      pp[i]+=d*dSdpi[i];
#endif
}

References kmom, and Momentum_Update().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ OMF2()

int OMF2	(	Complex *	u11t,
		Complex *	u12t,
		Complex_f *	u11t_f,
		Complex_f *	u12t_f,
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk4m,
		double *	dk4p,
		float *	dk4m_f,
		float *	dk4p_f,
		double *	dSdpi,
		double *	pp,
		int *	iu,
		int *	id,
		Complex *	gamval,
		Complex_f *	gamval_f,
		int *	gamin,
		Complex	jqq,
		float	beta,
		float	akappa,
		int	stepl,
		float	dt,
		double *	ancg,
		int *	itot,
		float	proby )

OMF integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.

Parameters

u11t,u12t	Double precision colour fields
u11t_f,u12t_f	Single precision colour fields
X0	Up/down partitioned pseudofermion field
X1	Holder for the partitioned fermion field, then the conjugate gradient output
Phi	Pseudofermion field
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
dk4m_f	\(\left(1+\gamma_0\right)e^{-\mu}\) float
dk4p_f	\(\left(1-\gamma_0\right)e^\mu\) float
dSdpi	The force
pp	Momentum field
iu,id	Lattice indices
gamin	Gamma indices
gamval	Double precision gamma matrices rescaled by kappa
gamval_f	Single precision gamma matrices rescaled by kappa
jqq	Diquark source
akappa	Hopping parameter
beta	Inverse gauge coupling
stepl	Steps per trajectory
dt	Step size
ancg	Counter for average conjugate gradient iterations
itot	Total average conjugate gradient iterations
proby	Termination probability for random trajectory length

Returns: Zero on success, integer error code otherwise

Definition at line 113 of file integrate.c.

{
   const double lambda=0.5-(pow(2.0*sqrt(326.0)+36.0,1.0/3.0)/12.0)+1.0/(6*pow(2.0*sqrt(326.0) + 36.0,1.0/3.0));
   //const double lambda=1.0/6.0;
   // const double lambda=0.5;
 
   //Gauge update by half dt
   const double dU = dt*0.5;
 
   //Momentum updates by lambda, 2lambda and (1-2lambda) in the middle
   const double dp= -lambda*dt;
   const double dp2= 2.0*dp;
   const double dpm= -(1.0-2.0*lambda)*dt;
   //Initial step forward for p
   //=======================
#ifdef _DEBUG
   printf("Evaluating force on rank %i\n", rank);
#endif
   Force(dSdpi, 1, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
         dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
   //Initial momentum update
   Momentum_Update(dp,dSdpi,pp);
 
   //Main loop for classical time evolution
   //======================================
   bool end_traj=false; int step =1;
   // for(int step = 1; step<=stepmax; step++){
   do{
#ifdef _DEBUG
      if(!rank)
         printf("Step: %d\n", step);
#endif
      //First gauge update
      Gauge_Update(dU,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for middle momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
      //Now do the middle momentum update
      Momentum_Update(dpm,dSdpi,pp);
 
      //Second gauge update
      Gauge_Update(dU,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for second momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
 
      //if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
      if(step==stepl){
         //Final momentum step
         Momentum_Update(dp,dSdpi,pp);
         *itot+=step;
         //Two force terms, so an extra factor of two in the average?
         //Or leave it as it was, to get the average CG iterations per trajectory rather than force
         *ancg/=step;
         end_traj=true;
         break;
      }
      else{
         //Since we apply the momentum at the start and end of a step we instead apply a double step here
         Momentum_Update(dp2,dSdpi,pp);
         step++;
      }
   }while(!end_traj);
   return 0;
}

References Force(), Gauge_Update(), Momentum_Update(), OMF2(), rank, and rescgg.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ OMF4()

int OMF4	(	Complex *	u11t,
		Complex *	u12t,
		Complex_f *	u11t_f,
		Complex_f *	u12t_f,
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk4m,
		double *	dk4p,
		float *	dk4m_f,
		float *	dk4p_f,
		double *	dSdpi,
		double *	pp,
		int *	iu,
		int *	id,
		Complex *	gamval,
		Complex_f *	gamval_f,
		int *	gamin,
		Complex	jqq,
		float	beta,
		float	akappa,
		int	stepl,
		float	dt,
		double *	ancg,
		int *	itot,
		float	proby )

Momentum updates

Outer updates depend on r1. We have two of these, doubled for between full steps

Middle updates depend on r3

Inner updates depend on r1 and r3

Gauge updates. These depend on r2 and r4

Outer gauge update depends on r2

Middle gauge update depends on r4

Inner gauge update depends on r2 and r4

Definition at line 183 of file integrate.c.

{
   //These values were lifted from openqcd-fastsum, and should probably be tuned for QC2D. They also probably never
   //will be...
   const double r1 = 0.08398315262876693;
   const double r2 = 0.2539785108410595;
   const double r3 = 0.6822365335719091;
   const double r4 = -0.03230286765269967;
   const double dpO= -r1*dt;
   const double dpO2= 2*dpO;
   const double dpM= -r3*dt;
   const double dpI= -(0.5-r1-r3)*dt;
 
   const double duO = dt*r2;
   const double duM = dt*r4;
   const double duI = dt*(1-2*(r2+r4));
 
   //Initial step forward for p
   //=======================
#ifdef _DEBUG
   printf("Evaluating force on rank %i\n", rank);
#endif
   Force(dSdpi, 1, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
         dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
   Momentum_Update(dpO,dSdpi,pp);
 
   //Main loop for classical time evolution
   //======================================
   bool end_traj=false; int step =1;
   // for(int step = 1; step<=stepmax; step++){
   do{
#ifdef _DEBUG
      if(!rank)
         printf("Step: %d\n", step);
#endif
      //First outer gauge update
      Gauge_Update(duO,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for first middle momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
      //Now do the first middle momentum update
      Momentum_Update(dpM,dSdpi,pp);
 
      //First middle gauge update
      Gauge_Update(duM,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for first inner momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
      //Now do the first inner momentum update
      Momentum_Update(dpI,dSdpi,pp);
 
      //Inner gauge update
      Gauge_Update(duI,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for second inner momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
      //Now do the second inner momentum update
      Momentum_Update(dpI,dSdpi,pp);
 
      //Second middle gauge update
      Gauge_Update(duM,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for second middle momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
      //Now do the second middle momentum update
      Momentum_Update(dpM,dSdpi,pp);
 
      //Second outer gauge update
      Gauge_Update(duO,pp,u11t,u12t,u11t_f,u12t_f);
 
      //Calculate force for outer momentum update
      Force(dSdpi, 0, rescgg,X0,X1,Phi,u11t,u12t,u11t_f,u12t_f,iu,id,gamval,gamval_f,gamin,dk4m,dk4p,\
            dk4m_f,dk4p_f,jqq,akappa,beta,ancg);
 
      //Outer momentum update depends on if we've finished the trajectory
      //if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
      if(step==stepl){
         //Final momentum step
         Momentum_Update(dpO,dSdpi,pp);
         *itot+=step;
 
         //Four force terms, so an extra factor of four in the average?
         //Or leave it as it was, to get the average CG iterations per trajectory rather than force
         *ancg/=step;
         end_traj=true;
         break;
      }
      else{
         //Since we apply the momentum at the start and end of a step we instead apply a double step here
         Momentum_Update(dpO2,dSdpi,pp);
         step++;
      }
   }while(!end_traj);
   return 0;
}

References Force(), Gauge_Update(), Momentum_Update(), OMF4(), rank, and rescgg.

Here is the call graph for this function:

Here is the caller graph for this function:

Functions

Detailed Description

Function Documentation

◆ Gauge_Update()

◆ Leapfrog()

◆ Momentum_Update()

◆ OMF2()

◆ OMF4()