su2hmc.c File Reference

An ecclectic collection of functions used in the HMC. More...

#include <assert.h>
#include <coord.h>
#include <matrices.h>
#include <par_mpi.h>
#include <random.h>
#include <string.h>
#include <su2hmc.h>

Include dependency graph for su2hmc.c:

Go to the source code of this file.

Functions
int	Init (int istart, int ibound, int iread, float beta, float fmu, float akappa, Complex_f ajq, Complex *u[2], Complex *ut[2], Complex_f *ut_f[2], Complex *gamval, Complex_f *gamval_f, int *gamin, double *dk[2], float *dk_f[2], unsigned int *iu, unsigned int *id)
	Initialises the system.
int	Hamilton (double *h, double *s, double res2, double *pp, Complex *X0, Complex *X1, Complex *Phi, Complex_f *ut[2], unsigned int *iu, unsigned int *id, Complex_f *gamval_f, int *gamin, float *dk[2], Complex_f jqq, float akappa, float beta, double *ancgh, int traj)
	Calculate the Hamiltonian.
int	C_gather (Complex_f *x, Complex_f *y, int n, unsigned int *table, unsigned int mu)
	Extracts all the single precision gauge links in the \(\mu\) direction only.
int	Z_gather (Complex *x, Complex *y, int n, unsigned int *table, unsigned int mu)
	Extracts all the double precision gauge links in the \(\mu\) direction only.
int	Fill_Small_Phi (int na, Complex *smallPhi, Complex *Phi)
int	UpDownPart (const int na, Complex *X0, Complex *R1)

Detailed Description

An ecclectic collection of functions used in the HMC.

Definition in file su2hmc.c.

Function Documentation

C_gather()

int C_gather ( Complex_f * x, Complex_f * y, int n, unsigned int * table, unsigned int mu )

inline

Extracts all the single precision gauge links in the \(\mu\) direction only.

Parameters

x	The output
y	The gauge field for a particular colour
n	Number of sites in the gauge field. This is typically kvol
table	Table containing information on nearest neighbours. Usually id or iu
mu	Direciton we're interested in extractng

Returns

Zero on success, integer error code otherwise

Definition at line 357 of file su2hmc.c.

{
   const char *funcname = "C_gather";
   //FORTRAN had a second parameter m giving the size of y (kvol+halo) normally
   //Pointers mean that's not an issue for us so I'm leaving it out
#ifdef __NVCC__
   cuC_gather(x,y,n,table,mu,dimBlock,dimGrid);
#else
#pragma omp parallel for simd aligned (x,y,table:AVX)
   for(int i=0; i<n; i++)
      x[i]=y[table[i*ndim+mu]*ndim+mu];
#endif
   return 0;
}

References Complex_f, and ndim.

Here is the caller graph for this function:

Fill_Small_Phi()

int Fill_Small_Phi ( int na, Complex * smallPhi, Complex * Phi )

inline

Copies necessary (2*4*kvol) elements of Phi into a vector variable

Parameters

na	flavour index
smallPhi	The partitioned output
Phi	The pseudofermion field

Returns

Zero on success, integer error code otherwise

Definition at line 385 of file su2hmc.c.

{
   /*Copies necessary (2*4*kvol) elements of Phi into a vector variable
    *
    * Globals:
    * =======
    * Phi:    The source array
    * 
    * Parameters:
    * ==========
    * int na: flavour index
    * Complex *smallPhi:     The target array
    *
    * Returns:
    * =======
    * Zero on success, integer error code otherwise
    */
   const char *funcname = "Fill_Small_Phi";
   //BIG and small phi index
#ifdef __NVCC__
   cuFill_Small_Phi(na,smallPhi,Phi,dimBlock,dimGrid);
#else
#pragma omp parallel for simd aligned(smallPhi,Phi:AVX) collapse(3)
   for(int i = 0; i<kvol;i++)
      for(int idirac = 0; idirac<ndirac; idirac++)
         for(int ic= 0; ic<nc; ic++)
            //   PHI_index=i*16+j*2+k;
            smallPhi[(i*ndirac+idirac)*nc+ic]=Phi[((na*kvol+i)*ngorkov+idirac)*nc+ic];
#endif
   return 0;
}

References Complex, kvol, nc, ndirac, and ngorkov.

Here is the caller graph for this function:

Hamilton()

int Hamilton	(	double *	h,
		double *	s,
		double	res2,
		double *	pp,
		Complex *	X0,
		Complex *	X1,
		Complex *	Phi,
		Complex_f *	ut[2],
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	gamval_f,
		int *	gamin,
		float *	dk[2],
		Complex_f	jqq,
		float	akappa,
		float	beta,
		double *	ancgh,
		int	traj )

Calculate the Hamiltonian.

Parameters

h	Hamiltonian
s	Action
res2	Limit for conjugate gradient
pp	Momentum field
X0	Up/down partitioned pseudofermion field
X1	Holder for the partitioned fermion field, then the conjugate gradient output
Phi	Pseudofermion field
ut	Gauge fields (single precision)
iu,id	Lattice indices
gamval_f	Single precision gamma matrices rescaled by kappa
gamin	Gamma indices
dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
jqq	Diquark source
akappa	Hopping parameter
beta	Inverse gauge coupling
ancgh	Conjugate gradient iterations counter
traj	Calling trajectory for error reporting

Returns

Zero on success. Integer Error code otherwise.

Definition at line 245 of file su2hmc.c.

                                                                                      {
   /*
    * @brief Calculate the Hamiltonian
    * 
    *
    * @param   h:                Hamiltonian
    * @param   s:                Action
    * @param   res2:             Limit for conjugate gradient
    * @param   X0:
    * @param   X1:
    * @param   Phi:
    * @param   ut[0],ut[1]:         Gauge fields
    * @param   ut[0],ut[1]:      Gauge fields
    * @param   iu,id:            Lattice indices
    * @param   gamval_f:         Gamma matrices
    * @param   gamin:            Gamma indices
    * @param   dk_f[0]:          $exp(-\mu)$ float
    * @param   dk_f[1]:          $exp(\mu)$ float
    * @param   jqq:              Diquark source
    * @param   akappa:           Hopping parameter
    * @param   beta:             Inverse gauge coupling
    * @param   ancgh:            Conjugate gradient iterations counter 
    *
    * @return  Zero on success. Integer Error code otherwise.
    */   
   const char *funcname = "Hamilton";
   //Iterate over momentum terms.
#ifdef __NVCC__
   double hp;
   int device=-1;
   cudaGetDevice(&device);
   cudaMemPrefetchAsync(pp,kmom*sizeof(double),device,NULL);
   cublasDnrm2(cublas_handle, kmom, pp, 1,&hp);
   hp*=hp;
#elif defined USE_BLAS
   double hp = cblas_dnrm2(kmom, pp, 1);
   hp*=hp;
#else
   double hp=0;
   for(int i = 0; i<kmom; i++)
      hp+=pp[i]*pp[i]; 
#endif
   hp*=0.5;
   double avplaqs, avplaqt;
   double hg = 0;
   //avplaq? isn't seen again here.
   Average_Plaquette(&hg,&avplaqs,&avplaqt,ut,iu,beta);
 
   double hf = 0; int itercg = 0;
#ifdef __NVCC__
   Complex *smallPhi;
   cudaMallocAsync((void **)&smallPhi,kferm2*sizeof(Complex),NULL);
#else
   Complex *smallPhi = aligned_alloc(AVX,kferm2*sizeof(Complex));
#endif
   //Iterating over flavours
   for(int na=0;na<nf;na++){
#ifdef __NVCC__
#ifdef _DEBUG
      cudaDeviceSynchronise();
#endif
      cudaMemcpyAsync(X1,X0+na*kferm2,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[0]);
#ifdef _DEBUG
      cudaDeviceSynchronise();
#endif
#else
      memcpy(X1,X0+na*kferm2,kferm2*sizeof(Complex));
#endif
      Fill_Small_Phi(na, smallPhi, Phi);
      if(Congradq(na,res2,X1,smallPhi,ut,iu,id,gamval_f,gamin,dk,jqq,akappa,&itercg))
         fprintf(stderr,"Trajectory %d\n", traj);
 
      *ancgh+=itercg;
#ifdef __NVCC__
      cudaMemcpyAsync(X0+na*kferm2,X1,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[0]);
#else
      memcpy(X0+na*kferm2,X1,kferm2*sizeof(Complex));
#endif
      Fill_Small_Phi(na, smallPhi,Phi);
#ifdef __NVCC__
      Complex dot;
      cublasZdotc(cublas_handle,kferm2,(cuDoubleComplex *)smallPhi,1,(cuDoubleComplex *) X1,1,(cuDoubleComplex *) &dot);
      hf+=creal(dot);
#elif defined USE_BLAS
      Complex dot;
      cblas_zdotc_sub(kferm2, smallPhi, 1, X1, 1, &dot);
      hf+=creal(dot);
#else
      //It is a dot product of the flattened arrays, could use
      //a module to convert index to coordinate array...
      for(int j=0;j<kferm2;j++)
         hf+=creal(conj(smallPhi[j])*X1[j]);
#endif
   }
#ifdef __NVCC__
   cudaFreeAsync(smallPhi,NULL);
#else
   free(smallPhi);
#endif
   //hg was summed over inside of Average_Plaquette.
#if(nproc>1)
   Par_dsum(&hp); Par_dsum(&hf);
#endif
   *s=hg+hf; *h=(*s)+hp;
#ifdef _DEBUG
   if(!rank)
      printf("hg=%.5e; hf=%.5e; hp=%.5e; h=%.5e\n", hg, hf, hp, *h);
#endif
   return 0;
}

References Average_Plaquette(), AVX, Complex, Complex_f, Congradq(), Fill_Small_Phi(), kferm2, kmom, nf, Par_dsum(), and rank.

Here is the call graph for this function:

Here is the caller graph for this function:

Init()

int Init	(	int	istart,
		int	ibound,
		int	iread,
		float	beta,
		float	fmu,
		float	akappa,
		Complex_f	ajq,
		Complex *	u[2],
		Complex *	ut[2],
		Complex_f *	ut_f[2],
		Complex *	gamval,
		Complex_f *	gamval_f,
		int *	gamin,
		double *	dk[2],
		float *	dk_f[2],
		unsigned int *	iu,
		unsigned int *	id )

Initialises the system.

Parameters

istart	Zero for cold, >1 for hot, <1 for none
ibound	Periodic boundary conditions
iread	Read configuration from file
beta	Inverse gauge coupling
fmu	Chemical potential
akappa	Hopping parameter
ajq	Diquark source
u	Gauge fields
ut	Trial gauge field
ut_f	Trial gauge field (single precision)
dk	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)^\mu\)
dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)^\mu\) float
iu,id	Up halo indices
gamin	Gamma matrix indices
gamval	Double precision gamma matrices rescaled by kappa
gamval_f	Single precision gamma matrices rescaled by kappa

Returns

Zero on success, integer error code otherwise

Definition at line 19 of file su2hmc.c.

                                                                                    {
   /*
    * Initialises the system
    *
    * Calls:
    * ======
    * Addrc, Check_addr, ran2, DHalo_swap_dir, Par_sread, Par_ranset, Reunitarise
    *
    * Globals:
    * =======
    * Complex gamval:      Gamma Matrices
    * Complex_f gamval_f:  Float Gamma matrices:
    *
    * Parameters:
    * ==========
    * int istart:          Zero for cold, >1 for hot, <1 for none
    * int ibound:          Periodic boundary conditions
    * int iread:           Read configuration from file
    * float beta:          beta
    * float fmu:           Chemical potential
    * float akappa:        Hopping parameter
    * Complex_f ajq:       Diquark source
    * Complex *u[0]:       First colour field
    * Complex *u[1]:       Second colour field
    * Complex *ut[0]:         First colour trial field
    * Complex *ut[1]:         Second colour trial field
    * Complex_f *ut_f[0]:  First float trial field
    * Complex_f *ut_f[1]:  Second float trial field
    * double   *dk[0]         $exp(-\mu)$
    * double   *dk[1]:     $exp(\mu)$
    * float    *dk_f[0]:      $exp(-\mu)$ float
    * float    *dk_f[1]:      $exp(\mu)$ float
    * unsigned int *iu:    Up halo indices
    * unsigned int *id:    Down halo indices
    *
    * Returns:
    * =======
    * Zero on success, integer error code otherwise
    */
   const char *funcname = "Init";
 
#ifdef _OPENMP
   omp_set_num_threads(nthreads);
#ifdef __INTEL_MKL__
   mkl_set_num_threads(nthreads);
#endif
#endif
   //First things first, calculate a few constants for coordinates
   Addrc(iu, id);
   //And confirm they're legit
   Check_addr(iu, ksize, ksizet, 0, kvol+halo);
   Check_addr(id, ksize, ksizet, 0, kvol+halo);
#ifdef _DEBUG
   printf("Checked addresses\n");
#endif
   double chem1=exp(fmu); double chem2 = 1/chem1;
   //CUDA this. Only limit will be the bus speed
#pragma omp parallel for simd //aligned(dk[0],dk[1]:AVX)
   for(int i = 0; i<kvol; i++){
      dk[1][i]=akappa*chem1;
      dk[0][i]=akappa*chem2;
   }
   //Anti periodic Boundary Conditions. Flip the terms at the edge of the time
   //direction
   if(ibound == -1 && pcoord[3+ndim*rank]==npt-1){
#ifdef _DEBUG
      printf("Implementing antiperiodic boundary conditions on rank %i\n", rank);
#endif
#pragma omp parallel for simd //aligned(dk[0],dk[1]:AVX)
      for(int k= kvol-1; k>=kvol-kvol3; k--){
         //int k = kvol - kvol3 + i;
         dk[1][k]*=-1;
         dk[0][k]*=-1;
      }
   }
   //These are constant so swap the halos when initialising and be done with it
   //May need to add a synchronisation statement here first
#if(npt>1)
   DHalo_swap_dir(dk[1], 1, 3, UP);
   DHalo_swap_dir(dk[0], 1, 3, UP);
#endif
   //Float versions
#ifdef __NVCC__
   cuReal_convert(dk_f[1],dk[1],kvol+halo,true,dimBlock,dimGrid);
   cuReal_convert(dk_f[0],dk[0],kvol+halo,true,dimBlock,dimGrid);
#else
#pragma omp parallel for simd //aligned(dk[0],dk[1],dk_f[0],dk_f[1]:AVX)
   for(int i=0;i<kvol+halo;i++){
      dk_f[1][i]=(float)dk[1][i];
      dk_f[0][i]=(float)dk[0][i];
   }
#endif
   //What row of each dirac/sigma matrix contains the entry acting on element i of the spinor
   int __attribute__((aligned(AVX))) gamin_t[4][4] =  {{3,2,1,0},{3,2,1,0},{2,3,0,1},{2,3,0,1}};
   //Gamma Matrices in Chiral Representation
   //See Appendix 8.1.2 of Montvay and Munster
   //_t is for temp. We copy these into the real gamvals later
#ifdef __NVCC__
   cudaMemcpy(gamin,gamin_t,4*4*sizeof(int),cudaMemcpyDefault);
#else
   memcpy(gamin,gamin_t,4*4*sizeof(int));
#endif
   //Each row of the dirac matrix contains only one non-zero entry, so that's all we encode here
   Complex __attribute__((aligned(AVX))) gamval_t[5][4] =  {{-I,-I,I,I},{-1,1,1,-1},{-I,I,I,-I},{1,1,1,1},{1,1,-1,-1}};
   //Each gamma matrix is rescaled by akappa by flattening the gamval array
#if defined USE_BLAS
   //Don't cuBLAS this. It is small and won't saturate the GPU. Let the CPU handle
   //it and just copy it later
   cblas_zdscal(5*4, akappa, gamval_t, 1);
#else
#pragma omp parallel for simd collapse(2) aligned(gamval,gamval_f:AVX)
   for(int i=0;i<5;i++)
      for(int j=0;j<4;j++)
         gamval_t[i][j]*=akappa;
#endif
 
 
#ifdef __NVCC__
   cudaMemcpy(gamval,gamval_t,5*4*sizeof(Complex),cudaMemcpyDefault);
   cuComplex_convert(gamval_f,gamval,20,true,dimBlockOne,dimGridOne);   
#else
   memcpy(gamval,gamval_t,5*4*sizeof(Complex));
   for(int i=0;i<5*4;i++)
      gamval_f[i]=(Complex_f)gamval[i];
#endif
 
#ifdef __CLOVER__
   int __attribute__((aligned(AVX))) sigin_t[7][4] =  {{0,1,2,3},{0,1,2,3},{1,0,3,2},{1,0,3,2},{1,0,3,2},{1,0,3,2},{0,1,2,3}};
   //The sigma matrices are the commutators of the gamma matrices. These are antisymmetric when you swap the indices
   //Zero is the identical index case.
   //1 is sigma_1,2
   //2 is sigma_1,3
   //3 is sigma_1,4
   //4 is sigma_2,3
   //5 is sigma_2,4
   //6 is sigma_3,4
   //TODO: Do we mutiply by 1/2 and kappa here or not? I say yes to be consistent with gamma. It means the factor of 2
   //in the sigma matrices get's droppeed here
   Complex __attribute__((aligned(AVX))) sigval_t[7][4] =  {{0,0,0,0},{-1,1,-1,1},{-I,I,-I,I},{1,1,-1,-1},{-1,-1,-1,-1},{-I,I,I,-I},{1,-1,-1,1}};
#if defined USE_BLAS
   cblas_zdscal(6*4, akappa, sigval_t+4*sizeof(Complex), 1);
#else
#pragma omp parallel for simd collapse(2) aligned(sigval,sigval_f:AVX)
   for(int i=1;i<7;i++)
      for(int j=0;j<4;j++)
         sigval_t[i][j]*=akappa;
#endif
 
#ifdef __NVCC__
   cudaMemcpy(sigval,sigval_t,7*4*sizeof(Complex),cudaMemcpyDefault);
   cuComplex_convert(sigval_f,sigval,28,true,dimBlockOne,dimGridOne);   
#else
   memcpy(sigval,sigval_t,7*4*sizeof(Complex));
   for(int i=0;i<7*4;i++)
      sigval_f[i]=(Complex_f)sigval[i];
#endif
#endif
 
   if(iread){
      if(!rank) printf("Calling Par_sread() for configuration: %i\n", iread);
      Par_sread(iread, beta, fmu, akappa, ajq,u[0],u[1],ut[0],ut[1]);
      Par_ranset(&seed,iread);
   }
   else{
      Par_ranset(&seed,iread);
      if(istart==0){
         //Initialise a cold start to zero
         //memset is safe to use here because zero is zero 
#pragma omp parallel for simd //aligned(ut[0]:AVX) 
         //Leave it to the GPU?
         for(int i=0; i<kvol*ndim;i++){
            ut[0][i]=1; ut[1][i]=0;
         }
      }
      else if(istart>0){
         //Ideally, we can use gsl_ranlux as the PRNG
#ifdef __RANLUX__
         for(int i=0; i<kvol*ndim;i++){
            ut[0][i]=2*(gsl_rng_uniform(ranlux_instd)-0.5+I*(gsl_rng_uniform(ranlux_instd)-0.5));
            ut[1][i]=2*(gsl_rng_uniform(ranlux_instd)-0.5+I*(gsl_rng_uniform(ranlux_instd)-0.5));
         }
         //If not, the Intel Vectorise Mersenne Twister
#elif (defined __INTEL_MKL__&&!defined USE_RAN2)
         //Good news, casting works for using a double to create random complex numbers
         vdRngUniform(VSL_RNG_METHOD_UNIFORM_STD_ACCURATE, stream, 2*ndim*kvol, ut[0], -1, 1);
         vdRngUniform(VSL_RNG_METHOD_UNIFORM_STD_ACCURATE, stream, 2*ndim*kvol, ut[1], -1, 1);
         //Last resort, Numerical Recipes' Ran2
#else
         for(int i=0; i<kvol*ndim;i++){
            ut[0][i]=2*(ran2(&seed)-0.5+I*(ran2(&seed)-0.5));
            ut[1][i]=2*(ran2(&seed)-0.5+I*(ran2(&seed)-0.5));
         }
#endif
      }
      else
         fprintf(stderr,"Warning %i in %s: Gauge fields are not initialised.\n", NOINIT, funcname);
 
#ifdef __NVCC__
      int device=-1;
      cudaGetDevice(&device);
      cudaMemPrefetchAsync(ut[0], ndim*kvol*sizeof(Complex),device,streams[0]);
      cudaMemPrefetchAsync(ut[1], ndim*kvol*sizeof(Complex),device,streams[1]);
#endif
      //Send trials to accelerator for reunitarisation
      Reunitarise(ut);
      //Get trials back
#ifdef __NVCC__
      cudaMemcpyAsync(u[0],ut[0],ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[0]);
      cudaMemPrefetchAsync(u[0], ndim*kvol*sizeof(Complex),device,streams[0]);
      cudaMemcpyAsync(u[1],ut[1],ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[1]);
      cudaMemPrefetchAsync(u[1], ndim*kvol*sizeof(Complex),device,streams[1]);
#else
      memcpy(u[0], ut[0], ndim*kvol*sizeof(Complex));
      memcpy(u[1], ut[1], ndim*kvol*sizeof(Complex));
#endif
   }
#ifdef _DEBUG
   printf("Initialisation Complete\n");
#endif
   return 0;
}

References Addrc(), AVX, Check_addr(), Complex, Complex_f, DHalo_swap_dir(), halo, ksize, ksizet, kvol, kvol3, ndim, npt, nthreads, Par_ranset(), Par_sread(), pcoord, ran2(), rank, Reunitarise(), seed, and UP.

Here is the call graph for this function:

Here is the caller graph for this function:

UpDownPart()

int UpDownPart ( const int na, Complex * X0, Complex * R1 )

inline

Definition at line 416 of file su2hmc.c.

                                                             {
#ifdef __NVCC__
   cuUpDownPart(na,X0,R1,dimBlock,dimGrid);
   cudaDeviceSynchronise();
#else
   //The only reason this was removed from the original function is for diagnostics
#pragma omp parallel for simd collapse(2) aligned(X0,R1:AVX)
   for(int i=0; i<kvol; i++)
      for(int idirac = 0; idirac < ndirac; idirac++){
         X0[((na*kvol+i)*ndirac+idirac)*nc]=R1[(i*ngorkov+idirac)*nc];
         X0[((na*kvol+i)*ndirac+idirac)*nc+1]=R1[(i*ngorkov+idirac)*nc+1];
      }
#endif
   return 0;
}

Z_gather()

int Z_gather ( Complex * x, Complex * y, int n, unsigned int * table, unsigned int mu )

inline

Extracts all the double precision gauge links in the \(\mu\) direction only.

Parameters

x	The output
y	The gauge field for a particular colour
n	Number of sites in the gauge field. This is typically kvol
table	Table containing information on nearest neighbours. Usually id or iu
mu	Direciton we're interested in extractng

Returns

Zero on success, integer error code otherwise

Definition at line 371 of file su2hmc.c.

{
   const char *funcname = "Z_gather";
   //FORTRAN had a second parameter m giving the size of y (kvol+halo) normally
   //Pointers mean that's not an issue for us so I'm leaving it out
#ifdef __NVCC__
   cuZ_gather(x,y,n,table,mu,dimBlock,dimGrid);
#else
#pragma omp parallel for simd aligned (x,y,table:AVX)
   for(int i=0; i<n; i++)
      x[i]=y[table[i*ndim+mu]*ndim+mu];
#endif
   return 0;
}

References Complex, and ndim.