Bosonic Observables

Collaboration diagram for Bosonic Observables:

Functions
int	Average_Plaquette (double *hg, double *avplaqs, double *avplaqt, Complex_f *ut[2], unsigned int *iu, float beta)
	Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.
int	SU2plaq (Complex_f *ut[2], Complex_f Sigma[2], unsigned int *iu, int i, int mu, int nu)
	Calculates the plaquette at site i in the \(\mu--\nu\) direction.
double	Polyakov (Complex_f *ut[2])
	Calculate the Polyakov loop (no prizes for guessing that one...).
void	cuAverage_Plaquette (double *hgs, double *hgt, Complex_f *u11t, Complex_f *u12t, unsigned int *iu, dim3 dimGrid, dim3 dimBlock)
	Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.
void	cuPolyakov (Complex_f *Sigma[2], Complex_f *ut[2], dim3 dimGrid, dim3 dimBlock)
	Calculate the Polyakov loop (no prizes for guessing that one...).
__device__ void	Device::cuSU2plaq (Complex_f *u11t, Complex_f *u12t, Complex_f *Sigma11, Complex_f *Sigma12, unsigned int *iu, const unsigned int i, const unsigned short mu, const unsigned short nu)
	Calculates the SU2 plaquette.
__global__ void	Kernels::Average_Plaquette (float *hgs_d, float *hgt_d, Complex_f *u11t, Complex_f *u12t, unsigned int *iu)
	Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.
__global__ void	Kernels::Polyakov (Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11t, Complex_f *u12t)
	Calculate the Polyakov loop (no prizes for guessing that one...).

Detailed Description

Function Documentation

Average_Plaquette() [1/2]

int Average_Plaquette	(	double *	hg,
		double *	avplaqs,
		double *	avplaqt,
		Complex_f *	ut[2],
		unsigned int *	iu,
		float	beta )

Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.

Parameters

[out]	hg	Gauge component of Hamilton
[out]	avplaqs	Average spacial Plaquette
[out]	avplaqt	Average Temporal Plaquette
[in]	ut	The trial fields
[in]	iu	Upper halo indices
[in]	beta	Inverse gauge coupling

Returns

Zero on success, integer error code otherwise

Postcondition

Contents of hg, avplaqs and avplaqt replaced with output

Definition at line 14 of file bosonic.c.

                                                                                                                   {
   const char funcname[] = "Average_Plaquette";
   /*There was a halo exchange here but moved it outside
     The FORTRAN code used several consecutive loops to get the plaquette
     Instead we'll just make the arrays variables and do everything in one loop
     Should work since in the FORTRAN Sigma11[i] only depends on i components  for example
     Since the \nu loop doesn't get called for \mu=0 we'll start at \mu=1
     */
#ifdef USE_GPU
   __managed__ double hgs = 0; __managed__ double hgt = 0;
   cuAverage_Plaquette(&hgs, &hgt, ut[0], ut[1], iu,dimGrid,dimBlock);
#else
   double hgs = 0; double hgt = 0;
   for(int mu=1;mu<ndim;mu++)
      for(int nu=0;nu<mu;nu++)
         //Don't merge into a single loop. Makes vectorisation easier?
         //Or merge into a single loop and dispense with the a arrays?
#pragma omp parallel for simd reduction(+:hgs,hgt)
         for(int i=0;i<kvol;i++){
            Complex_f Sigma[2];
            SU2plaq(ut,Sigma,iu,i,mu,nu);
            switch(mu){
               //Time component
               case(ndim-1):  hgt -= creal(Sigma[0]);
                              break;
                              //Space component
               default: hgs -= creal(Sigma[0]);
                        break;
            }
         }
#endif
#if(nproc>1)
   Par_dsum(&hgs); Par_dsum(&hgt);
#endif
   *avplaqs=-hgs/(3.0*gvol); *avplaqt=-hgt/(gvol*3.0);
   *hg=(hgs+hgt)*beta;
#ifdef _DEBUG
   if(!rank)
      printf("hgs=%e  hgt=%e  hg=%e\n", hgs, hgt, *hg);
#endif
   return 0;
}

References Complex_f, creal, cuAverage_Plaquette(), dimBlock, dimGrid, gvol, kvol, ndim, Par_dsum(), rank, and SU2plaq().

Here is the call graph for this function:

Here is the caller graph for this function:

Average_Plaquette() [2/2]

__global__ void Kernels::Average_Plaquette	(	float *	hgs_d,
		float *	hgt_d,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu )

Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.

Parameters

[out]	hgs_d,hgt_d	Gauge component of Hamilton
[in]	u11t,u12t	Gauge fields
[in]	iu	Upper halo indices

Definition at line 59 of file cubosonic.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;
      Complex_f Sigma11=0; Complex_f Sigma12=0;
      //TODO: Check if μ and ν loops inside of site loop is faster. I suspect it is due to memory locality.
      for(unsigned int i=threadId;i<kvol;i+=bsize*gsize){
         float hg_c[2];
         hg_c[0]=0; hg_c[1]=0;
 
         for(unsigned short mu=1;mu<ndim;mu++)
            for(unsigned short nu=0;nu<mu;nu++){
               //This is threadsafe as the μ and ν loops are not distributed across threads
               cuSU2plaq(u11t,u12t,&Sigma11,&Sigma12,iu,i,mu,nu);
               switch(mu){
                  //Time component
                  case(ndim-1):
                     hg_c[0] -= creal(Sigma11);
                     break;
                     //Space component
                  default:
                     hg_c[1] -=  creal(Sigma11);
                     break;
               }
            }
         hgt_d[i]=hg_c[0]; hgs_d[i]=hg_c[1];
      }
   }

References Complex_f, creal, Device::cuSU2plaq(), kvol, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

cuAverage_Plaquette()

void cuAverage_Plaquette	(	double *	hgs,
		double *	hgt,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu,
		dim3	dimGrid,
		dim3	dimBlock )

Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.

Calling wrappers.

Parameters

[out]	hgs,hgt	Gauge component of Hamilton
[in]	u11t,u12t	Gauge fields
[in]	iu	Upper halo indices
[in]	dimGrid	CUDA grid dimensions
[in]	dimBlock	CUDA block dimensions

Postcondition

Contents of hgs and hgt replaced with the results

Calculates the gauge action using new (how new?) lookup table Follows a routine called qedplaq in some QED3 code.

Definition at line 120 of file cubosonic.cu.

                                                                                                                                           {
   // float *hgs_d, *hgt_d;
   int device=-1;
   cudaGetDevice(&device);
   float *hgs_d, *hgt_d;
   //Thrust want things in a weird format for the reduction, thus we oblige
   cudaMallocAsync((void **)&hgs_d,kvol*sizeof(float),NULL);
   thrust::device_ptr<float> hgs_T = thrust::device_pointer_cast(hgs_d);
   cudaMallocAsync((void **)&hgt_d,kvol*sizeof(float),NULL);
   thrust::device_ptr<float> hgt_T = thrust::device_pointer_cast(hgt_d);
 
   Kernels::Average_Plaquette<<<dimGrid,dimBlock,0,NULL>>>(hgs_d, hgt_d, u11t, u12t, iu);
   cudaDeviceSynchronise();
 
   *hgs= (double)thrust::reduce(hgs_T,hgs_T+kvol,(float)0);
   *hgt= (double)thrust::reduce(hgt_T,hgt_T+kvol,(float)0);
   cudaDeviceSynchronise();
 
   cudaFreeAsync(hgs_d,streams[0]); cudaFreeAsync(hgt_d,streams[1]);
}

References Kernels::Average_Plaquette(), Complex_f, cudaDeviceSynchronise, dimBlock, dimGrid, kvol, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuPolyakov()

void cuPolyakov	(	Complex_f *	Sigma[2],
		Complex_f *	ut[2],
		dim3	dimGrid,
		dim3	dimBlock )

Calculate the Polyakov loop (no prizes for guessing that one...).

Parameters

[out]	Sigma	Components of the Polyakov loop
[in]	ut	The gauge fields
[in]	dimGrid	CUDA grid dimensions
[in]	dimBlock	CUDA block dimensions

Postcondition

Contents of Sigma replaced with the Polyakov loop values

Definition at line 140 of file cubosonic.cu.

                                                                                   {
   int device=-1;
   cudaGetDevice(&device);
   cudaMallocManaged((void **)&Sigma[0],kvol3*sizeof(Complex_f),cudaMemAttachGlobal);
#ifdef _DEBUG
   cudaMallocManaged((void **)&Sigma[1],kvol3*sizeof(Complex_f),cudaMemAttachGlobal);
#else
   cudaMallocAsync((void **)&Sigma[1],kvol3*sizeof(Complex_f),streams[0]);
#endif
   //Extract the time component from each site and save in corresponding Sigma
   cublasCcopy(cublas_handle,kvol3, (cuComplex *)(ut[0])+3*kvolHalo, 1, (cuComplex *)Sigma[0], 1);
   cublasCcopy(cublas_handle,kvol3, (cuComplex *)(ut[1])+3*kvolHalo, 1, (cuComplex *)Sigma[1], 1);
 
   cudaDeviceSynchronise();
   Kernels::Polyakov<<<dimGrid,dimBlock>>>(Sigma[0],Sigma[1],ut[0],ut[1]);
   //cudaMemPrefetchAsync(Sigma[0],kvol3*sizeof(Complex_f),cudaCpuDeviceId,streams[0]);
#ifdef _DEBUG
   cudaFree(Sigma[1]);
#else
   cudaFreeAsync(Sigma[1],streams[1]);
#endif
   cudaDeviceSynchronise();
}

References Complex_f, cublas_handle, cudaDeviceSynchronise, dimBlock, dimGrid, kvol3, kvolHalo, Kernels::Polyakov(), and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuSU2plaq()

__device__ void Device::cuSU2plaq	(	Complex_f *	u11t,
		Complex_f *	u12t,
		Complex_f *	Sigma11,
		Complex_f *	Sigma12,
		unsigned int *	iu,
		const unsigned int	i,
		const unsigned short	mu,
		const unsigned short	nu )

Calculates the SU2 plaquette.

Parameters

[in]	u11t,u12t	Gauge fields
[out]	Sigma11,Sigma12	Plaquette entries
[in]	iu	Site indices in the up direction
[in]	i	Site
[in]	mu,nu	Plaquette direction

Postcondition

Contents of Sigma11 and Sigma12 replaced with plaquettes.

Definition at line 26 of file cubosonic.cu.

                                                                                {
      const unsigned int uidm = iu[i+kvol*mu]; 
      unsigned int ind=i+kvol*mu;
      //Need a second index in the nu direction for the first step
      unsigned int indn=uidm+kvol*nu;
      *Sigma11=u11t[ind]*u11t[indn]-u12t[ind]*conj(u12t[indn]);
      *Sigma12=u11t[ind]*u12t[indn]+u12t[ind]*conj(u11t[indn]);
 
      const int uidn = iu[i+kvol*nu]; 
      ind=uidn+kvol*mu;
      Complex_f a11=*Sigma11*conj(u11t[ind])+*Sigma12*conj(u12t[ind]);
      Complex_f a12=-*Sigma11*u12t[ind]+*Sigma12*u11t[ind];
 
      ind=i+kvol*nu;
      *Sigma11=a11*conj(u11t[ind])+a12*conj(u12t[ind]);
      *Sigma12=-a11*u12t[ind]+a12*u11t[ind];
      return;
   }

References Complex_f, conj(), and kvol.

Here is the call graph for this function:

Here is the caller graph for this function:

Polyakov() [1/2]

__global__ void Kernels::Polyakov	(	Complex_f *	Sigma11,
		Complex_f *	Sigma12,
		Complex_f *	u11t,
		Complex_f *	u12t )

Calculate the Polyakov loop (no prizes for guessing that one...).

Parameters

[out]	Sigma11,Sigma12	Components of the Polyakov loop
[in]	u11t,u12t	The gauge fields

Definition at line 97 of file cubosonic.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<kvol3;i+=gsize*bsize){
         Complex_f Sig[2]; Sig[0]=Sigma11[i]; Sig[1]=Sigma12[i];
         Complex_f u[2];
         for(unsigned int it=1;it<ksizet;it++){
            const unsigned int indexu=it*kvol3+i;
            u[0]=u11t[indexu+3*kvol];u[1]=u12t[indexu+3*kvol];
            Complex_f a11=Sig[0]*u[0]-Sig[1]*conj(u[1]);
            //Instead of having to store a second buffer just assign it directly
            Sig[1]=Sig[0]*u[1]+Sig[1]*conj(u[0]);
            Sig[0]=a11;
         }
         Sigma11[i]=Sig[0]; Sigma12[i]=Sig[1];
      }
   }

References Complex_f, conj(), ksizet, kvol, and kvol3.

Here is the call graph for this function:

Here is the caller graph for this function:

Polyakov() [2/2]

double Polyakov

(

Complex_f *

ut[2]

)

Calculate the Polyakov loop (no prizes for guessing that one...).

Parameters

[in]

ut

The gauge fields

Returns

Double corresponding to the polyakov loop

Definition at line 74 of file bosonic.c.

                                 {
   const char funcname[] = "Polyakov";
   double poly = 0;
   Complex_f *Sigma[2];
#ifdef USE_GPU
   cuPolyakov(Sigma,ut,dimGrid,dimBlock);
#else
   Sigma[0] = (Complex_f *)aligned_alloc(AVX,kvol3*sizeof(Complex_f));
   Sigma[1] = (Complex_f *)aligned_alloc(AVX,kvol3*sizeof(Complex_f));
 
   //Extract the time component from each site and save in corresponding Sigma
   memcpy(Sigma[0],ut[0]+3*kvolHalo,kvol3*sizeof(Complex_f));
   memcpy(Sigma[1],ut[1]+3*kvolHalo,kvol3*sizeof(Complex_f));
   /* Some Fortran commentary
      Changed this routine.
      ut[0] and ut[1] now defined as normal ie (kvol+halo,4).
      Copy of Sigma[0] and Sigma[1] is changed so that it copies
      in blocks of ksizet.
      Variable indexu also used to select correct element of ut[0] and ut[1] 
      in loop 10 below.
 
      Change the order of multiplication so that it can
      be done in parallel. Start at t=1 and go up to t=T:
      previously started at t+T and looped back to 1, 2, ... T-1
      Buffers
      There is a dependency. Can only parallelise the inner loop
      */
#pragma unroll
   for(int it=1;it<ksizet;it++)
#pragma omp parallel for simd
      for(int i=0;i<kvol3;i++){
         //Seems a bit more efficient to increment indexu instead of reassigning
         //it every single loop
         int indexu=it*kvol3+i;
         Complex_f   a11=Sigma[0][i]*ut[0][indexu+kvol*3]-Sigma[1][i]*conj(ut[1][indexu+kvol*3]);
         //Instead of having to store a second buffer just assign it directly
         Sigma[1][i]=Sigma[0][i]*ut[1][indexu+kvol*3]+Sigma[1][i]*conj(ut[0][indexu+kvol*3]);
         Sigma[0][i]=a11;
      }
#endif
 
   //Multiply this partial loop with the contributions of the other cores in the
   //Time-like dimension
   //
   //Par_tmul does nothing if there is only a single processor in the time direction. So we only compile
   //its call if it is required
#if (npt>1)
#ifdef __NVCC_
#error Par_tmul is not yet implimented in CUDA as Sigma[1] is device only memory
#endif
#ifdef _DEBUG
   printf("Multiplying with MPI\n");
#endif
   Par_tmul(Sigma[0], Sigma[1]);
   //end of #if(npt>1)
#endif
   /*Now all cores have the value for the complete Polyakov line at all spacial sites
     We need to globally sum over spacial processors but not across time as these
     are duplicates. So we zero the value for all but t=0
     This is (according to the FORTRAN code) a bit of a hack
     I will expand on this hack and completely avoid any work
     for this case rather than calculating everything just to set it to zero
     */
   if(!pcoord[3+rank*ndim])
#pragma omp parallel for simd reduction(+:poly)
      for(unsigned int i=0;i<kvol3;i++)
         poly+=creal(Sigma[0][i]);
#ifdef USE_GPU
   cudaFree(Sigma[0]);
#else
   free(Sigma[0]); 
   free(Sigma[1]);
#endif
 
#if(nproc>1)
   Par_dsum(&poly);
#endif
   poly/=gvol3;
   return poly;   
}

References AVX, Complex_f, conj(), creal, cuPolyakov(), dimBlock, dimGrid, gvol3, ksizet, kvol, kvol3, kvolHalo, ndim, Par_dsum(), pcoord, and rank.

Here is the call graph for this function:

Here is the caller graph for this function:

SU2plaq()

int SU2plaq	(	Complex_f *	ut[2],
		Complex_f	Sigma[2],
		unsigned int *	iu,
		int	i,
		int	mu,
		int	nu )

Calculates the plaquette at site i in the \(\mu--\nu\) direction.

Parameters

[in]	ut	Trial fields
[out]	Sigma	Plaquette components
[in]	i	Lattice site
[in]	iu	Upper halo indices
[in]	mu,nu	Plaquette direction. Note that mu and nu can be negative to facilitate calculating plaquettes for Clover terms. No sanity checks are conducted on them in this routine.

Returns

Zero on success, integer error code otherwise

Postcondition

Plaquettes written into Sigma

References Complex, Complex_f, dimBlock, dimGrid, hd, hu, nadj, ndim, ndirac, and streams.

Here is the caller graph for this function: