Molecular dynamics

Collaboration diagram for Molecular dynamics:

Topics
	Clover Force

Functions
int	Gauge_Update (const double d, double *pp, Complex *ut[2], Complex_f *ut_f[2])
	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 *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, const Complex jqq, const float beta, const float akappa, const float c_sw, const int stepl, const float dt, double *ancg, int *itot, const 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 *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, const Complex jqq, const float beta, const float akappa, const float c_sw, const int stepl, const float dt, double *ancg, int *itot, const float proby)
	OMF second order five step integrator.
int	OMF4 (Complex *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, const Complex jqq, const float beta, const float akappa, const float c_sw, const int stepl, const float dt, double *ancg, int *itot, const float proby)
	OMF fourth order eleven step integrator.
void	cuGauge_Update (const double d, double *pp, Complex *ut[2], dim3 dimGrid, dim3 dimBlock)
	CUDA wrapper for the gauge update during the integration step of the HMC.
void	Force_s (double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], unsigned int *iu, const unsigned short gamin[16], const float akappa, const unsigned short mu)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
void	Force_t (double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk[2], unsigned int *iu, const unsigned short gamin[16], float akappa)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
int	Force (double *dSdpi, const bool iflag, double res1, Complex *X0, Complex *X1, Complex *Phi, Complex *ut[2], Complex_f *ut_f[2], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], const Complex_f jqq, const float akappa, const float beta, const float c_sw, double *ancg)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
int	Gauge_force (double *dSdpi, Complex_f *ut[2], unsigned int *iu, unsigned int *id, float beta)
	Calculates the gauge force due to the Wilson Action at each intermediate time.
int	Congradq (int na, double res, Complex *X1, Complex *r, Complex *ud[2], Complex_f *ut[2], Complex_f *clover_f[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float c_sw, int *itercg)
	Matrix Inversion via Conjugate Gradient (up/down flavour partitioning). Solves \((M^\dagger)Mx=\Phi\) Implements up/down partitioning The matrix multiplication step is done at mixed precision, while the update is done at double.
int	Congradp (int na, double res, Complex *Phi, Complex *xi, Complex *ud[2], Complex_f *ut[2], Complex_f *clover_f[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float c_sw, int *itercg)
	Matrix Inversion via Conjugate Gradient (no up/down flavour partitioning). Solves \((M^\dagger)Mx=\Phi\) The matrix multiplication step is done at single precision, while the update is done at double.
void	cuGauge_force (Complex_f *ut[2], double *dSdpi, float beta, unsigned int *iu, unsigned int *id, dim3 dimGrid, dim3 dimBlock)
	Calculate the gauge contribution to the force.
void	cuForce (double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk[2], unsigned int *iu, const unsigned short gamin[16], float akappa, dim3 dimGrid, dim3 dimBlock)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
__global__ void	Kernels::Plus_staple (const int mu, const int nu, unsigned int *iu, Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11t, Complex_f *u12t)
	Calculates the staple in the positive \(\mu\) direction.
__global__ void	Kernels::Minus_staple (const int mu, const int nu, unsigned int *iu, unsigned int *id, Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11sh, Complex_f *u12sh, Complex_f *u11t, Complex_f *u12t)
	Calculates the staple in the negative \(\mu\) direction.
__global__ void	Kernels::cuGaugeForce (int mu, Complex_f *Sigma11, Complex_f *Sigma12, double *dSdpi, Complex_f *u11t, Complex_f *u12t, float beta)
	Calculates the gauge force due to the Wilson Action at each intermediate time.
__global__ void	Kernels::cuForce_s (double *dSdpi, Complex_f *u11t, Complex_f *u12t, Complex_f *X1, Complex_f *X2, Complex_f gamval[20], unsigned int *iu, const unsigned short gamin[16], float akappa, const unsigned short mu)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
__global__ void	Kernels::cuForce_t (double *dSdpi, Complex_f *u11t, Complex_f *u12t, Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk4m, float *dk4p, unsigned int *iu, const unsigned short gamin[16], float akappa)
	Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.
__global__ void	Kernels::cuGauge_Update (const double d, double *pp, Complex *u11t, Complex *u12t, int mu)
	Gauge update for the integration step of the HMC.

Detailed Description

Function Documentation

Congradp()

int Congradp	(	int	na,
		double	res,
		Complex *	Phi,
		Complex *	xi,
		Complex *	ud[2],
		Complex_f *	ut[2],
		Complex_f *	clover_f[nc],
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		double *	dk[2],
		float *	dk_f[2],
		Complex_f	jqq,
		float	akappa,
		float	c_sw,
		int *	itercg )

Matrix Inversion via Conjugate Gradient (no up/down flavour partitioning). Solves \((M^\dagger)Mx=\Phi\) The matrix multiplication step is done at single precision, while the update is done at double.

Parameters

[in]	na	Flavour index
[in]	res	Limit for conjugate gradient
[in]	Phi	Pseudofermion field.
[in,out]	xi	Returned as \((M^\dagger M)^{-1} \Phi\)
[in]	ut,ud	Double/float Gauge fields
[in]	iu,id	Upper/Lower halo indices
[in]	gamval,gamval_f	double float Gamma matrices rescaled by kappa
[in]	gamin	Dirac indices
[in]	clover_f	Array of clover fields
[in]	sigval,sigval_f	Double/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	dk,dk_f	Double/float \(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]	jqq	Diquark source
[in]	akappa	Hopping Parameter
[in]	c_sw	Clover coefficient.
[in]	itercg	Counts the iterations of the conjugate gradient

Returns

0 on success, integer error code otherwise

Postcondition

Contents of Phi overwritten

TODO: Set this as an argument How much does the residue have to shrink by before we do a double precision update

\(x2 = (M^\dagger M)p \)

\(\alpha_d= p* (M^\dagger M)p\)

\(alpha=\frac{\alpha_n} {\alpha_d}=\frac{r\cdot r}{p(M^\dagger M)p}\)

\(x+\alpha p\)

Definition at line 735 of file congrad.c.

                                                        {
   const char funcname[] = "Congradp";
   //Return value
   int ret_val=0;
   const double resid = res*res;
#ifdef _DEBUGCG
#warning "CG Debugging"
   char *endline = "\n";
#else
   char *endline = "\r";
#endif

   const float d_prec=1.0f/128.0f;
 
   //These were evaluated only in the first loop of niterx so we'll just do it outside of the loop.
   alignas(8) double alphan=1.0;
   //These alpha and beta terms should be double, but that causes issues with BLAS. Instead we declare
   //them Complex and work with the real part (especially for \alpha_d)
   //Give initial values Will be overwritten if niterx>0
   alignas(8) double betad = 1.0; alignas(8) Complex_f alphad=0; alignas(16)Complex alpha = 1; //Alignment needed for cuBLAS
 
   Complex_f *p_f, *r_f, *x1_f, *x2_f, *xi_f;
   Complex *p, *r, *x1, *x2, *clover[nc];
   P_allocate_f(&p_f,&r_f,&x1_f,&x2_f,&xi_f);
   P_allocate(&p, &r, &x1, &x2,clover);
 
   ComplexConvert(r_f,Phi+na*kferm,kferm,true,1);
   ComplexConvert(xi_f,xi,kferm,true,1);
   if(c_sw){
      ComplexConvert(clover_f[0],clover[0],6*kvol,false,1);
      ComplexConvert(clover_f[1],clover[1],6*kvol,false,1);
   }
   //Instead of copying element-wise in a loop, use memcpy.
#ifdef USE_GPU
   //Get r in single precision
   //Get xi  in single precision
   cudaMemcpy(r,Phi+na*kferm,kferm*sizeof(Complex),cudaMemcpyDefault);
#if (nproc>1)//strided memcpy
   for(unsigned short j=0;j<nc*ngorkov;j++){
      cudaMemcpyAsync(p+j*kvolHalo,xi+j*kvol,kvol*sizeof(Complex),cudaMemcpyDefault,streams[j]);
      cudaMemcpyAsync(p_f+j*kvolHalo,xi_f+j*kvol,kvol*sizeof(Complex_f),cudaMemcpyDefault,streams[j]);
   }
#else
   cudaMemcpyAsync(p,xi,kferm*sizeof(Complex),cudaMemcpyDefault,streams[0]);
   cudaMemcpyAsync(p_f,xi_f,kferm*sizeof(Complex_f),cudaMemcpyDefault,streams[1]);
#endif
   cudaDeviceSynchronise();
#else
   memcpy(r,Phi+na*kferm,kferm*sizeof(Complex));
   //These strided loops for MPI halos are a total pain in the arse...
   for(unsigned short j=0;j<nc*ngorkov;j++){
      memcpy(p+j*kvolHalo,xi+j*kvol,kvol*sizeof(Complex));
      memcpy(p_f+j*kvolHalo,xi_f+j*kvol,kvol*sizeof(Complex_f));
   }
#endif
 
   double betan=1;double beta_max=FLT_MAX; bool do_dp=true;
   for((*itercg)=0; (*itercg)<niterc; (*itercg)++){
      //Don't overwrite on first run. 
      if(do_dp){
#ifdef _DEBUGCG
         if(!rank)
            printf("Going to double precision on iteration %d.  %sbetan %e\talpha %e.  %s",
                  *itercg,endline,betan,creal(alpha),endline);
#elifdef _DEBUG
         if(!rank)
            printf("\nGoing to double precision on iteration %d. %sbetan %e\talpha %e.\n",
                  *itercg,endline,betan,creal(alpha));
#endif
         ComplexConvert(r_f,r,kferm,false,1);
         //TODO: Banking on converting the halo too being faster than multiple launches
         ComplexConvert(p_f,p,kvol,false,nc*ngorkov);
#ifdef USE_GPU
         //Update the residue vector, but not on the first call.
         if(*itercg)
            cuMixed_Sumto((double *)xi,(float *)xi_f,2*kferm,dimGrid,dimBlock);
         //Bring everything into double precision
         //Reset xi_f to zero.
         cudaDeviceSynchronise();
         cudaMemsetAsync(xi_f,0,kferm*sizeof(Complex_f),streams[4]);
#else
         //Update the residue vector, but not on the first call.
         if(*itercg)
#pragma omp parallel for simd aligned(xi,xi_f:AVX)
            for(unsigned int i=0;i<kferm;i++)
               xi[i]+=(Complex)xi_f[i];
         memset(xi_f,0,kferm*sizeof(Complex_f));
#endif
         //No need to synchronise here.  The memcpy in Dslash is blocking
         Dslash(x1,p,ud,iu,id,gamval,gamin,dk,jqq,akappa);
         if(c_sw)
            ByClover(x1,p,clover,sigval,akappa,sigin,false);
         Dslashd(x2,x1,ud,iu,id,gamval,gamin,dk,jqq,akappa);
         if(c_sw)
            ByClover(x2,x1,clover,sigval,akappa,sigin,true);
#ifdef   USE_GPU
         cudaDeviceSynchronise();
#endif

         if(*itercg){
            alpha=0;
#ifdef USE_GPU
#if (nproc>1)//Strided
            for(unsigned short j=0;j<nc*ngorkov;j++){
               alignas(16) Complex alpha_t=0;
               cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)p+j*kvolHalo,1,(cuDoubleComplex *)x2+j*kvol,1,(cuDoubleComplex *)&alpha_t);
               alpha+=alpha_t
            }
#else
            cublasZdotc(cublas_handle,kferm,(cuDoubleComplex *)p,1,(cuDoubleComplex *)x2,1,(cuDoubleComplex *)&alpha);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ngorkov;j++){
               Complex alpha_t=0;
               cblas_zdotc_sub(kvol, p+j*kvolHalo, 1, x2+j*kvol, 1, &alpha_t);
               alpha+=alpha_t;
            }
#else
#pragma omp parallel for simd collapse(2) aligned(p,x2:AVX) reduction(+:alpha)
            for(unsigned short j=0;j<nc*ngorkov;j++)
               for(unsigned int i=0; i<kvol; i++)
                  alpha+=conj(p[i+j*kvolHalo])*x2[i+j*kvol];
#endif
            //For now I'll cast it into a float for the reduction.  Each rank only sends and writes
            //to the real part so this is fine
#if (nproc>1)
            Par_dsum((double *)&alpha);
#endif
            alpha=alphan/creal(alpha);
#ifdef USE_GPU
#if (nproc>1)
            for(unsigned short j=0;j<nc*ngorkov;j++)
               cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&alpha,(cuDoubleComplex *)p+j*kvolHalo,1,(cuDoubleComplex *)xi+j*kvol,1);
#else
            cublasZaxpy(cublas_handle,kferm,(cuDoubleComplex *)&alpha,(cuDoubleComplex *)p,1,(cuDoubleComplex *)xi,1);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ngorkov;j++)
               cblas_zaxpy(kvol, &alpha, p+j*kvolHalo, 1, xi+j*kvol, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(p,xi:AVX)
            for(unsigned short j=0;j<nc*ngorkov;j++)
               for(unsigned int i=0; i<kvol; i++)
                  xi[i+j*kvol]+=alpha*p[i+j*kvolHalo];
#endif
         }
 
#ifdef   USE_GPU
         Complex alpha_m=(Complex)(-alpha);
         cublasZaxpy(cublas_handle, kferm,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)x2,1,(cuDoubleComplex *)r,1);
         double betan_d;
         cublasDznrm2(cublas_handle,kferm,(cuDoubleComplex *)r,1,&betan_d);
         betan=betan_d*betan_d;
#elif defined USE_BLAS
         const Complex alpha_m = (Complex)(-alpha);
         cblas_zaxpy(kferm, &alpha_m, x2, 1, r, 1);
         //Undo the negation for the BLAS routine
         double betan_d = cblas_dznrm2(kferm, r,1);
         //Gotta square it to "undo" the norm
         betan = betan_d*betan_d;
#else
         betan=0;
#pragma omp parallel for simd aligned(r,x2:AVX) reduction(+:betan)
         for(unsigned int i=0; i<kferm; i++){
            r[i]-=alpha*x2[i];
            betan += conj(r[i])*r[i];
         }
         double betan_d=sqrt(betan);
#endif
         //And...  reduce.
#if (nproc>1)
         Par_dsum(&betan);
#endif
         betan_d=sqrt(betan);
         //Update beta_max.  Mandatory for double precision.
         beta_max=betan_d;
#ifdef _DEBUG
         if(! rank) printf("DP Iter(CG)=%i\tbeta_n=%e\talpha=%e\n", *itercg, betan, creal(alpha));
         fflush(stdout);
#endif
         alignas(16) const Complex beta = (*itercg) ?    betan/betad :   0;
         betad=betan; alphan=betan;
#ifdef USE_GPU
         alpha_m=1;
#if (nproc>1)
         for(unsigned short j=0;j<nc*ngorkov;j++){
            cublasZdscal(cublas_handle,kvol,(double *)&beta,(cuDoubleComplex *)p+j*kvolHalo,1);
            cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)r+j*kvol,1,(cuDoubleComplex *)p+j*kvolHalo,1);
         }
#else
         cublasZdscal(cublas_handle,kferm,(double *)&beta,(cuDoubleComplex *)p,1);
         cublasZaxpy(cublas_handle,kferm,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)r,1,(cuDoubleComplex *)p,1);
#endif
#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
         const Complex a = 1.0;
         //There is cblas_? axpby in the MKL and AMD though, set a = 1 and b = \beta.
         //If we get a small enough \beta_n before hitting the iteration cap we break
         for(unsigned short j=0;j<nc*ngorkov;j++)
            cblas_zaxpby(kvol, &a, r+j*kvol, 1, &beta,  p+j*kvolHalo, 1);
#elif defined USE_BLAS
         const Complex a = 1.0;
         for(unsigned short j=0;j<nc*ngorkov;j++){
            cblas_zscal(kvol,&beta,p+j*kvolHalo,1);
            cblas_zaxpy(kvol,&a,r+j*kvol,1,p+j*kvolHalo,1);
         }
#else
#pragma omp parallel for simd collapse(2) aligned(r,p:AVX)
         for(unsigned short j=0;j<nc*ngorkov;j++)
            for(unsigned int i=0; i<kvol; i++)
               p[i+j*kvolHalo]=r[i+j*kvol]+beta*p[i+j*kvolHalo];
#endif
         ComplexConvert(p_f,p,kvol,true,nc*ngorkov);
         ComplexConvert(r_f,r,kferm,true,1);
#ifdef _DEBUGCG
         if(! rank) printf("Double precision.  Iter(CG)=%i\tbeta_n=%e\talpha=%e\n", *itercg, betan, creal(alpha));
         fflush(stdout);
#endif
         if(betan<resid){
            (*itercg)++;
#ifdef _DEBUG
            if(!rank) printf("\nIter(CG)=%i\tResidue:  %e\tTolerance:  %e\n", *itercg, betan, resid);
#endif
            ret_val=0;  break;
         }
         else
            do_dp=false;
      }
      else{
         //x2=(M^\dagger)x1=(M^\dagger)Mp
         Dslash_f(x1_f,p_f,ut,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
         if(c_sw)
            ByClover_f(x1_f,p_f,clover_f,sigval_f,akappa,sigin,false);
         Dslashd_f(x2_f,x1_f,ut,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
         if(c_sw)
            ByClover_f(x2_f,x1_f,clover_f,sigval_f,akappa,sigin,true);
#ifdef USE_GPU
         cudaDeviceSynchronise();
#endif
         //We can't evaluate \alpha on the first niterx because we need to get \beta_n.
         if(*itercg){
            //x*.x
            alphad=0;
#ifdef USE_GPU
#if(nproc>1)//strided
            for(unsigned short j=0;j<nc*ngorkov;j++){
               alignas(8) Complex_f alpha_t=0;
               cublasCdotc(cublas_handle,kvol,(cuComplex *)p_f+j*kvolHalo,1,(cuComplex *)x2_f+j*kvol,1,(cuComplex *)&alpha_t);
               alphad+=alpha_t;
            }
#else
            cublasCdotc(cublas_handle,kferm,(cuComplex *)p_f,1,(cuComplex *)x2_f,1,(cuComplex *)&alphad);
#endif
#elifdef USE_BLAS
            for(unsigned short j=0;j<nc*ngorkov;j++){
               alignas(8) Complex_f alpha_t=0;
               cblas_cdotc_sub(kvol,p_f+j*kvolHalo,1,x2_f+j*kvol,1,&alpha_t);
               alphad+=alpha_t;
            }
#else
#pragma omp parallel for simd collapse(2) aligned(p_f,x2_f:AVX) reduction(+:alphad)
            for(unsigned short j=0;j<nc*ngorkov;j++)
               for(unsigned int i = 0; i<kvol; i++)
                  alphad+=conj(p_f[i+j*kvolHalo])*x2_f[i+j*kvol];
#endif
#if (nproc>1)
            Par_fsum((float *)&alphad);
#endif
            //alpha=(r.r)/p(M^\dagger)Mp
            alpha=alphan/creal(alphad);
            //       Complex_f alpha_f = (Complex_f)alpha;
            //x+\alpha p
#ifdef USE_BLAS
            alignas(8) Complex_f alpha_f=(float)alpha;
#ifdef USE_GPU
#if (nproc>1) //strided
            for(unsigned short j=0;j<nc*ngorkov;j++)
               cublasCaxpy(cublas_handle,kvol,(cuComplex*) &alpha_f,(cuComplex*) p_f+j*kvolHalo,1,(cuComplex*) xi_f+j*kvol,1);
#else
            cublasCaxpy(cublas_handle,kferm,(cuComplex*) &alpha_f,(cuComplex*) p_f,1,(cuComplex*) xi_f,1);
#endif
#else
            for(unsigned short j=0;j<nc*ngorkov;j++)
               cblas_caxpy(kvol, (Complex_f*)&alpha_f,(Complex_f*)p_f+j*kvolHalo, 1, (Complex_f*)xi_f+j*kvol, 1);
#endif
#else
#pragma omp parallel for simd collapse(2) aligned(xi_f,p_f:AVX)
            for(unsigned short j=0;j<nc*ngorkov;j++)
               for(unsigned int i = 0; i<kvol; i++)
                  xi_f[i+j*kvol]+=alpha*p_f[i+j*kvolHalo];
#endif
         }
 
         //r=\alpha(M^\dagger)Mp and \beta_n=r*.r
         alignas(8) float betan_f=0;
#if defined USE_BLAS
         alignas(16) Complex_f alpha_m=(Complex_f)(-alpha);
#ifdef USE_GPU
         cublasCaxpy(cublas_handle,kferm, (cuComplex *)&alpha_m,(cuComplex *) x2_f, 1,(cuComplex *) r_f, 1);
         //cudaDeviceSynchronise();
         //r*.r
         cublasScnrm2(cublas_handle,kferm,(cuComplex *) r_f,1,(float *)&betan_f);
#else
         cblas_caxpy(kferm,(Complex_f*) &alpha_m,(Complex_f*) x2_f, 1,(Complex_f*) r_f, 1);
         //r*.r
         betan_f = cblas_scnrm2(kferm, (Complex_f*)r_f,1);
#endif
         //Gotta square it to "undo" the norm
         betan=betan_f*betan_f;
#else
         //Just like Congradq, this loop could be unrolled but will need a reduction to deal with the betan 
         //addition.
         betan = 0;
         //If we get a small enough \beta_n before hitting the iteration cap we break
#pragma omp parallel for simd aligned(x2_f,r_f:AVX) reduction(+:betan)
         for(unsigned int i = 0; i<kferm;i++){
            r_f[i]-=alpha*x2_f[i];
            betan+=conj(r_f[i])*r_f[i];
         }
         betan_f=sqrt(betan);
#endif
         //This is basically just congradq at the end. Check there for comments
#if (nproc>1)
         Par_dsum(&betan);
#endif
         //Update beta_max if needed. On paper congrad is monotonically decreasing
         betan_f=sqrt(betan);
         beta_max = (betan_f>beta_max) ?betan_f : beta_max;
#ifdef _DEBUG
         if(!rank) printf("Iter (CG) = %i beta_n= %e alpha= %e%s", *itercg, betan, creal(alpha),endline);
#endif
         if(betan_f<beta_max*d_prec){
#ifdef _DEBUG
            if(!rank)
               printf("Residue %e is less than %e times %e. %s",betan_f,d_prec,beta_max,endline);
#endif
            do_dp=true;
         }
         else if(betan<resid){
            //Started counting from zero so add one to make it accurate
            (*itercg)++;
#ifdef _DEBUG
            if(!rank) printf("\nIter (CG) = %i resid = %e toler = %e\n", *itercg, betan, resid);
#endif
            //Final iteration should be in double precision.
            do_dp=true;
         }
         else if(*itercg==niterc-1){
            if(!rank) fprintf(stderr, "Warning %i in %s: Exceeded iteration limit %i beta_n=%e\n",
                  ITERLIM, funcname, niterc, betan);
            ret_val=ITERLIM;  break;
         }
         //Note that beta below is not the global beta and scoping is used to avoid conflict between them
         alignas(8) Complex beta = (*itercg) ? betan/betad : 0;
         betad=betan; alphan=betan;
         //BLAS for p=r+\betap doesn't exist in standard BLAS. This is NOT an axpy case as we're multiplying y by 
         //\beta instead of x.
         //There is cblas_zaxpby in the MKL though, set a = 1 and b = \beta.
#ifdef USE_BLAS
         alignas(8) Complex_f beta_f = (Complex_f)beta;
         alignas(8) Complex_f a = 1.0;
#ifdef USE_GPU
#if (nproc>1) //strided
         for(unsigned short j=0;j<nc*ngorkov;j++){
            cublasCscal(cublas_handle,kvol,(cuComplex *)&beta_f,(cuComplex *)p_f+j*kvolHalo,1);
            cublasCaxpy(cublas_handle,kvol,(cuComplex *)&a,(cuComplex *)r_f+j*kvol,1,(cuComplex *)p_f+j*kvolHalo,1);
         }
#else
         cublasCscal(cublas_handle,kferm,(cuComplex *)&beta_f,(cuComplex *)p_f,1);
         cublasCaxpy(cublas_handle,kferm,(cuComplex *)&a,(cuComplex *)r_f,1,(cuComplex *)p_f,1);
#endif
         cudaDeviceSynchronise();
#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
         for(unsigned short j=0;j<nc*ngorkov;j++)
            cblas_caxpby(kvol, &a, r_f+j*kvol, 1, &beta_f,  p_f+j*kvolHalo, 1);
#else
         for(unsigned short j=0;j<nc*ngorkov;j++){
            cblas_cscal(kvol,&beta_f,p_f+j*kvolHalo,1);
            cblas_caxpy(kvol,&a,r_f+j*kvol,1,p_f+j*kvolHalo,1);
         }
#endif
#else
#pragma omp parallel for simd aligned(r_f,p_f:AVX)
         for(unsigned short j=0;j<nc*ngorkov;j++)
            for(unsigned int i=0; i<kvol; i++)
               p_f[i+j*kvolHalo]=r_f[i+j*kvol]+beta*p_f[i+j*kvolHalo];
#endif
      }
   }
#ifdef USE_GPU
   cudaDeviceSynchronise();
#endif
   P_free_f(&p_f,&r_f,&x1_f,&x2_f,&xi_f);
   P_free(&p,&r,&x1,&x2,clover);
   return ret_val;
}

References ByClover(), ByClover_f(), Complex, Complex_f, ComplexConvert(), conj(), creal, cublas_handle, cudaDeviceSynchronise, cuMixed_Sumto(), dimBlock, dimGrid, Dslash(), Dslash_f(), Dslashd(), Dslashd_f(), ITERLIM, kferm, kvol, kvolHalo, nc, ngorkov, niterc, P_allocate(), P_allocate_f(), P_free(), P_free_f(), Par_dsum(), Par_fsum(), rank, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

Congradq()

int Congradq	(	int	na,
		double	res,
		Complex *	X1,
		Complex *	r,
		Complex *	ud[2],
		Complex_f *	ut[2],
		Complex_f *	clover_f[nc],
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		double *	dk[2],
		float *	dk_f[2],
		Complex_f	jqq,
		float	akappa,
		float	c_sw,
		int *	itercg )

Matrix Inversion via Conjugate Gradient (up/down flavour partitioning). Solves \((M^\dagger)Mx=\Phi\) Implements up/down partitioning The matrix multiplication step is done at mixed precision, while the update is done at double.

Parameters

[in]	na	Flavour index
[in]	res	Limit for conjugate gradient
[in,out]	X1	Pseudofermion field \(\Phi\) initially, returned as \((M^\dagger M)^{-1} \Phi\)
[in]	r	Partition of \(\Phi\) being used. Gets recycled as the residual vector
[in]	ud,ut	Double/float Trial colour fields
[in]	iu,id	Upper/lower halo indices
[in]	gamval,gamval_f	Double/float gamma matrices rescaled by kappa
[in]	gamin	What element of the spinor is multiplied by row idirac each gamma matrix?
[in]	clover_f	Array of clover fields
[in]	sigval,sigval_f	Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	dk,dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]	jqq	Diquark source
[in]	akappa	Hopping Parameter
[in]	c_sw	Clover coefficient. If non-zero calculate the clover contribution
[in]	itercg	Counts the iterations of the conjugate gradient

Returns

0 on success, integer error code otherwise

Postcondition

Contents of X1 and r overwritten

TODO: Set this as an argument How much does the residue have to shrink by before we do a double precision update

The \(\kappa^2\) factor is needed to normalise the fields correctly \(j_{qq}\) is the diquark condensate and is global scope.

The alpha and beta terms should be double, but that causes issues with BLAS pointers. Instead we declare them complex and work with the real part (especially for \(\alpha_d\)) Give initial values Will be overwritten if niterx>0

\(x2 = (M^\dagger M)p \)

\(\alpha_d= p* (M^\dagger M+c_\text{SW} \sum\limits_{\mu\ne\nu}\frac{1}{2}\sigma_{\mu\nu}F_{\mu\nu}+J^2)p\)

\(\alpha=\frac{\alpha_n}{\alpha_d}=\frac{r\cdot r}{p(M^\dagger M+c_\text{SW} \sum\limits_{\mu\ne\nu}\frac{1}{2}\sigma_{\mu\nu}F_{\mu\nu}+J^2)p}\)

\(x+\alpha p\)

\(x2 = (M^\dagger M)p \)

\(x2 = (M^\dagger M+J^2)p\)

First iteration \( \vec{r}_0=\vec{p}_0\)
\(\alpha_d= p* (M^\dagger M+c_\text{SW} \sum\limits_{\mu\ne\nu}\frac{1}{2}\sigma_{\mu\nu}F_{\mu\nu}+J^2)p\)

\(\alpha=\frac{\alpha_n}{\alpha_d}=\frac{r\cdot r}{p(M^\dagger M+c_\text{SW} \sum\limits_{\mu\ne\nu}\frac{1}{2}\sigma_{\mu\nu}F_{\mu\nu}+J^2)p}\)

\(x+\alpha p\)

\(r_{n+1} = r_n-\alpha(M^\dagger M)p_n\) and \(\beta_n=r^\dagger r\)

Definition at line 278 of file congrad.c.

                                                        {
   const char funcname[] = "Congradq";
   int ret_val=0;
   const double resid = res*res;
#ifdef _DEBUGCG
#warning "CG Debugging"
   char *endline = "\n";
#else
   char *endline = "\r";
#endif

   const float d_prec=1.0f/128.0f;
   const Complex_f fac_f = conj(jqq)*jqq*akappa*akappa;
   //These were evaluated only in the first loop of niterx so we'll just do it outside of the loop.
   //n suffix is numerator, d is denominator
   alignas(16) double alphan=1;
   //Alignment needed for cuBLAS
   alignas(16) double betad = 1.0; alignas(8) Complex_f alphad=0; alignas(16) Complex alpha = 1;
   //Because we're dealing with flattened arrays here we can call cblas safely without the halo
 
   Complex_f *p_f, *x1_f, *x2_f, *r_f, *X1_f;
   Complex   *p, *x1, *x2, *clover[2];
   Q_allocate_f(&p_f,&x1_f,&x2_f,&r_f,&X1_f);
   Q_allocate(&p,&x1,&x2,clover);
 
   //Instead of copying element-wise in a loop, use memcpy.
   //Get X1 in single precision
   //Since X1 has a halo and X1_f does not we have to do this manually
#if (nproc>1)
   for(unsigned short j=0;j<nc*ndirac;j++)
      ComplexConvert(X1_f+j*kvol,X1+j*kvolHalo,kvol,true,1);
#else
   ComplexConvert(X1_f,X1,kferm2,true,1);
#endif
   ComplexConvert(r_f,r,kferm2,true,1);
   //And clover in double
   if(c_sw){
      ComplexConvert(clover_f[0],clover[0],6*kvol,false,1);
      ComplexConvert(clover_f[1],clover[1],6*kvol,false,1);
   }
#ifdef USE_GPU
   //Ensure conversion is done
   cudaDeviceSynchronise();
   //Needs to be strided
#if (nproc>1)
   for(unsigned int j=0;j<nc*ndirac;j++)
      cudaMemcpyAsync(p_f+j*kvolHalo, X1_f+j*kvol, kvol*sizeof(Complex_f),cudaMemcpyDefault,streams[j]);
   cudaDeviceSynchronise();
#else
   cudaMemcpy(p_f, X1_f, kferm2*sizeof(Complex_f),cudaMemcpyDefault);
#endif
#else
#pragma omp parallel for simd aligned(X1_f,r_f,X1,r:AVX)
   for(unsigned int j=0;j<nc*ndirac;j++)
      memcpy(p_f+j*kvolHalo, X1_f+j*kvol, kvol*sizeof(Complex_f));
#endif
 
   alignas(16) double betan=1;double beta_max=FLT_MAX; bool do_dp=true;
   for(*itercg=0; *itercg<niterc; (*itercg)++){
      if(do_dp){
#ifdef _DEBUGCG
         if(!rank)
            printf("Going to double precision on iteration %d. betan %e\talpha %e.\n",
                  *itercg,betan,creal(alpha));
#elifdef _DEBUG
         if(!rank)
            printf("\nGoing to double precision on iteration %d. betan %e\talpha %e.\n",
                  *itercg,betan,creal(alpha));
#endif
         ComplexConvert(r_f,r,kferm2,false,1);
         ComplexConvert(p_f,p,kvol,false,nc*ndirac);
#ifdef USE_GPU
         //Update the residue vector, but not on the first call.
         //TODO: Check for multi-gpu. I fear this will get messy
         if(*itercg)
            cuMixed_Sumto((double *)X1,(float *)X1_f,2*kferm2,dimGrid,dimBlock);
         //Bring everything into double precision
         //Reset X1_f to zero.
         cudaMemset(X1_f,0,kferm2*sizeof(Complex_f));
#else
         //Update the residue vector, but not on the first call.
         if(*itercg)
#pragma omp parallel for simd collapse(2) aligned(X1,X1_f:AVX)
            for(unsigned short j=0;j<nc*ndirac;j++)
               for(unsigned int i=0;i<kvol;i++){
                  X1[i+j*kvolHalo]+=(Complex)X1_f[i+j*kvol];
               }
         memset(X1_f,0,kferm2*sizeof(Complex_f));
#endif
         //No need to synchronise here. The memcpy in Hdslash is blocking
         Hdslash(x1,p,ud,iu,id,gamval,gamin,dk,akappa);
         //Clover contribution
         if(c_sw)
            HbyClover(x1,p,clover,sigval,akappa,sigin,false);
         Hdslashd(x2,x1,ud,iu,id,gamval,gamin,dk,akappa);
         //Clover contribution
         if(c_sw)
            HbyClover(x2,x1,clover,sigval,akappa,sigin,true);
#ifdef   USE_GPU
         cudaDeviceSynchronise();
#endif
         if(fac_f!=0){
            alignas(16) const double fac=(double)fac_f;
#ifdef   USE_GPU
            //Multiple ranks means we need striding
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++) 
               cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&fac,(cuDoubleComplex *)p+j*kvolHalo,1,(cuDoubleComplex *)x2+j*kvol,1);
#else
            //Single GPU, no halos, so make it one call to speed things up
            cublasZaxpy(cublas_handle,kferm2,(cuDoubleComplex *)&fac,(cuDoubleComplex *)p,1,(cuDoubleComplex *)x2,1);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ndirac;j++) 
               cblas_zaxpy(kvol, &fac, p+j*kvolHalo, 1, x2+j*kvol, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(p,x2:AVX)
            for(unsigned short j=0;j<nc*ndirac;j++) 
               for(unsigned int i=0; i<kvol; i++)
                  x2[i+j*kvol]+=fac*p[i+j*kvolHalo];
#endif
         }

         if(*itercg){
            alpha=0;
#ifdef USE_GPU
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++){
               Complex alpha_t=0;
               cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)p+j*kvolHalo,1,(cuDoubleComplex *)x2+j*kvol,1,(cuDoubleComplex *)&alpha_t);
               alpha+=alpha_t;
            }
#else
            cublasZdotc(cublas_handle,kferm2,(cuDoubleComplex *)p,1,(cuDoubleComplex *)x2,1,(cuDoubleComplex *)&alpha);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ndirac;j++){
               Complex alpha_t=0;
               cblas_zdotc_sub(kvol, p+j*kvolHalo, 1, x2+j*kvol, 1, &alpha_t);
               alpha+=creal(alpha_t);
            }
#else
#pragma omp parallel for simd collapse(2) aligned(p,x2:AVX) reduction(+:alpha)
            for(unsigned short j=0;j<nc*ndirac;j++) 
               for(unsigned int i=0; i<kvol; i++)
                  alpha+=conj(p[i+j*kvolHalo])*x2[i+j*kvol];
#endif
            //For now I'll cast it into a float for the reduction. Each rank only sends and writes
            //to the real part so this is fine
#if (nproc>1)
            Par_dsum((double *)&alpha);
#endif
            alpha=alphan/creal(alpha);
#ifdef USE_GPU
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++)
               cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&alpha,(cuDoubleComplex *)p+j*kvolHalo,1,(cuDoubleComplex *)X1+j*kvolHalo,1);
#else
            cublasZaxpy(cublas_handle,kferm2,(cuDoubleComplex *)&alpha,(cuDoubleComplex *)p,1,(cuDoubleComplex *)X1,1);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ndirac;j++)
               cblas_zaxpy(kvol, &alpha, p+j*kvolHalo, 1, X1+j*kvolHalo, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(p,X1:AVX)
            for(unsigned short j=0;j<nc*ndirac;j++)
               for(unsigned int i=0; i<kvol; i++)
                  X1[i+j*kvolHalo]+=alpha*p[i+j*kvolHalo];
#endif
         }
 
#ifdef   USE_GPU
         alignas(16) Complex alpha_m=(Complex)(-alpha);
         cublasZaxpy(cublas_handle, kferm2,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)x2,1,(cuDoubleComplex *)r,1);
         alignas(16) double betan_d;
         cublasDznrm2(cublas_handle,kferm2,(cuDoubleComplex *)r,1,&betan_d);
         betan=betan_d*betan_d;
#elif defined USE_BLAS
         const Complex alpha_m = (Complex)(-alpha);
         cblas_zaxpy(kferm2, &alpha_m, x2, 1, r, 1);
         //Undo the negation for the BLAS routine
         double betan_d = cblas_dznrm2(kferm2, r,1);
         //Gotta square it to "undo" the norm
         betan = betan_d*betan_d;
#else
         betan=0;
#pragma omp parallel for simd aligned(r_f,x2_f:AVX) reduction(+:betan) 
         for(unsigned int i=0; i<kferm2; i++){
            r[i]-=alpha*x2[i];
            betan += conj(r[i])*r[i];
         }
         alignas(16) double betan_d=sqrt(betan);
#endif
         //And... reduce.
#if (nproc>1)
         Par_dsum(&betan);
#endif
         betan_d=sqrt(betan);
         //Update beta_max. Mandatory for double precision.
         beta_max=betan_d;
#ifdef _DEBUG
         if(!rank) printf("DP Iter(CG)=%i\tbeta_n=%e\talpha=%e\n", *itercg, betan, creal(alpha));
         fflush(stdout);
#endif
         alignas(16) const Complex beta = (*itercg) ?  betan/betad : 0;
         betad=betan; alphan=betan;
#ifdef USE_GPU
         alpha_m=1;
         //Strided multi-gpu
#if (nproc>1)
         for(unsigned short j=0;j<nc*ndirac;j++){
            cublasZdscal(cublas_handle,kvol,(double *)&beta,(cuDoubleComplex *)p+j*kvolHalo,1);
            cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)r+j*kvol,1,(cuDoubleComplex *)p+j*kvolHalo,1);
         }
         //And single GPU
#else
         cublasZdscal(cublas_handle,kferm2,(double *)&beta,(cuDoubleComplex *)p,1);
         cublasZaxpy(cublas_handle,kferm2,(cuDoubleComplex *)&alpha_m,(cuDoubleComplex *)r,1,(cuDoubleComplex *)p,1);
#endif
#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
         const Complex a = 1.0;
         //There is cblas_?axpby in the MKL and AMD though, set a = 1 and b = \beta.
         //If we get a small enough \beta_n before hitting the iteration cap we break
         for(unsigned short j=0;j<nc*ndirac;j++)
            cblas_zaxpby(kvol, &a, r+j*kvol, 1, &beta,  p+j*kvolHalo, 1);
#elifdef USE_BLAS
         const Complex a = 1.0;
         for(unsigned short j=0;j<nc*ndirac;j++){
            cblas_zscal(kvol,&beta,p+j*kvolHalo,1);
            cblas_zaxpy(kvol,&a,r+j*kvol,1,p+j*kvolHalo,1);
         }
#else 
#pragma omp parallel for simd collapse(2) aligned(r,p:AVX)
         for(unsigned short j=0;j<nc*ndirac;j++)
            for(unsigned int i=0; i<kvol; i++)
               p[i+j*kvolHalo]=r[i+j*kvol]+beta*p[i+j*kvolHalo];
#endif
         ComplexConvert(p_f,p,kvol,true,nc*ndirac);
         ComplexConvert(r_f,r,kferm2,true,1);
#ifdef _DEBUGCG
         if(!rank) printf("Double precision. Iter(CG)=%i\tbeta_n=%e\talpha=%e\n", *itercg, betan, creal(alpha));
         fflush(stdout);
#endif
         if(betan<resid){ 
            (*itercg)++;
#ifdef _DEBUG
            if(!rank) printf("\nIter(CG)=%i\tResidue: %e\tTolerance: %e\n", *itercg, betan, resid);
#endif
            ret_val=0;  break;
         }
         else
            do_dp=false;
      }
      else{
         //No need to synchronise here. The memcpy in Hdslash is blocking
         Hdslash_f(x1_f,p_f,ut,iu,id,gamval_f,gamin,dk_f,akappa);
         //Clover contribution
         if(c_sw)
            HbyClover_f(x1_f,p_f,clover_f,sigval_f,akappa,sigin,false);
         Hdslashd_f(x2_f,x1_f,ut,iu,id,gamval_f,gamin,dk_f,akappa);
         //Clover contribution
         if(c_sw)
            HbyClover_f(x2_f,x1_f,clover_f,sigval_f,akappa,sigin,true);
#ifdef   USE_GPU
         cudaDeviceSynchronise();
#endif
         //No point adding zero a couple of hundred times if the diquark source is zero
         if(fac_f!=0){
#ifdef   USE_GPU
            //Strided multi-gpu
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++) 
               cublasCaxpy(cublas_handle,kvol,(cuComplex *)&fac_f,(cuComplex *)p_f+j*kvolHalo,1,(cuComplex *)x2_f+j*kvol,1);
            //Single GPU
#else
            cublasCaxpy(cublas_handle,kferm2,(cuComplex *)&fac_f,(cuComplex *)p_f,1,(cuComplex *)x2_f,1);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ndirac;j++) 
               cblas_caxpy(kvol, &fac_f, p_f+j*kvolHalo, 1, x2_f+j*kvol, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(p_f,x2_f:AVX)
            for(unsigned short j=0;j<nc*ndirac;j++) 
               for(unsigned int i=0; i<kvol; i++)
                  x2_f[i+j*kvol]+=fac_f*p_f[i+j*kvolHalo];
#endif
         }
         if(*itercg){
            alphad=0;
#ifdef USE_GPU
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++){
               Complex alpha_t=0;
               cublasCdotc(cublas_handle,kvol,(cuComplex *)p_f+j*kvolHalo,1,(cuComplex *)x2_f+j*kvol,1,(cuComplex *)&alpha_t);
               alphad+=alpha_t;
            }
#else
            cublasCdotc(cublas_handle,kferm2,(cuComplex *)p_f,1,(cuComplex *)x2_f,1,(cuComplex *)&alphad);
#endif
#elif defined USE_BLAS
            for(unsigned short j=0;j<nc*ndirac;j++){
               Complex_f alpha_t=0;
               cblas_cdotc_sub(kvol, p_f+j*kvolHalo, 1, x2_f+j*kvol, 1, &alpha_t);
               alphad+=creal(alpha_t);
            }
#else
#pragma omp parallel for simd aligned(p_f,x2_f:AVX) reduction(alphad:+)
            for(unsigned short j=0;j<nc*ndirac;j++)
               for(unsigned int i=0; i<kvol; i++)
                  alphad+=conj(p_f[i+j*kvolHalo])*x2_f[i+j*kvol];
#endif
            //For now I'll cast it into a float for the reduction. Each rank only sends and writes
            //to the real part so this is fine
#if (nproc>1)
            Par_fsum((float *)&alphad);
#endif
            alpha=alphan/creal(alphad);
#ifdef USE_GPU
            alignas(8) Complex_f alpha_f = (Complex_f)alpha;
#if (nproc>1)
            for(unsigned short j=0;j<nc*ndirac;j++)
               cublasCaxpy(cublas_handle,kvol,(cuComplex *)&alpha_f,(cuComplex *)p_f+j*kvolHalo,1,(cuComplex *)X1_f+j*kvol,1);
#else
            cublasCaxpy(cublas_handle,kferm2,(cuComplex *)&alpha_f,(cuComplex *)p_f,1,(cuComplex *)X1_f,1);
#endif
#elif defined USE_BLAS
            Complex_f alpha_f = (Complex_f)alpha;
            for(unsigned short j=0;j<nc*ndirac;j++)
               cblas_caxpy(kvol, &alpha_f, p_f+j*kvolHalo, 1, X1_f+j*kvol, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(X1_f,p_f:AVX)
            for(unsigned short j=0;j<nc*ndirac;j++)
               for(unsigned int i=0; i<kvol; i++)
                  X1_f[i+j*kvol]+=alpha*p_f[i+j*kvolHalo];
#endif
         }        
         // And no Halos here so nice and easy
#ifdef   USE_GPU
         alignas(8) __managed__ Complex_f alpha_m=(Complex_f)(-alpha);
         cublasCaxpy(cublas_handle, kferm2,(cuComplex *)&alpha_m,(cuComplex *)x2_f,1,(cuComplex *)r_f,1);
         alignas(8) float betan_f;
         cublasScnrm2(cublas_handle,kferm2,(cuComplex *)r_f,1,&betan_f);
         betan = betan_f*betan_f;
#elif defined USE_BLAS
         const Complex_f alpha_m = (Complex_f)(-alpha);
         cblas_caxpy(kferm2, &alpha_m, x2_f, 1, r_f, 1);
         //Undo the negation for the BLAS routine
         float betan_f = cblas_scnrm2(kferm2, r_f,1);
         //Gotta square it to "undo" the norm
         betan = betan_f*betan_f;
#else
         betan=0;
#pragma omp parallel for simd aligned(r_f,x2_f:AVX) reduction(+:betan) 
         for(unsigned int i=0; i<kferm2; i++){
            r_f[i]-=alpha*x2_f[i];
            betan += conj(r_f[i])*r_f[i];
         }
         alignas(8) float betan_f=sqrt(betan);
#endif
         //And... reduce.
#if (nproc>1)
         Par_dsum(&betan);
#endif
         betan_f=sqrt(betan);
         //Update beta_max if needed. On paper congrad is monotonically decreasing
         beta_max = (betan_f>beta_max) ?betan_f : beta_max;
 
#ifdef _DEBUG
         if(!rank) printf("Iter(CG)=%i\tbeta_n=%e\talpha=%e%s", *itercg, betan, creal(alpha),endline);
         fflush(stdout);
#endif
         if(betan_f<beta_max*d_prec){ 
#ifdef _DEBUG
            if(!rank)
               printf("\nResidue %e is less than %e times %e.%s",betan_f,d_prec,beta_max,endline);
#endif
            do_dp=true;
         }
         else if(betan<resid){
#ifdef _DEBUG
            if(!rank)
               printf("\nBetan %e is less than target residue %e.\n",betan,resid);
#endif
            do_dp=true;
            //break;
         }
         else if(*itercg==niterc-1){
            if(!rank) fprintf(stderr, "\nWarning %i in %s: Exceeded iteration limit %i beta_n=%e\n",\
                  ITERLIM, funcname, *itercg, betan);
            ret_val=ITERLIM;  break;
         }
         //Here we evaluate beta=(r_{k+1}.r_{k+1})/(r_k.r_k) and then shuffle our indices down the line.
         //On the first iteration we define beta to be zero.
         //Note that beta below is not the global beta and scoping is used to avoid conflict between them
         Complex beta = (*itercg) ?  betan/betad : 0;
         betad=betan; alphan=betan;
         //BLAS for p=r+\betap doesn't exist in standard BLAS. This is NOT an axpy case as we're multiplying y by
         //\beta instead of x.
#ifdef USE_GPU
         alignas(8) Complex_f beta_f=(Complex_f)beta;
         alpha_m = 1.0;
#if (nproc>1)
         for(unsigned short j=0;j<nc*ndirac;j++){
            cublasCscal(cublas_handle,kvol,(cuComplex *)&beta_f,(cuComplex *)p_f+j*kvolHalo,1);
            cublasCaxpy(cublas_handle,kvol,(cuComplex *)&alpha_m,(cuComplex *)r_f+j*kvolHalo,1,(cuComplex *)p_f+j*kvolHalo,1);
         }
#else
         cublasCscal(cublas_handle,kferm2,(cuComplex *)&beta_f,(cuComplex *)p_f,1);
         cublasCaxpy(cublas_handle,kferm2,(cuComplex *)&alpha_m,(cuComplex *)r_f,1,(cuComplex *)p_f,1);
#endif
#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
         Complex_f a = 1.0;
         Complex_f beta_f=(Complex_f)beta;
         //There is cblas_?axpby in the MKL and AMD though, set a = 1 and b = \beta.
         //If we get a small enough \beta_n before hitting the iteration cap we break
         for(unsigned short j=0;j<nc*ndirac;j++)
            cblas_caxpby(kvol, &a, r_f+j*kvol, 1, &beta_f,  p_f+j*kvolHalo, 1);
#elif defined USE_BLAS
         Complex_f beta_f=(Complex_f)beta;
         Complex_f a = 1.0;
         for(unsigned short j=0;j<nc*ndirac;j++){
            cblas_cscal(kvol,&beta_f,p_f+j*kvolHalo,1);
            cblas_caxpy(kvol,&a,r_f+j*kvol,1,p_f+j*kvolHalo,1);
         }
#else 
#pragma omp parallel for simd collapse(2) aligned(p_f,r_f:AVX)
         for(unsigned short j=0;j<nc*ndirac;j++)
            for(unsigned int i=0; i<kvol; i++)
               p_f[i+j*kvolHalo]=r_f[i+j*kvol]+beta*p_f[i+j*kvolHalo];
#endif
      }
   }
   Q_free_f(&p_f,&x1_f,&x2_f,&r_f,&X1_f);
   Q_free(&p,&x1,&x2,clover);
   return ret_val;
}

References Complex, Complex_f, ComplexConvert(), conj(), creal, cublas_handle, cudaDeviceSynchronise, cuMixed_Sumto(), dimBlock, dimGrid, HbyClover(), HbyClover_f(), Hdslash(), Hdslash_f(), Hdslashd(), Hdslashd_f(), ITERLIM, kferm2, kvol, kvolHalo, nc, ndirac, niterc, Par_dsum(), Par_fsum(), Q_allocate(), Q_allocate_f(), Q_free(), Q_free_f(), rank, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuForce()

void cuForce	(	double *	dSdpi,
		Complex_f *	ut[2],
		Complex_f *	X1,
		Complex_f *	X2,
		Complex_f	gamval[20],
		float *	dk[2],
		unsigned int *	iu,
		const unsigned short	gamin[16],
		float	akappa,
		dim3	dimGrid,
		dim3	dimBlock )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	ut	Gauge fields
[in]	X1	Inverted field
[in]	X2	\(MX_1\)
[in]	gamval	Double/float precision gamma matrices rescaled by \(\kappa\)
[in]	gamin	Gamma indices
[in]	iu	Lattice indices
[in]	dk	\(e^{-\mu}\) and \(e^\mu\)
[in]	akappa	Hopping parameter
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition

Fermion force added onto dSdpi

Definition at line 373 of file cuforce.cu.

                                                {
   const char *funcname = "Force";
   //X1=(M†M)^{1} Phi
   // Transpose_z(X1,ndirac*nc,kvol); Transpose_z(X2,ndirac*nc,kvol);
   cudaDeviceSynchronise();
#pragma unroll
   for(unsigned short mu=0;mu<3;mu++){
      Kernels::cuForce_s<<<dimGrid,dimBlock,0,streams[mu]>>>(dSdpi,ut[0],ut[1],X1,X2,gamval,iu,gamin,akappa,mu);
   }
   //Set stream for time direction
   unsigned short mu=3;
   Kernels::cuForce_t<<<dimGrid,dimBlock,0,streams[mu]>>>(dSdpi,ut[0],ut[1],X1,X2,gamval,dk[0],dk[1],iu,gamin,akappa);
   cudaDeviceSynchronise();
}

References Complex_f, cudaDeviceSynchronise, Kernels::cuForce_s(), Kernels::cuForce_t(), dimBlock, and dimGrid.

Here is the call graph for this function:

Here is the caller graph for this function:

cuForce_s()

__global__ void Kernels::cuForce_s	(	double *	dSdpi,
		Complex_f *	u11t,
		Complex_f *	u12t,
		Complex_f *	X1,
		Complex_f *	X2,
		Complex_f	gamval[20],
		unsigned int *	iu,
		const unsigned short	gamin[16],
		float	akappa,
		const unsigned short	mu )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	u11t,u12t	Float precision colour fields
[in]	X1	Inverted field
[in]	X2	\(MX_1\)
[in]	gamval	Gamma matrices rescaled by \(\kappa\)
[in]	iu	Lattice indices
[in]	gamin	Gamma indices
[in]	akappa	Hopping parameter
[in]	mu	Force direction

Postcondition

Force added to dSdpi

Definition at line 146 of file cuforce.cu.

                                                                                                {
      const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
      const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
      const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
      const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
      const unsigned int gthreadId= blockId * bsize+bthreadId;
 
      for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
         const unsigned int ind=i+kvolHalo*mu;
         const Complex_f u11s=u11t[ind]; const Complex_f u12s=u12t[ind];
         const unsigned int uid = iu[i+kvol*mu];
         //Similarly to Hdslash we always see idirac*nc so we do that here too.
         for(unsigned short idirac=0;idirac<nc*ndirac;idirac+=nc){
            Complex_f X1s[nc];    Complex_f X1su[nc];
            Complex_f X2s[nc];    Complex_f X2su[nc];
 
            X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
            X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
            X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
            X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
 
            //       Need to be double to avoid accumulation errors
            double dSdpis[3];
            //Careful!! cant use ind here as dSdpi has no halo!
            dSdpis[0]=dSdpi[i+kvol*mu];
            //Multiplying by i and taking the real component is the same as taking the negative imaginary component
            //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
            //generators by hand
            dSdpis[0]+=-akappa*(
                  conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                  +conj(X1s[1])*(u11s*X2su[0]+u12s*X2su[1])
                  +conj(X1su[0])*(u12s*X2s[0]-conj(u11s)*X2s[1])
                  +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1])).imag();
 
            dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
            dSdpis[1]+=akappa*(
                  (conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                   +conj(X1s[1])*(-u11s*X2su[0]-u12s*X2su[1])
                   +conj(X1su[0])*(-u12s*X2s[0]-conj(u11s)*X2s[1])
                   +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1]))).real();
 
            dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
            dSdpis[2]+=-akappa*(
                  conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
                  +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
                  +conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
                  +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1])).imag();
 
            const unsigned short gindex=mu*ndirac+(idirac>>1);
            const Complex_f gamval_c=gamval[gindex];
            //Rescaling gind by nc
            const unsigned short gind = gamin[gindex]<<1;   
            X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
            X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
 
            //If you are asked to rederive the force from Montvay and Munster you'll notice that it should be kappa*gamma
            //but below is only gamma. We rescaled gamma by kappa already when we defined it so that's where it has gone
            dSdpis[0]+=-(gamval_c*
                  (conj(X1s[0])* (-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                   +conj(X1s[1])* (u11s *X2su[0]+u12s *X2su[1])
                   +conj(X1su[0])* (-u12s *X2s[0] +conj(u11s)*X2s[1])
                   +conj(X1su[1])*(u11s *X2s[0] +conj(u12s)*X2s[1]))).imag();
            dSdpi[i+kvol*mu]=dSdpis[0];
 
            dSdpis[1]+=(gamval_c*
                  (conj(X1s[0])* (-conj(u12s)*X2su[0] +conj(u11s)*X2su[1])
                   +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1])
                   +conj(X1su[0])* (u12s *X2s[0]+conj(u11s)*X2s[1])
                   +conj(X1su[1])* (-u11s *X2s[0]+conj(u12s)*X2s[1]))).real();
            dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
 
            dSdpis[2]+=-(gamval_c*
                  (conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
                   +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
                   +conj(X1su[0])*(conj(u11s)*X2s[0]+u12s *X2s[1])
                   +conj(X1su[1])*(conj(u12s)*X2s[0]-u11s *X2s[1]))).imag();
            dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
         }
      }
   }

References Complex_f, conj(), kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

cuForce_t()

__global__ void Kernels::cuForce_t	(	double *	dSdpi,
		Complex_f *	u11t,
		Complex_f *	u12t,
		Complex_f *	X1,
		Complex_f *	X2,
		Complex_f	gamval[20],
		float *	dk4m,
		float *	dk4p,
		unsigned int *	iu,
		const unsigned short	gamin[16],
		float	akappa )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	u11t,u12t	Float precision colour fields
[in]	X1	Inverted field
[in]	X2	\(MX_1\)
[in]	gamval	Gamma matrices rescaled by \(\kappa\)
[in]	dk4m,dk4p	\(e^{-\mu}\) and \(e^\mu\)
[in]	iu	Lattice indices
[in]	gamin	Gamma indices
[in]	akappa	Hopping parameter

Postcondition

Force added to dSdpi

Definition at line 243 of file cuforce.cu.

                                                                                                 {
      const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
      const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
      const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
      const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
      const unsigned int gthreadId= blockId * bsize+bthreadId;
 
      const unsigned short mu=3;
      for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
         const unsigned int ind=i+kvolHalo*mu;
         const Complex_f u11s=u11t[ind];  const Complex_f u12s=u12t[ind];
         const float dk4ms=dk4m[i]; const float dk4ps=dk4p[i];
         //Up indices
         const unsigned int uid = iu[i+kvol*mu];
         //Similarly to Hdslash we always see idirac*nc so we do that here too.
         for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
            Complex_f X1s[nc];    Complex_f X1su[nc];
            Complex_f X2s[nc];    Complex_f X2su[nc];
 
            X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
            X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
            X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
            X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
 
            //       Need to be double to avoid accumulation errors
            double dSdpis[3];
            dSdpis[0]=dSdpi[i+kvol*mu];
            //Multiplying by i and taking the real component is the same as taking the negative imaginary component
            //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
            //generators by hand
            dSdpis[0]+=-(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                     +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
                  +dk4ps*(conj(X1su[0])*(+u12s*X2s[0]-conj(u11s)*X2s[1])
                     +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1]))).imag();
 
            dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
            dSdpis[1]+=(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                     +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1]))
                  +dk4ps*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
                     +conj(X1su[1])*( u11s *X2s[0]-conj(u12s)*X2s[1]))).real();
 
            dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
            dSdpis[2]+=-(dk4ms* (conj(X1s[0])* (u11s *X2su[0]+u12s *X2su[1])
                     +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
                  +dk4ps*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
                     +conj(X1su[1])* (-conj(u12s)*X2s[0]+u11s *X2s[1]))).imag();
 
            const unsigned short gindex=mu*ndirac+(idirac>>1);
            //Rescaling gind by nc
            const unsigned short gind = gamin[gindex]<<1;   
            X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
            X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
 
            dSdpis[0]+=-(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                     +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
                  -dk4ps*(conj(X1su[0])* (u12s *X2s[0]-conj(u11s)*X2s[1])
                     +conj(X1su[1])*(-u11s *X2s[0]-conj(u12s)*X2s[1]))).imag();
            dSdpi[i+kvol*mu]=dSdpis[0];
 
            dSdpis[1]+=(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                     +conj(X1s[1])*(-u11s*X2su[0]-u12s *X2su[1]))
                  -dk4ps*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
                     +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1]))).real();
            dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
 
            dSdpis[2]+=-(dk4ms*(conj(X1s[0])*(u11s*X2su[0] +u12s *X2su[1])
                     +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
                  -dk4ps*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
                     +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1]))).imag();
            dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
         }
      }
   }

References Complex_f, conj(), kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

cuGauge_force()

void cuGauge_force	(	Complex_f *	ut[2],
		double *	dSdpi,
		float	beta,
		unsigned int *	iu,
		unsigned int *	id,
		dim3	dimGrid,
		dim3	dimBlock )

Calculate the gauge contribution to the force.

Parameters

[in]	ut	Gauge fields
[out]	dSdpi	Force
[in]	beta	Inverse gauge coupling
[in]	iu,id	Upper/lower indices
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition

Contents of dSdpi overwritten

Definition at line 319 of file cuforce.cu.

                                                                                                                           {
   const char funcname[] = "Gauge_force";
   int device=-1;
   cudaGetDevice(&device);
   Complex_f *Sigma[ndim][2], *ush[ndim][2];
   for(unsigned short i=0;i<ndim;i++){
#ifdef _DEBUG
      cudaMallocManaged((void **)&Sigma[i][0],kvol*sizeof(Complex_f),cudaMemAttachGlobal);
      cudaMallocManaged((void **)&Sigma[i][1],kvol*sizeof(Complex_f),cudaMemAttachGlobal);
      cudaMallocManaged((void **)&ush[i][0],kvolHalo*sizeof(Complex_f),cudaMemAttachGlobal);
      cudaMallocManaged((void **)&ush[i][1],kvolHalo*sizeof(Complex_f),cudaMemAttachGlobal);
#else
      cudaMallocAsync((void **)&Sigma[i][0],kvol*sizeof(Complex_f),streams[i]);
      cudaMallocAsync((void **)&Sigma[i][1],kvol*sizeof(Complex_f),streams[i]);
      cudaMallocAsync((void **)&ush[i][0],kvolHalo*sizeof(Complex_f),streams[i]);
      cudaMallocAsync((void **)&ush[i][1],kvolHalo*sizeof(Complex_f),streams[i]);
#endif
   }
   for(unsigned short mu=0; mu<ndim; mu++){
      cudaMemsetAsync(Sigma[mu][0],0, kvol*sizeof(Complex_f),streams[mu]);
      cudaMemsetAsync(Sigma[mu][1],0, kvol*sizeof(Complex_f),streams[mu]);
      for(unsigned short nu=0; nu<ndim; nu++)
         if(nu!=mu){
            //The @f$-\nu@f$ Staple
            Kernels::Plus_staple<<<dimGrid,dimBlock,0,streams[mu]>>>(mu, nu, iu, Sigma[mu][0], Sigma[mu][1],ut[0],ut[1]);
            Kernels::Gather<<<dimGrid,dimBlock,0,streams[mu]>>>(ush[mu][0], ut[0], kvol, id, nu);
            Kernels::Gather<<<dimGrid,dimBlock,0,streams[mu]>>>(ush[mu][1], ut[1], kvol, id, nu);
 
#if(nproc>1)
            //Prefetch to the CPU for until we get NCCL working
            //cudaMemPrefetchAsync(ush[0], kvolHalo*sizeof(Complex_f),cudaCpuDeviceId,streams[0]);
            //cudaMemPrefetchAsync(ush[1], kvolHalo*sizeof(Complex_f),cudaCpuDeviceId,streams[1]);
            CHalo_swap_dir(ush[mu][0], 1, mu, DOWN); CHalo_swap_dir(ush[mu][1], 1, mu, DOWN);
            //cudaMemPrefetchAsync(ush[0]+kvol, halo*sizeof(Complex_f),device,streams[0]);
            //cudaMemPrefetchAsync(ush[1]+kvol, halo*sizeof(Complex_f),device,streams[1]);
#endif
            //Next up, the @f$-\nu@f$ staple
            Kernels::Minus_staple<<<dimGrid,dimBlock,0,streams[mu]>>>(mu, nu, iu, id,Sigma[mu][0],Sigma[mu][1],\
                  ush[mu][0],ush[mu][1],ut[0],ut[1]);
         }
      //Now get the gauge force acting in the @f$\mu@f$ direction
      Kernels::cuGaugeForce<<<dimGrid,dimBlock,0,streams[mu]>>>(mu,Sigma[mu][0],Sigma[mu][1],dSdpi,ut[0],ut[1],beta);
   }
   for(unsigned short i=0;i<ndim;i++){
#ifdef _DEBUG
      cudaFree(Sigma[i][0]); cudaFree(Sigma[i][1]);
      cudaFree(ush[i][0]); cudaFree(ush[i][1]);
#else
      cudaFreeAsync(Sigma[i][0],streams[i]); cudaFreeAsync(Sigma[i][1],streams[i]);
      cudaFreeAsync(ush[i][0],streams[i]); cudaFreeAsync(ush[i][1],streams[i]);
#endif
   }
   cudaDeviceSynchronise();
}

References CHalo_swap_dir(), Complex_f, cudaDeviceSynchronise, Kernels::cuGaugeForce(), dimBlock, dimGrid, DOWN, Kernels::Gather(), kvol, kvolHalo, Kernels::Minus_staple(), ndim, Kernels::Plus_staple(), and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuGauge_Update() [1/2]

void cuGauge_Update	(	const double	d,
		double *	pp,
		Complex *	ut[2],
		dim3	dimGrid,
		dim3	dimBlock )

CUDA wrapper for the gauge update during the integration step of the HMC.

Parameters

[in]	d	Gauge step size
[in]	pp	Momentum field
[in,out]	ut	: Double precision gauge fields
[in]	dimGrid,dimBlock	CUDA Grid/Block dimensions

Postcondition

ut updated in place

Definition at line 298 of file cusu2hmc.cu.

                                                                                            {
   for(int mu=0;mu<ndim;mu++)
      Kernels::cuGauge_Update<<<dimGrid,dimBlock,0,streams[mu]>>>(d,pp,ut[0],ut[1],mu);
   cudaDeviceSynchronise();
}

References Complex, cudaDeviceSynchronise, Kernels::cuGauge_Update(), dimBlock, dimGrid, ndim, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuGauge_Update() [2/2]

__global__ void Kernels::cuGauge_Update	(	const double	d,
		double *	pp,
		Complex *	u11t,
		Complex *	u12t,
		int	mu )

Gauge update for the integration step of the HMC.

Parameters

[in]	d	Gauge step size
[in]	pp	Momentum field
[in,out]	u11t,u12t	Double precision gauge fields
[in]	mu	Lattice direction being updated

Postcondition

ut and ut_f updated in place

Definition at line 176 of file cusu2hmc.cu.

                                                                                               {
   const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
   const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
   const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
   const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
   const unsigned int gthreadId= blockId * bsize+bthreadId;
   for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
      //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
      const unsigned int ind = i+kvol*mu;
      double AAA = d*sqrt(pp[ind]*pp[ind]\
            +pp[i+kvol*(1*ndim+mu)]*pp[i+kvol*(1*ndim+mu)]\
            +pp[i+kvol*(2*ndim+mu)]*pp[i+kvol*(2*ndim+mu)]);
      double CCC = cos(AAA);
      double SSS = d*sin(AAA)/AAA;
      Complex a11 = CCC+I*SSS*pp[i+kvol*(2*ndim+mu)];
      Complex a12 = pp[i+kvol*(1*ndim+mu)]*SSS + I*SSS*pp[ind];
      //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[ind];
      u11t[ind] = a11*b11-a12*conj(u12t[ind]);
      u12t[ind] = a11*u12t[ind]+a12*conj(b11);
   }
}

References Complex, conj(), I, kvol, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

cuGaugeForce()

__global__ void Kernels::cuGaugeForce	(	int	mu,
		Complex_f *	Sigma11,
		Complex_f *	Sigma12,
		double *	dSdpi,
		Complex_f *	u11t,
		Complex_f *	u12t,
		float	beta )

Calculates the gauge force due to the Wilson Action at each intermediate time.

Parameters

[out]	dSdpi	The force
[in]	u11t,u12t	Gauge fields
[in]	Sigma11,Sigma12	Staples
[in]	mu	Generator index
[in]	beta	Inverse gauge coupling

Postcondition

Gauge force written to dSdpi.

Definition at line 89 of file cuforce.cu.

                                                                                                                                          {
      const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
      const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
      const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
      const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
      for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
         const unsigned int ind = i+kvolHalo*mu;
         Complex_f a11 = u11t[ind]*Sigma12[i]+u12t[ind]*conj(Sigma11[i]);
         Complex_f a12 = u11t[ind]*Sigma11[i]+conj(u12t[ind])*Sigma12[i];
         //Not worth splitting into different streams, before we get ideas...
         dSdpi[i+kvol*mu]=beta*a11.imag();
         dSdpi[i+kvol*(1*ndim+mu)]=beta*a11.real();
         dSdpi[i+kvol*(2*ndim+mu)]=beta*a12.imag();
      }
   }

References Complex_f, conj(), kvol, kvolHalo, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

Force()

int Force	(	double *	dSdpi,
		const bool	iflag,
		double	res1,
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		Complex *	ut[2],
		Complex_f *	ut_f[2],
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		double *	dk[2],
		float *	dk_f[2],
		const Complex_f	jqq,
		const float	akappa,
		const float	beta,
		const float	c_sw,
		double *	ancg )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	iflag	Invert before evaluating the force. 0 to invert, one not to. Blame FORTRAN...
[in]	res1	Conjugate gradient residue
[in]	X0	Up/down partitioned pseudofermion field
[in]	X1	Holder for the partitioned fermion field, then the inverted dield
[in]	Phi	Pseudofermion field
[in]	ut,ut_f	Double/float precision colour fields
[in]	iu,id	Lattice indices
[in]	gamval,gamval_f	Double/float precision gamma matrices rescaled by \(\kappa\)
[in]	gamin	Gamma indices
[in]	sigval,sigval_f	Double/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	dk,dk_f	Double/float \(e^{-\mu}\) and \(e^\mu\)
[in]	jqq	Diquark source
[in]	akappa	Hopping parameter
[in]	beta	Inverse gauge coupling
[in]	c_sw	Clover coefficient. If non-zero calculate the clover contribution
[in]	ancg	Counter for conjugate gradient iterations

Returns

Zero on success, integer error code otherwise

Postcondition

Force added to dSdpi

Todo

If Gauge_force has not been called then set dSdpi to zero. The rest of the force uses += and as a result will add onto the previous answer instead of giving just the fermion force at this step

\((X1=(M\dagger M)^{-1} \Phi\)

Definition at line 228 of file force.c.

                                                                                                                           {
   const char funcname[] = "Force";
#ifdef USE_GPU
   int device=-1;
   cudaGetDevice(&device);
#endif
#ifndef NO_GAUGE
   Gauge_force(dSdpi,ut_f,iu,id,beta);
#else
#endif
   if(!akappa)
      return 0;
   //X1=(M†M)^{1} Phi
   int itercg=1;
   Complex_f *clover[2];
#ifdef USE_GPU
   Complex_f *X1_f, *X2_f;
   cudaMallocAsync((void **)&X1_f,kferm2Halo*sizeof(Complex_f),streams[1]);
   cudaMallocAsync((void **)&X2_f,kferm2Halo*sizeof(Complex_f),streams[0]);
#else
   Complex_f *X1_f= (Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
   Complex_f *X2_f= (Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
#endif
   if(c_sw)
   Clover(clover,ut_f,iu,id);
 
   for(int na = 0; na<nf; na++){
#ifdef USE_GPU
#if(nproc>1) //Strided
      for(unsigned short j=0;j<nc*idirac;j++)
         cudaMemcpyAsync(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[j]);
#else
      cudaMemcpyAsync(X1,X0+na*kferm2,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,NULL);
#endif
#else
      for(unsigned short j=0;j<nc*ndirac;j++)
         memcpy(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex));
#endif
      if(!iflag){
         int itercg=1;
#ifdef USE_GPU
         Complex *smallPhi;
         cudaMallocAsync((void **)&smallPhi,kferm2*sizeof(Complex),streams[0]);
#else
         Complex *smallPhi = (Complex *)aligned_alloc(AVX,kferm2*sizeof(Complex)); 
#endif
         Fill_Small_Phi(na, smallPhi, Phi);
         Congradq(na,res1,X1,smallPhi,ut,ut_f,clover,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
               jqq,akappa,c_sw,&itercg);
#ifdef USE_GPU
         cudaFreeAsync(smallPhi,streams[0]);
#else
         free(smallPhi);
#endif
         *ancg+=itercg;
#ifdef USE_GPU
         alignas(16) const Complex blasa=2.0; alignas(16) const double blasb=-1.0;
         cublasZdscal(cublas_handle,kferm2,&blasb,(cuDoubleComplex *)(X0+na*kferm2),1);
#if(nproc>1) //strided
         for(unsigned short j=0;j<nc*ndirac;j++)
            cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&blasa,(cuDoubleComplex *)X1+j*kvolHalo,1,(cuDoubleComplex *)X0+na*kferm2+j*kvol,1);
#else
         cublasZaxpy(cublas_handle,kferm2,(cuDoubleComplex *)&blasa,(cuDoubleComplex *)X1,1,(cuDoubleComplex *)(X0+na*kferm2),1);
#endif
#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
         const Complex blasa=2.0; const Complex blasb=-1.0;
         //This is not a general BLAS Routine. BLIS and MKl support it
         //CUDA and GSL does not support it
         cblas_zaxpby(kferm2, &blasa, X1, 1, &blasb, X0+na*kferm2, 1); 
#elifdef USE_BLAS
         const Complex blasa=2.0; const double blasb=-1.0;
         cblas_zdscal(kferm2,blasb,X0+na*kferm2,1);
         for(unsigned short j=0;j<nc*ndirac;j++)
            cblas_zaxpy(kvol,&blasa,X1+j*kvolHalo,1,X0+na*kferm2+j*kvol,1);
#else
#pragma omp parallel for simd collapse(2) aligned(X0,X1:AVX)
         for(int idirac=0;idirac<ndirac;idirac++){
            for(int i=0;i<kvol;i++)
               X0[i+kvol*(0+nc*(idirac+ndirac*na))]=
                  2*X1[i+kvolHalo*(0+idirac*c)]-X0[i+kvol*(0+nc*(idirac+ndirac*na))];
            X0[i+kvol*(1+nc*(idirac+ndirac*na))]=
               2*X1[i+kvolHalo*(1+idirac*c)]-X0[i+kvol*(1+nc*(idirac+ndirac*na))];
         }
#endif
      }
#ifdef USE_GPU
      cudaDeviceSynchronise();
#endif
      //Since it has to be stridded in MPI, we have to pass kvol and nc*ndirac instead of kferm2
      ComplexConvert(X1_f,X1,kvol,true,nc*ndirac);
      Hdslash_f(X2_f,X1_f,ut_f,iu,id,gamval_f,gamin,dk_f,akappa);
      if(c_sw)
      HbyClover_f(X2_f,X1_f,clover,sigval_f,akappa,sigin,false);
      //NOTE: This was orginally two. But was changed as a test for Claude so the two appears inside Force_s and Force_t
      //It may need to be reverted back later
      alignas(8) const float blasd=1.0;
#ifdef USE_GPU
      cudaDeviceSynchronise();
#if(nproc>1)
      for(unsigned short j=0;j<nc*ndirac;j++)
         cublasCsscal(cublas_handle,kvol, &blasd, (cuComplex *)X2_f+j*kvolHalo, 1);
#else
      cublasCsscal(cublas_handle,kferm2, &blasd, (cuComplex *)X2_f, 1);
#endif
#elif defined USE_BLAS
      for(unsigned short j=0;j<nc*ndirac;j++)
         cblas_csscal(kvol, blasd, X2_f+j*kvolHalo, 1);
#else
#pragma unroll
#pragma omp parallel for simd collapse(2) aligned(X2_f:AVX)
      for(unsigned short j=0;j<nc*ndirac;j++)
         for(unsigned int i=0;i<kvol;i++)
            X2_f[i+j*kvolHalo]*=1;
#endif
#if(npx>1)
      CHalo_swap_dir(X1_f,nc*ndirac,0,DOWN);
      CHalo_swap_dir(X2_f,nc*ndirac,0,DOWN);
#endif
#if(npy>1)
      CHalo_swap_dir(X1_f,nc*ndirac,1,DOWN);
      CHalo_swap_dir(X2_f,nc*ndirac,1,DOWN);
#endif
#if(npz>1)
      CHalo_swap_dir(X1_f,nc*ndirac,2,DOWN);
      CHalo_swap_dir(X2_f,nc*ndirac,2,DOWN);
#endif
#if(npt>1)
      CHalo_swap_dir(X1_f,nc*ndirac,3,DOWN);
      CHalo_swap_dir(X2_f,nc*ndirac,3,DOWN);
#endif
 
      // The original FORTRAN Comment:
      //    dSdpi=dSdpi-Re(X1*(d(Mdagger)dp)*X2) -- Yikes!
      //   we're gonna need drugs for this one......
      //
      //  Makes references to X1(.,.,iu(i,mu)) AND X2(.,.,iu(i,mu))
      //  as a result, need to swap the DOWN halos in all dirs for
      //  both these arrays, each of which has 8 cpts
      //
#ifdef USE_GPU
      cuForce(dSdpi,ut_f,X1_f,X2_f,gamval_f,dk_f,iu,gamin,akappa,dimGrid,dimBlock);
      cudaDeviceSynchronise();
#else
      //Thankfully the CUDA version is much neater so we're using that style going forwards
      for(unsigned short mu=0;mu<ndim-1;mu++)
         Force_s(dSdpi,ut_f,X1_f,X2_f,gamval_f,iu,gamin,akappa,mu);
      Force_t(dSdpi,ut_f,X1_f,X2_f,gamval_f,dk_f,iu,gamin,akappa);
#endif
         if(c_sw){
      //Clover_Force(dSdpi,ut_f,X1_f,X2_f,sigval_f,sigin,iu,id,akappa);
      Clov_Force(dSdpi,ut_f,X1_f,X2_f,sigval_f,sigin,iu,id,akappa);
      Clover_free(clover);
      }
   }
#ifdef USE_GPU
   cudaFreeAsync(X1_f,streams[0]); cudaFreeAsync(X2_f,streams[1]);
#else
   free(X1_f); free(X2_f);
#endif
   return 0;
}

References AVX, CHalo_swap_dir(), Clov_Force(), Clover(), Clover_free(), Complex, Complex_f, ComplexConvert(), Congradq(), cublas_handle, cudaDeviceSynchronise, cuForce(), dimBlock, dimGrid, DOWN, Fill_Small_Phi(), Force_s(), Force_t(), Gauge_force(), HbyClover_f(), Hdslash_f(), kferm2, kferm2Halo, kvol, kvolHalo, nc, ndim, ndirac, nf, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

Force_s()

void Force_s	(	double *	dSdpi,
		Complex_f *	ut[2],
		Complex_f *	X1,
		Complex_f *	X2,
		Complex_f	gamval[20],
		unsigned int *	iu,
		const unsigned short	gamin[16],
		const float	akappa,
		const unsigned short	mu )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	ut	Float precision colour fields
[in]	X1	Inverted field
[in]	X2	\(MX_1\)
[in]	gamval	Gamma matrices rescaled by \(\kappa\)
[in]	iu	Lattice indices
[in]	gamin	Gamma indices
[in]	akappa	Hopping parameter
[in]	mu	Force direction

Postcondition

Force added to dSdpi

Definition at line 81 of file force.c.

                                                                                                   {
#pragma omp parallel for simd
   for(unsigned int i=0;i<kvol;i++){
      const unsigned int ind=i+kvolHalo*mu;
      const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
      const unsigned int uid = iu[i+kvol*mu];
      //Similarly to Hdslash we always see idirac*nc so we do that here too.
      for(unsigned short idirac=0;idirac<nc*ndirac;idirac+=nc){
         Complex_f X1s[nc];    Complex_f X1su[nc];
         Complex_f X2s[nc];    Complex_f X2su[nc];
 
         X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
         X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
         X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
         X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
 
         float dSdpis[3];
         dSdpis[0]=dSdpi[i+kvol*mu];
         //Multiplying by i and taking the real component is the same as taking the negative imaginary component
         //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
         //generators by hand
         dSdpis[0]+=-akappa*cimag(
               conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
               +conj(X1s[1])*(u11s*X2su[0]+u12s*X2su[1])
               +conj(X1su[0])*(u12s*X2s[0]-conj(u11s)*X2s[1])
               +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1]));
 
         dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
         dSdpis[1]+=akappa*creal(
               (conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                +conj(X1s[1])*(-u11s*X2su[0]-u12s*X2su[1])
                +conj(X1su[0])*(-u12s*X2s[0]-conj(u11s)*X2s[1])
                +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1])));
 
         dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
         dSdpis[2]+=-akappa*cimag(
               conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
               +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
               +conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
               +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1]));
 
         const unsigned short gindex=mu*ndirac+(idirac>>1);
         const Complex_f gamval_c=gamval[gindex];
         //Rescaling gind by nc
         const unsigned short gind = gamin[gindex]<<1;   
         X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
         X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
 
         //If you are asked to rederive the force from Montvay and Munster you'll notice that it should be kappa*gamma
         //but below is only gamma. We rescaled gamma by kappa already when we defined it so that's where it has gone
         dSdpis[0]+=-cimag(gamval_c*
               (conj(X1s[0])* (-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                +conj(X1s[1])* (u11s *X2su[0]+u12s *X2su[1])
                +conj(X1su[0])* (-u12s *X2s[0] +conj(u11s)*X2s[1])
                +conj(X1su[1])*(u11s *X2s[0] +conj(u12s)*X2s[1])));
         dSdpi[i+kvol*mu]=dSdpis[0];
 
         dSdpis[1]+=creal(gamval_c*
               (conj(X1s[0])* (-conj(u12s)*X2su[0] +conj(u11s)*X2su[1])
                +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1])
                +conj(X1su[0])* (u12s *X2s[0]+conj(u11s)*X2s[1])
                +conj(X1su[1])* (-u11s *X2s[0]+conj(u12s)*X2s[1])));
         dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
 
         dSdpis[2]+=-cimag(gamval_c*
               (conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
                +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
                +conj(X1su[0])*(conj(u11s)*X2s[0]+u12s *X2s[1])
                +conj(X1su[1])*(conj(u12s)*X2s[0]-u11s *X2s[1])));
         dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
      }
   }
   return;
}

References cimag, Complex_f, conj(), creal, kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

Force_t()

void Force_t	(	double *	dSdpi,
		Complex_f *	ut[2],
		Complex_f *	X1,
		Complex_f *	X2,
		Complex_f	gamval[20],
		float *	dk[2],
		unsigned int *	iu,
		const unsigned short	gamin[16],
		float	akappa )

Calculates the force \(\frac{dS}{d\pi}\) at each intermediate time.

Parameters

[in,out]	dSdpi	The force
[in]	ut	Float precision colour fields
[in]	X1	Inverted field
[in]	X2	\(MX_1\)
[in]	gamval	Gamma matrices rescaled by \(\kappa\)
[in]	dk	\(e^{-\mu}\) and \(e^\mu\)
[in]	iu	Lattice indices
[in]	gamin	Gamma indices
[in]	akappa	Hopping parameter

Postcondition

Force added to dSdpi

Definition at line 156 of file force.c.

                                                                                  {
 
   const unsigned short mu=3;
#pragma omp parallel for simd
   for(unsigned int i=0;i<kvol;i++){
      const unsigned int ind=i+kvolHalo*mu;
      const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
      //TODO: The only diffrence with these is that the sign flips for the temporal components
      //       Can we figure out a way of doing this without having to read in a large array. 
      //       Will result in a conditional inside a CUDA loop. If i>kvol3
      const float dks[2] = {dk[0][i],dk[1][i]};
      //Up indices
      const unsigned int uid = iu[i+kvol*mu];
      //Similarly to Hdslash we always see idirac*nc so we do that here too.
      for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
         Complex_f X1s[nc];    Complex_f X1su[nc];
         Complex_f X2s[nc];    Complex_f X2su[nc];
 
         X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
         X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
         X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
         X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
 
         float dSdpis[3];
         dSdpis[0]=dSdpi[i+kvol*mu];
         //Multiplying by i and taking the real component is the same as taking the negative imaginary component
         //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
         //generators by hand
         dSdpis[0]+=-cimag(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                  +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
               +dks[1]*(conj(X1su[0])*(+u12s*X2s[0]-conj(u11s)*X2s[1])
                  +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1])));
 
         dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
         dSdpis[1]+=creal(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                  +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1]))
               +dks[1]*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
                  +conj(X1su[1])*( u11s *X2s[0]-conj(u12s)*X2s[1])));
 
         dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
         dSdpis[2]+=-cimag(dks[0]* (conj(X1s[0])* (u11s *X2su[0]+u12s *X2su[1])
                  +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
               +dks[1]*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
                  +conj(X1su[1])* (-conj(u12s)*X2s[0]+u11s *X2s[1])));
 
         const unsigned short gindex=mu*ndirac+(idirac>>1);
         //Rescaling gind by nc
         const unsigned short gind = gamin[gindex]<<1;   
         X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
         X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
 
         dSdpis[0]+=-cimag(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                  +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
               -dks[1]*(conj(X1su[0])* (u12s *X2s[0]-conj(u11s)*X2s[1])
                  +conj(X1su[1])*(-u11s *X2s[0]-conj(u12s)*X2s[1])));
         dSdpi[i+kvol*mu]=dSdpis[0];
 
         dSdpis[1]+=creal(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
                  +conj(X1s[1])*(-u11s*X2su[0]-u12s *X2su[1]))
               -dks[1]*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
                  +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1])));
         dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
 
         dSdpis[2]+=-cimag(dks[0]*(conj(X1s[0])*(u11s*X2su[0] +u12s *X2su[1])
                  +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
               -dks[1]*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
                  +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1])));
         dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
      }
   }
}

References cimag, Complex_f, conj(), creal, kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

Gauge_force()

int Gauge_force	(	double *	dSdpi,
		Complex_f *	ut[2],
		unsigned int *	iu,
		unsigned int *	id,
		float	beta )

Calculates the gauge force due to the Wilson Action at each intermediate time.

Parameters

[out]	dSdpi	The force
[in]	ut	Gauge fields
[in]	iu,id	Lattice indices
[in]	beta	Inverse gauge coupling

Returns

Zero on success, integer error code otherwise

Postcondition

Contents of dSdpi replaced by gauge force

Definition at line 11 of file force.c.

                                                                                              {
   const char funcname[] = "Gauge_force";
 
   //We define zero halos for debugging
   // #ifdef _DEBUG
   //    memset(ut[0][kvol], 0, ndim*halo*sizeof(Complex_f));  
   //    memset(ut[1][kvol], 0, ndim*halo*sizeof(Complex_f));  
   // #endif
   //Was a trial field halo exchange here at one point.
#ifdef USE_GPU
   cuGauge_force(ut,dSdpi,beta,iu,id,dimGrid,dimBlock);
   cudaDeviceSynchronise();
#else
   Complex_f *Sigma[2], *ush[2];
   Sigma[0] = (Complex_f *)aligned_alloc(AVX,kvol*sizeof(Complex_f)); 
   Sigma[1]= (Complex_f *)aligned_alloc(AVX,kvol*sizeof(Complex_f)); 
   ush[0] = (Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f)); 
   ush[1] = (Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f)); 
   //Holders for directions
   for(int mu=0; mu<ndim; mu++){
      memset(Sigma[0],0, kvol*sizeof(Complex_f));
      memset(Sigma[1],0, kvol*sizeof(Complex_f));
      for(int nu=0; nu<ndim; nu++)
         if(nu!=mu){
            //The +ν Staple
#pragma omp parallel for simd //aligned(ut[0],ut[1],Sigma[0],Sigma[1],iu:AVX)
            for(int i=0;i<kvol;i++){
               int uidm = iu[mu*kvol+i];
               int uidn = iu[nu*kvol+i];
               Complex_f   a11=ut[0][uidm+kvolHalo*nu]*conj(ut[0][uidn+kvolHalo*mu])+\
                               ut[1][uidm+kvolHalo*nu]*conj(ut[1][uidn+kvolHalo*mu]);
               Complex_f   a12=-ut[0][uidm+kvolHalo*nu]*ut[1][uidn+kvolHalo*mu]+\
                               ut[1][uidm+kvolHalo*nu]*ut[0][uidn+kvolHalo*mu];
               Sigma[0][i]+=a11*conj(ut[0][i+kvolHalo*nu])+a12*conj(ut[1][i+kvolHalo*nu]);
               Sigma[1][i]+=-a11*ut[1][i+kvolHalo*nu]+a12*ut[0][i+kvolHalo*nu];
            }
            C_gather(ush[0], ut[0], kvol, id, nu);
            C_gather(ush[1], ut[1], kvol, id, nu);
#if(nproc>1)
            CHalo_swap_dir(ush[0], 1, mu, DOWN); CHalo_swap_dir(ush[1], 1, mu, DOWN);
#endif
            //Next up, the -ν staple
#pragma omp parallel for simd //aligned(ut[0],ut[1],ush[0],ush[1],Sigma[0],Sigma[1],iu,id:AVX)
            for(int i=0;i<kvol;i++){
               int uidm = iu[mu*kvol+i];
               int didn = id[nu*kvol+i];
               //uidm is correct here
               Complex_f a11=conj(ush[0][uidm])*conj(ut[0][didn+kvolHalo*mu])-\
                             ush[1][uidm]*conj(ut[1][didn+kvolHalo*mu]);
               Complex_f a12=-conj(ush[0][uidm])*ut[1][didn+kvolHalo*mu]-\
                             ush[1][uidm]*ut[0][didn+kvolHalo*mu];
               Sigma[0][i]+=a11*ut[0][didn+kvolHalo*nu]-a12*conj(ut[1][didn+kvolHalo*nu]);
               Sigma[1][i]+=a11*ut[1][didn+kvolHalo*nu]+a12*conj(ut[0][didn+kvolHalo*nu]);
            }
         }
#pragma omp parallel for simd //aligned(ut[0],ut[1],Sigma[0],Sigma[1],dSdpi:AVX)
      for(int i=0;i<kvol;i++){
         const unsigned int ind = i+kvolHalo*mu;
         Complex_f a11 = ut[0][ind]*Sigma[1][i]+ut[1][ind]*conj(Sigma[0][i]);
         Complex_f a12 = ut[0][ind]*Sigma[0][i]+conj(ut[1][ind])*Sigma[1][i];
 
         dSdpi[i+kvol*mu]=(double)(beta*cimag(a11));
         dSdpi[i+kvol*(1*ndim+mu)]=(double)(beta*creal(a11));
         dSdpi[i+kvol*(2*ndim+mu)]=(double)(beta*cimag(a12));
      }
   }
   free(ush[0]); free(ush[1]); free(Sigma[0]); free(Sigma[1]);
#endif
   return 0;
}

References AVX, C_gather(), CHalo_swap_dir(), cimag, Complex_f, conj(), creal, cudaDeviceSynchronise, cuGauge_force(), dimBlock, dimGrid, DOWN, kvol, kvolHalo, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

Gauge_Update()

int Gauge_Update	(	const double	d,
		double *	pp,
		Complex *	ut[2],
		Complex_f *	ut_f[2] )

Gauge update for the integration step of the HMC.

Parameters

[in]	d	Gauge step size
[in]	pp	Momentum field
[in,out]	ut	Double precision gauge fields
[in,out]	ut_f	Single precision gauge fields

Returns

Zero on success, integer error code otherwise.

Postcondition

ut and ut_f updated in place

Definition at line 37 of file integrate.c.

                                                                               {
   char funcname[] = "Gauge_Update"; 
#ifdef USE_GPU
   cuGauge_Update(d,pp,ut,dimGrid,dimBlock);
#else
#pragma omp parallel for simd collapse(2) aligned(pp:AVX) 
   for(unsigned int i=0;i<kvol;i++)
      for(unsigned short 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
          */
         const unsigned int ind = i+kvol*mu;
         double AAA = d*sqrt(pp[ind]*pp[ind]\
               +pp[i+kvol*(1*ndim+mu)]*pp[i+kvol*(1*ndim+mu)]\
               +pp[i+kvol*(2*ndim+mu)]*pp[i+kvol*(2*ndim+mu)]);
         double CCC = cos(AAA);
         double SSS = d*sin(AAA)/AAA;
         Complex a11 = CCC+I*SSS*pp[i+kvol*(2*ndim+mu)];
         Complex a12 = pp[i+kvol*(1*ndim+mu)]*SSS + I*SSS*pp[ind];
         //b11 and b12 are ut[0] and ut[1] terms, so we'll use ut[1] directly
         //but use b11 for ut[0] to prevent RAW dependency
         Complex b11 = ut[0][i+kvolHalo*mu];
         ut[0][i+kvolHalo*mu] = a11*b11-a12*conj(ut[1][i+kvolHalo*mu]);
         ut[1][i+kvolHalo*mu] = a11*ut[1][i+kvolHalo*mu]+a12*conj(b11);
      }
#endif
   Reunitarise(ut);
   //Get trial fields from accelerator for halo exchange
   Trial_Exchange(ut,ut_f);
   return 0;
}

References Complex, Complex_f, conj(), cuGauge_Update(), dimBlock, dimGrid, I, kvol, kvolHalo, 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 *	ut[2],
		Complex_f *	ut_f[2],
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk[2],
		float *	dk_f[2],
		double *	dSdpi,
		double *	pp,
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		const Complex	jqq,
		const float	beta,
		const float	akappa,
		const float	c_sw,
		const int	stepl,
		const float	dt,
		double *	ancg,
		int *	itot,
		const 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

[in,out]	ut	Double precision colour fields
[in,out]	ut_f	Single precision colour fields
[in]	X0	Up/down partitioned pseudofermion field
[in]	X1	Holder for the partitioned fermion field, then the conjugate gradient output
[in]	Phi	Pseudofermion field
[in]	dk	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]	dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
[in,out]	dSdpi	The force
[in,out]	pp	Momentum field
[in]	iu,id	Lattice indices
[in]	gamin	Gamma indices
[in]	gamval	Double precision gamma matrices rescaled by kappa
[in]	gamval_f	Single precision gamma matrices rescaled by kappa
[in]	sigval,sigval_f	Double/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	c_sw	Clover coefficient
[in]	jqq	Diquark source
[in]	akappa	Hopping parameter
[in]	beta	Inverse gauge coupling
[in]	stepl	Steps per trajectory
[in]	dt	Step size
[in]	ancg	Counter for average conjugate gradient iterations
[in]	itot	Total average conjugate gradient iterations
[in]	proby	Termination probability for random trajectory length

Returns

Zero on success, integer error code otherwise

Postcondition

ut, ut_f, pp and dSdpi updated in place throughout the integration

References Complex, and Complex_f.

Here is the caller graph for this function:

Minus_staple()

__global__ void Kernels::Minus_staple	(	const int	mu,
		const int	nu,
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	Sigma11,
		Complex_f *	Sigma12,
		Complex_f *	u11sh,
		Complex_f *	u12sh,
		Complex_f *	u11t,
		Complex_f *	u12t )

Calculates the staple in the negative \(\mu\) direction.

Parameters

[in]	mu	\(\mu\) direction
[in]	nu	\(\nu\) direction
[in]	iu,id	Upper/lower indices
[in,out]	Sigma11,Sigma12	Staple output
[in]	u11sh,u12sh	Gauge fields in \(\mu\) direction only
[in]	u11t,u12t	Gauge fields

Postcondition

Staples added to Sigma11 and Sigma12

Definition at line 55 of file cuforce.cu.

                                                                              {
      const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
      const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
      const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
      const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
      for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
         const unsigned int uidm = iu[mu*kvol+i];
         const unsigned int didn = id[nu*kvol+i];
         //uidm is correct here
         unsigned int ind=didn+kvolHalo*mu;
         Complex_f u11s=u11t[ind]; Complex_f u12s=u12t[ind];
         Complex_f a11=conj(u11sh[uidm])*conj(u11s)-\
                       u12sh[uidm]*conj(u12s);
         Complex_f a12=-conj(u11sh[uidm])*u12s-\
                       u12sh[uidm]*u11s;
         ind=didn+kvolHalo*nu;
         u11s=u11t[ind]; u12s=u12t[ind];
         Sigma11[i]+=a11*u11s-a12*conj(u12s);
         Sigma12[i]+=a11*u12s+a12*conj(u11s);
      }
   }

References Complex_f, conj(), kvol, and kvolHalo.

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

[in]	d	Step size
[in,out]	pp	Momentum field
[in]	dSdpi	Force field

Returns

Zero on success, integer error code otherwise.

Postcondition

pp update in place

Definition at line 71 of file integrate.c.

{
   const char funcname[] = "Momentum_Update";
#ifdef USE_GPU
   cublasDaxpy(cublas_handle,kmom, &d, dSdpi, 1, pp, 1);
   cudaDeviceSynchronise();
#elif defined USE_BLAS
   cblas_daxpy(kmom, d, dSdpi, 1, pp, 1);
#else
#pragma omp parallel for simd
   for(int i=0;i<kmom;i++)
      pp[i]+=d*dSdpi[i];
#endif
   return 0;
}

References cublas_handle, cudaDeviceSynchronise, and kmom.

Here is the caller graph for this function:

OMF2()

int OMF2	(	Complex *	ut[2],
		Complex_f *	ut_f[2],
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk[2],
		float *	dk_f[2],
		double *	dSdpi,
		double *	pp,
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		const Complex	jqq,
		const float	beta,
		const float	akappa,
		const float	c_sw,
		const int	stepl,
		const float	dt,
		double *	ancg,
		int *	itot,
		const float	proby )

OMF second order five step integrator.

Parameters

[in,out]	ut	Double precision colour fields
[in,out]	ut_f	Single precision colour fields
[in]	X0	Up/down partitioned pseudofermion field
[in]	X1	Holder for the partitioned fermion field, then the conjugate gradient output
[in]	Phi	Pseudofermion field
[in]	dk	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]	dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
[in,out]	dSdpi	The force
[in,out]	pp	Momentum field
[in]	iu,id	Lattice indices
[in]	gamin	Gamma indices
[in]	gamval	Double precision gamma matrices rescaled by kappa
[in]	gamval_f	Single precision gamma matrices rescaled by kappa
[in]	sigval,sigval_f	Double/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	c_sw	Clover coefficient
[in]	jqq	Diquark source
[in]	akappa	Hopping parameter
[in]	beta	Inverse gauge coupling
[in]	stepl	Steps per trajectory
[in]	dt	Step size
[in]	ancg	Counter for average conjugate gradient iterations
[in]	itot	Total average conjugate gradient iterations
[in]	proby	Termination probability for random trajectory length

Returns

Zero on success, integer error code otherwise

Postcondition

ut, ut_f, pp and dSdpi updated in place throughout the integration

Definition at line 141 of file integrate.c.

   {
      const char funcname[] = "OMF2";
      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
      if(!rank){
         double ave_dSdpi[3];
         Force_debug(ave_dSdpi,dSdpi);
         printf("Before evaluating force on rank %i, dSdpi[0] %e dSdpi[1] %e dSdpi[2] %e\n", rank,ave_dSdpi[0],ave_dSdpi[1],ave_dSdpi[2]);
         double ave_pp[3];
         Force_debug(ave_pp,pp);
         printf("Average momentum pp[0] %e pp[1] %e pp[2] %e\n", ave_pp[0],ave_pp[1],ave_pp[2]);
      }
#endif
      Force(dSdpi, 1, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,beta,c_sw,ancg);
#ifdef _DEBUG
      if(!rank){
         double ave_dSdpi[3];
         Force_debug(ave_dSdpi,dSdpi);
         printf("Average force on rank %i, dSdpi[0] %e dSdpi[1] %e dSdpi[2] %e\n", rank,ave_dSdpi[0],ave_dSdpi[1],ave_dSdpi[2]);
      }
#endif
      //Initial momentum update
      Momentum_Update(dp,dSdpi,pp);
#ifdef _DEBUG
      if(!rank){
         double ave_pp[3];
         Force_debug(ave_pp,pp);
         printf("Average initial momentum on rank %i, pp[0] %e pp[1] %e pp[2] %e\n", rank,ave_pp[0],ave_pp[1],ave_pp[2]);
      }
#endif
 
      //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,ut,ut_f);
 
         //Calculate force for middle momentum update
         Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,beta,c_sw,ancg);
#ifdef _DEBUG
         if(!rank){
            double ave_dSdpi[3];
            Force_debug(ave_dSdpi,dSdpi);
            printf("First force update on rank %i, dSdpi[0] %e dSdpi[1] %e dSdpi[2] %e\n", rank,ave_dSdpi[0],ave_dSdpi[1],ave_dSdpi[2]);
         }
#endif
         //Now do the middle momentum update
         Momentum_Update(dpm,dSdpi,pp);
#ifdef _DEBUG
         if(!rank){
            double ave_pp[3];
            Force_debug(ave_pp,pp);
            printf("Average middle momentum on rank %i, pp[0] %e pp[1] %e pp[2] %e\n", rank,ave_pp[0],ave_pp[1],ave_pp[2]);
         }
#endif
 
         //Second gauge update
         Gauge_Update(dU,pp,ut,ut_f);
 
         //Calculate force for second momentum update
         Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,beta,c_sw,ancg);
#ifdef _DEBUG
         if(!rank){
            double ave_dSdpi[3];
            Force_debug(ave_dSdpi,dSdpi);
            printf("Second force update on rank %i, dSdpi[0] %e dSdpi[1] %e dSdpi[2] %e\n", rank,ave_dSdpi[0],ave_dSdpi[1],ave_dSdpi[2]);
         }
#endif
 
         //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);
#ifdef _DEBUG
            if(!rank){
               double ave_pp[3];
               Force_debug(ave_pp,pp);
               printf("Average final momentum on rank %i, pp[0] %e pp[1] %e pp[2] %e\n", rank,ave_pp[0],ave_pp[1],ave_pp[2]);
            }
#endif
            *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);
#ifdef _DEBUG
            if(!rank){
               double ave_pp[3];
               Force_debug(ave_pp,pp);
               printf("Average intermediate momentum on rank %i, pp[0] %e pp[1] %e pp[2] %e\n", rank,ave_pp[0],ave_pp[1],ave_pp[2]);
            }
#endif
            step++;
         }
      }while(!end_traj);
         return 0;
      }

References Complex, Complex_f, Force(), Force_debug(), Gauge_Update(), Momentum_Update(), rank, and rescgg.

Here is the call graph for this function:

Here is the caller graph for this function:

OMF4()

int OMF4	(	Complex *	ut[2],
		Complex_f *	ut_f[2],
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		double *	dk[2],
		float *	dk_f[2],
		double *	dSdpi,
		double *	pp,
		unsigned int *	iu,
		unsigned int *	id,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		const unsigned short	gamin[16],
		Complex *	sigval,
		Complex_f *	sigval_f,
		unsigned short *	sigin,
		const Complex	jqq,
		const float	beta,
		const float	akappa,
		const float	c_sw,
		const int	stepl,
		const float	dt,
		double *	ancg,
		int *	itot,
		const float	proby )

OMF fourth order eleven step integrator.

Parameters

[in,out]	ut	Double precision colour fields
[in,out]	ut_f	Single precision colour fields
[in]	X0	Up/down partitioned pseudofermion field
[in]	X1	Holder for the partitioned fermion field, then the conjugate gradient output
[in]	Phi	Pseudofermion field
[in]	dk	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]	dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
[in,out]	dSdpi	The force
[in,out]	pp	Momentum field
[in]	iu,id	Lattice indices
[in]	gamin	Gamma indices
[in]	gamval	Double precision gamma matrices rescaled by kappa
[in]	gamval_f	Single precision gamma matrices rescaled by kappa
[in]	sigval,sigval_f	Double/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	c_sw	Clover coefficient
[in]	jqq	Diquark source
[in]	akappa	Hopping parameter
[in]	beta	Inverse gauge coupling
[in]	stepl	Steps per trajectory
[in]	dt	Step size
[in]	ancg	Counter for average conjugate gradient iterations
[in,out]	itot	Total average conjugate gradient iterations
[in]	proby	Termination probability for random trajectory length

Returns

Zero on success, integer error code otherwise

Postcondition

ut, ut_f, pp and dSdpi updated in place throughout the integration

References Complex, Complex_f, dimBlock, and dimGrid.

Here is the caller graph for this function:

Plus_staple()

__global__ void Kernels::Plus_staple	(	const int	mu,
		const int	nu,
		unsigned int *	iu,
		Complex_f *	Sigma11,
		Complex_f *	Sigma12,
		Complex_f *	u11t,
		Complex_f *	u12t )

Calculates the staple in the positive \(\mu\) direction.

Parameters

[in]	mu	\(\mu\) direction
[in]	nu	\(\nu\) direction
[in]	iu	Upper indices
[in,out]	Sigma11,Sigma12	Staple output
[in]	u11t,u12t	Gauge fields

Postcondition

Staples added to Sigma11 and Sigma12

Definition at line 23 of file cuforce.cu.

                                                                                                                                                     {
      const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
      const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
      const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
      const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
      for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
         const unsigned int uidm = iu[mu*kvol+i];
         unsigned int indn=uidm+kvolHalo*nu;
         const unsigned int uidn = iu[nu*kvol+i];
         unsigned int indm=uidn+kvolHalo*mu;
         Complex_f   a11=u11t[indn]*conj(u11t[indm])+\
                         u12t[indn]*conj(u12t[indm]);
         Complex_f   a12=-u11t[indn]*u12t[indm]+\
                         u12t[indn]*u11t[indm];
         indn=i+kvolHalo*nu;
         Sigma11[i]+=a11*conj(u11t[indn])+a12*conj(u12t[indn]);
         Sigma12[i]+=-a11*u12t[indn]+a12*u11t[indn];
      }
   }

References Complex_f, conj(), kvol, and kvolHalo.

Here is the call graph for this function:

Here is the caller graph for this function: