Collaboration diagram for Helper functions:

Topics
	Site Indexing

Macros
#define	MIN(x, y)
	Return the max of two values.
#define	MAX(x, y)
	Return the min of two values.

Functions
void	Transpose_z (Complex *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
void	Transpose_c (Complex_f *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
void	Transpose_d (double *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
void	Transpose_f (float *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
void	Transpose_I (int *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
void	Transpose_U (unsigned int *out, const int, const int)
	In place transpose used to convert from AoS to SoA memory layout.
double	cureduce_sum_d (double *input, const unsigned int n, const unsigned short stream)
	Sum all terms in an array of doubles.
float	cureduce_sum_f (float *input, const unsigned int n, const unsigned short stream)
	Sum all terms in an array of floats.
void	cuTranspose_z (Complex *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuTranspose_c (Complex_f *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuTranspose_d (double *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuTranspose_f (float *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuTranspose_I (int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuTranspose_U (unsigned int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
	In place transpose used to convert from AoS to SoA memory layout.
void	cuMixed_Sumto (double d, float f, const unsigned int n, const dim3 dimGrid, const dim3 dimBlock)
	Add a single to a double value, and save the output in the double array For complex valued arrays, one may cast the complex<double> and complex<float> arrays to double and float arrays, and use 2N for the array length instead.
int	Init (const int istart, const int ibound, const int iread, const float beta, const float fmu, const float akappa, const Complex_f ajq, const float c_sw, Complex u[2], Complex ut[2], Complex_f ut_f[2], Complex gamval[20], Complex_f gamval_f[20], unsigned short gamin[16], double dk[2], float dk_f[2], unsigned int iu, unsigned int *id)
	Initialises the system.
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)
	Copies necessary (24kvol) elements of Phi into a vector variable.
int	UpDownPart (const unsigned int na, Complex X0, Complex R1)
	Up/Down partitioning of the pseudofermion field.
int	Reunitarise (Complex *ut[2])
	Reunitarises u11t and u12t as in conj(u11t[i])u11t[i]+conj(u12t[i])u12t[i]=1.
int	ComplexConvert (Complex_f a, Complex b, const unsigned int len, const bool dtof, const unsigned short stride)
	takes an array of complex float and double precision numbers and converts the precision
void	Init_CUDA (Complex u11t, Complex u12t, Complex gamval[20], Complex_f gamval_f[20], unsigned short gamin[16], double dk4m, double dk4p, unsigned int iu, unsigned int id)
	Initialise CUDA cuInit was taken already by CUDA (unsurprisingly).
void	cuFill_Small_Phi (const unsigned int na, Complex smallPhi, Complex Phi, dim3 dimBlock, dim3 dimGrid)
	Copies necessary (24kvol) elements of Phi into a vector variable.
void	cuComplex_convert (Complex_f a, Complex b, const unsigned int len, const bool dtof, dim3 dimBlock, dim3 dimGrid)
	takes an array of complex float and double precision numbers and converts the precision
void	cuReal_convert (float a, double b, const unsigned int len, const bool dtof, dim3 dimBlock, dim3 dimGrid)
	takes an array of real-valued float and double precision numbers and converts the precision
void	cuUpDownPart (const unsigned int na, Complex X0, Complex R1, dim3 dimBlock, dim3 dimGrid)
	Up/Down partitioning of the pseudofermion field.
void	cuReunitarise (Complex *ut[2], dim3 dimGrid, dim3 dimBlock)
	Reunitarises u11t and u12t as in conj(u11t[i])u11t[i]+conj(u12t[i])u12t[i]=1.
void	blockInit (int x, int y, int z, int t, dim3 dimBlock, dim3 dimGrid)
	Initialises the CUDA grid and block size for a given lattice.
template<typename T>
__device__ __forceinline__ T	conj (const T &z)
	Complex Conjugation.
void	Q_allocate_f (Complex_f p_f, Complex_f x1_f, Complex_f x2_f, Complex_f r_f, Complex_f **X1_f)
	Allocates memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	Q_allocate (Complex p, Complex x1, Complex *x2, Complex clover[2])
	Allocates double precision memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	P_allocate_f (Complex_f p_f, Complex_f r_f, Complex_f x1_f, Complex_f x2_f, Complex_f **xi_f)
	Allocates memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	P_allocate (Complex p, Complex r, Complex x1, Complex x2, Complex *clover[2])
	Allocates double precision memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	Q_free_f (Complex_f p_f, Complex_f x1_f, Complex_f x2_f, Complex_f r_f, Complex_f **X1_f)
	Frees memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	Q_free (Complex p, Complex x1, Complex *x2, Complex clover[2])
	Frees double precision memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	P_free_f (Complex_f p_f, Complex_f r_f, Complex_f x1_f, Complex_f x2_f, Complex_f **xi_f)
	Frees memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
void	P_free (Complex p, Complex r, Complex x1, Complex x2, Complex *clover[2])
	Frees memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.
template<typename T>
__global__ void	Kernels::Gather (T x, T y, const unsigned int n, unsigned int *table, const unsigned short mu)
	Extracts all the single precision gauge links in the \(\mu\) direction only.
template<typename T, unsigned int bsize>
__device__ void	Device::warpReduce_sum (volatile T *sdata, const unsigned int tid)
	Performs a warp reduction for sum.
template<typename T>
__global__ void	Kernels::Transpose (T out, const T in, const int fast_in, const int fast_out)
	Swaps the order of the gauge field so that it is now SoA instead of AoS and it is nice and coalesced in memory.
__global__ void	Kernels::Mixed_Sumto (double d, float f, const unsigned int n)
	Sums a float array into a double array.
template<typename T, unsigned int bsize>
__global__ void	Kernels::reduce_sum (T g_in_data, T g_out_data, const unsigned int n)
	Performs a block reduction for sum.
__global__ void	Kernels::Real_convert (float a, double b, const unsigned int len, const bool dtof)
	takes an array of real float and double precision numbers and converts the precision
__global__ void	Kernels::cuFill_Small_Phi (const unsigned int na, Complex smallPhi, Complex Phi)
	Copies necessary (24kvol) elements of Phi into a vector variable.
__global__ void	Kernels::cuUpDownPart (const unsigned int na, Complex X0, Complex R1)
	Up/Down partitioning of the pseudofermion field.
template<typename T>
__global__ void	Kernels::cuReunitarise (complex< T > u11t, complex< T > u12t)
	Reunitarises u11t and u12t as in conj(u11t[i])u11t[i]+conj(u12t[i])u12t[i]=1.

Detailed Description

Macro Definition Documentation

◆ MAX

#define MAX	(		x,
			y )

Value:

(x>y?x:y)

Return the min of two values.

Parameters

[in] x,y values being checked

Definition at line 21 of file cusu2hmc.cu.

◆ MIN

#define MIN	(		x,
			y )

Value:

(x<y?x:y)

Return the max of two values.

Parameters

[in] x,y values being checked

Definition at line 14 of file cusu2hmc.cu.

Function Documentation

◆ blockInit()

void blockInit	(	int	x,
		int	y,
		int	z,
		int	t,
		dim3 *	dimBlock,
		dim3 *	dimGrid )

Initialises the CUDA grid and block size for a given lattice.

Parameters

[in]	x,y,z,t	Lattice dimensions
[out]	dimGrid,dimBlock	CUDA grid/block size

Postcondition: dimGrid and dimBlock initialised

Definition at line 205 of file cusu2hmc.cu.

                                                                         {
 
   const char *funcname = "blockInit";
 
   int device=-1; cudaGetDevice(&device);
   cudaDeviceProp prop; cudaGetDeviceProperties(&prop, device);
   //Threads per block
   int tpb=prop.maxThreadsPerBlock/8;
   //Warp size
   int tpw=prop.warpSize;
   int bx=1;
   //Set bx to be the largest power of 2 less than x that fits in a block
   while(bx<=x/2 && bx<tpb)
      bx*=2;
   int by=1;
   //Set by to be the largest power of 2 less than y such that bx*by fits in a block
   while(by<=y/2 && bx*by<tpb)
      by*=2;
 
   if(bx*by>=128){
      *dimBlock=dim3(bx,by);
      //If the block size neatly divides the lattice size we can create
      //extra blocks safely
      int res= ((nx*ny)/(bx*by) > 1) ? (nx*ny)/(bx*by) :1;
      //    int res = 1;
      *dimGrid=dim3(nz,nt,res);
   }
   else{
      int bz=1;
      //Set by to be the largest power of 2 less than y such that bx*by fits in an optimal block
      while(bz<z/2 && bx*by*bz<tpb)
         bz*=2;
      *dimBlock=dim3(bx,by,bz);
 
      //If we have an awkward block size then flag it.
      if(bx*by*bz%tpw!=0)
         fprintf(stderr,"Alert %i in %s: Suboptimal block size for warp size %d. bx=%d by=%d bz=%d\n",
               BLOCKALERT, funcname, tpw, bx, by,bz);
      int res= ((nx*ny)/(bx*by) > 1) ? (nx*ny)/(bx*by) :1;
      *dimGrid=dim3(z/bz,nt,res);
   }
   printf("Block: (%d,%d,%d)\tGrid: (%d,%d,%d)\n",dimBlock->x,dimBlock->y,dimBlock->z,dimGrid->x,dimGrid->y,dimGrid->z);
}

References BLOCKALERT, dimBlock, dimGrid, nt, nx, ny, and nz.

Here is the caller graph for this function:

◆ 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

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

Returns: Zero on success, integer error code otherwise

Postcondition: Contents of x replaced with output

Definition at line 291 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
#pragma omp parallel for simd aligned (x,y,table:AVX)
   for(unsigned int i=0; i<n; i++)
      x[i]=y[table[i+kvol*mu]+kvolHalo*mu];
   return 0;
}

References Complex_f, kvol, and kvolHalo.

Here is the caller graph for this function:

◆ ComplexConvert()

int ComplexConvert	(	Complex_f *	a,
		Complex *	b,
		const unsigned int	len,
		const bool	dtof,
		const unsigned short	stride )

takes an array of complex float and double precision numbers and converts the precision

Parameters

[in,out]	a	Float array
[in,out]	b	Double array
[in]	len	Number of elements to convert per stride. Striding needed to handle halo terms
[in]	dtof	If true, convert double to float. Otherwise convert float to double
[in]	stride	For terms with a halo, we need to convert in blocks of len separated by (len+halo)

Returns: Zero on success, integer error code otherwise

Postcondition: Depending on the value of dtof, either the contents of a or b are overwritten with those of the other array in the opposite precision.

Definition at line 420 of file coord.c.

                                                                                                                  {
   const char funcname[] = "ComplexConvert";
   switch(stride){
      case(0):
         fprintf(stderr,"Error %i in %s: Stride of %d is not valid.\nExiting...\n\n",STRDERROR,funcname,stride);
#if (nproc>1)
         MPI_Abort(comm,STRDERROR);
#else
         exit(STRDERROR);
#endif
         break;
      case(1):
#ifdef USE_GPU
         cuComplex_convert(a,b,len*stride,dtof,dimBlock,dimGrid);
#else
         if(dtof)
#pragma omp parallel for simd aligned(a,b:AVX)
            for(unsigned int i=0;i<len*stride;i++)
               a[i]=(Complex_f)b[i];
         else
#pragma omp parallel for simd aligned(a,b:AVX)
            for(unsigned int i=0;i<len*stride;i++)
               b[i]=(Complex)a[i];
#endif
         break;
      default:
         for(unsigned short j=0;j<stride;j++){
#ifdef USE_GPU
            cuComplex_convert(a+j*(len+halo),b+j*(len+halo),len,dtof,dimBlock,dimGrid);
#else
            if(dtof)
#pragma omp parallel for simd aligned(a,b:AVX)
               for(unsigned int i=0;i<len;i++)
                  a[i+j*(len+halo)]=(Complex_f)b[i+j*(len+halo)];
            else
#pragma omp parallel for simd aligned(a,b:AVX)
               for(unsigned int i=0;i<len;i++)
                  b[i+j*(len+halo)]=(Complex)a[i+j*(len+halo)];
#endif
         }
         break;
   }
   return 0;
}

References Complex, Complex_f, cuComplex_convert(), dimBlock, dimGrid, halo, and STRDERROR.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ conj()

template<typename T>

__device__ __forceinline__ T conj ( const T & z )

Complex Conjugation.

Parameters

z	Number to be conjugated

Returns: Conjugate of z

Definition at line 33 of file cusu2hmc.cu.

                                             {
   return T(z.real(),-z.imag());
}

Here is the caller graph for this function:

◆ cuComplex_convert()

void cuComplex_convert	(	Complex_f *	a,
		Complex *	b,
		const unsigned int	len,
		const bool	dtof,
		dim3	dimBlock,
		dim3	dimGrid )

takes an array of complex float and double precision numbers and converts the precision

Parameters

[in,out]	a	Float array
[in,out]	b	Double array
[in]	len	Number of elements to convert
[in]	dtof	If true, convert double to float. Otherwise convert float to double
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition: Depending on the value of dtof, either the contents of a or b are overwritten with those of the other array in the opposite precision.

Definition at line 284 of file cusu2hmc.cu.

                                                                                                                      {
   const char *funcname = "cuComplex_convert";
   Kernels::Real_convert<<<dimGrid,dimBlock>>>((float *)a,(double *)b,2*len,dtof);
}

References Complex, Complex_f, dimBlock, dimGrid, and Kernels::Real_convert().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuFill_Small_Phi() [1/2]

void cuFill_Small_Phi	(	const unsigned int	na,
		Complex *	smallPhi,
		Complex *	Phi,
		dim3	dimBlock,
		dim3	dimGrid )

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

Parameters

[in]	na	flavour index
[out]	smallPhi	The partitioned output
[in]	Phi	The pseudofermion field
[in]	dimGrid,dimBlock	CUDA grid/block size

Definition at line 288 of file cusu2hmc.cu.

                                                                                                          {
   Kernels::cuFill_Small_Phi<<<dimGrid,dimBlock>>>(na,smallPhi,Phi);
}

References Complex, Kernels::cuFill_Small_Phi(), dimBlock, and dimGrid.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuFill_Small_Phi() [2/2]

__global__ void Kernels::cuFill_Small_Phi	(	const unsigned int	na,
		Complex *	smallPhi,
		Complex *	Phi )

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

Parameters

na	flavour index
smallPhi	The partitioned output
Phi	The pseudofermion field

Postcondition: Result written into smallPhi

Definition at line 76 of file cusu2hmc.cu.

{
   //BIG and small phi index
   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)
      for(unsigned short idirac = 0; idirac<ndirac; idirac++)
         for(unsigned short ic= 0; ic<nc; ic++)
            //   PHI_index=i*16+j*2+k;
            smallPhi[i + kvol * (ic + nc * idirac)] = Phi[i + kvol * (ic + nc * (idirac + ngorkov * na))];
}

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

Here is the caller graph for this function:

◆ cuMixed_Sumto()

void cuMixed_Sumto	(	double *	d,
		float *	f,
		const unsigned int	n,
		const dim3	dimGrid,
		const dim3	dimBlock )

Add a single to a double value, and save the output in the double array For complex valued arrays, one may cast the complex<double> and complex<float> arrays to double and float arrays, and use 2N for the array length instead.

Parameters

[in,out]	d	Double array
[in]	f	float array
[in]	n	Array lengths
[in]	dimGrid,dimBlock	CUDA grid/block

Postcondition: d now contains the sum. (i.e. result is stored in place)

Definition at line 769 of file cumatrices.cu.

                                                                                                   {
   Kernels::Mixed_Sumto<<<dimGrid,dimBlock>>>(d,f,n);
   return;
}

References dimBlock, dimGrid, and Kernels::Mixed_Sumto().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuReal_convert()

void cuReal_convert	(	float *	a,
		double *	b,
		const unsigned int	len,
		const bool	dtof,
		dim3	dimBlock,
		dim3	dimGrid )

takes an array of real-valued float and double precision numbers and converts the precision

Parameters

[in,out]	a	Float array
[in,out]	b	Double array
[in]	len	Number of elements to convert
[in]	dtof	If true, convert double to float. Otherwise convert float to double
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition: Depending on the value of dtof, either the contents of a or b are overwritten with those of the other array in the opposite precision.

Definition at line 280 of file cusu2hmc.cu.

                                                                                                              {
   const char *funcname = "cuComplex_convert";
   Kernels::Real_convert<<<dimGrid,dimBlock>>>(a,b,len,dtof);
}

References dimBlock, dimGrid, and Kernels::Real_convert().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cureduce_sum_d()

double cureduce_sum_d	(	double *	input,
		const unsigned int	n,
		const unsigned short	stream )

Sum all terms in an array of doubles.

Parameters

[in]	input	Input array
[in]	n	Number of terms
[in]	stream	What stream to use (useful for simultaneous reductions)

Returns: Sum of all terms in input

Definition at line 587 of file cumatrices.cu.

                                                                                      {
   const unsigned int bsize=256;
   unsigned int gsize=(n + (2 * bsize) - 1) / (2 * bsize); 
   double *cachein, *cacheout;
   cudaMallocAsync(&cacheout,gsize*sizeof(double),streams[stream]);
   Kernels::reduce_sum<double,bsize><<<gsize,bsize,bsize*sizeof(double),streams[stream]>>>(input,cacheout,n);
   while(gsize>1){
      cudaMallocAsync(&cachein,gsize*sizeof(double),streams[stream]);
      cudaMemcpyAsync(cachein,cacheout,gsize*sizeof(double),cudaMemcpyDefault,streams[stream]);
      cudaFreeAsync(cacheout,streams[stream]);
      gsize>>=1;
      cudaMallocAsync(&cacheout,gsize*sizeof(double),streams[stream]);
      Kernels::reduce_sum<double,bsize><<<gsize,bsize,bsize*sizeof(double),streams[stream]>>>(cachein,cacheout,gsize);
      cudaFreeAsync(cachein,streams[stream]);
   }
   double output=0;
   cudaStreamSynchronize(streams[stream]);
   cudaMemcpyAsync(&output,cacheout,sizeof(double),cudaMemcpyDefault,streams[stream]);
   cudaStreamSynchronize(streams[stream]);
   cudaFreeAsync(cacheout,streams[stream]);
   return output;
}

References Kernels::reduce_sum(), and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cureduce_sum_f()

float cureduce_sum_f	(	float *	input,
		const unsigned int	n,
		const unsigned short	stream )

Sum all terms in an array of floats.

Parameters

[in]	input	Input array
[in]	n	Number of terms
[in]	stream	What stream to use (useful for simultaneous reductions)

Returns: Sum of all terms in input

References Complex, Complex_f, dimBlock, and dimGrid.

◆ cuReunitarise() [1/2]

void cuReunitarise	(	Complex *	ut[2],
		dim3	dimGrid,
		dim3	dimBlock )

Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.

If you're looking at the FORTRAN code be careful. There are two header files for the /trial/ header. One with u11 u12 (which was included here originally) and the other with u11t and u12t.

Parameters

[in,out]	ut	Trial fields to be reunitarised
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition: ut replaced with reunitarised gauge fields

Definition at line 294 of file cusu2hmc.cu.

                                                               {
   Kernels::cuReunitarise<<<dimGrid,dimBlock>>>(ut[0],ut[1]);
   cudaDeviceSynchronise();
}

References Complex, cudaDeviceSynchronise, Kernels::cuReunitarise(), dimBlock, and dimGrid.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuReunitarise() [2/2]

template<typename T>

__global__ void Kernels::cuReunitarise	(	complex< T > *	u11t,
		complex< T > *	u12t )

Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.

If you're looking at the FORTRAN code be careful. There are two header files for the /trial/ header. One with u11 u12 (which was included here originally) and the other with u11t and u12t.

Parameters

[in,out] u11t,u12t Trial fields to be reunitarised

Postcondition: u11t and u12t replaced with reunitarised gauge fields

Definition at line 129 of file cusu2hmc.cu.

                                                                  {
   /*
    * Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1
    *
    * If you're looking at the FORTRAN code be careful. There are two header files
    * for the /trial/ header. One with u11 u12 (which was included here originally)
    * and the other with u11t and u12t.
    *
    * Globals:
    * =======
    * u11t, u12t
    *
    * Returns:
    * ========
    * Zero on success, integer error code otherwise
    */
   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*ndim; i+=gsize*bsize){
      //Declaring anorm inside the loop will hopefully let the compiler know it
      //is safe to vectorise aggessively
      double anorm=sqrt(conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]).real();
      //    Exception handling code. May be faster to leave out as the exit prevents vectorisation.
      //    if(anorm==0){
      //       fprintf(stderr, "Error %i in %s on rank %i: anorm = 0 for μ=%i and i=%i.\nExiting...\n\n",
      //             DIVZERO, funcname, rank, mu, i);
      //       MPI_Finalise();
      //       exit(DIVZERO);
      //    }
      u11t[i]/=anorm;
      u12t[i]/=anorm;
   }
}

References conj(), kvol, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_c()

void cuTranspose_c	(	Complex_f *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 732 of file cumatrices.cu.

                                                                                                                  {
   Complex_f *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(Complex_f));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(Complex_f),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaDeviceSynchronise();
   cudaFree(holder);
}

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_d()

void cuTranspose_d	(	double *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 740 of file cumatrices.cu.

                                                                                                               {
   double *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(double));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(double),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaFree(holder);
}

References dimBlock, dimGrid, and Kernels::Transpose().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_f()

void cuTranspose_f	(	float *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 747 of file cumatrices.cu.

                                                                                                              {
   float *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(float));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(float),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaFree(holder);
}

References dimBlock, dimGrid, and Kernels::Transpose().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_I()

void cuTranspose_I	(	int *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 754 of file cumatrices.cu.

                                                                                                            {
   int *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(int));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(int),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaFree(holder);
}

References dimBlock, dimGrid, and Kernels::Transpose().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_U()

void cuTranspose_U	(	unsigned int *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 761 of file cumatrices.cu.

                                                                                                                     {
   unsigned int *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(unsigned int));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(unsigned int),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaFree(holder);
}

References dimBlock, dimGrid, and Kernels::Transpose().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuTranspose_z()

void cuTranspose_z	(	Complex *	out,
		const int	fast_in,
		const int	fast_out,
		const dim3	dimGrid,
		const dim3	dimBlock )

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index
[in]	dimGrid	CUDA grid layout
[in]	dimBlock	CUDA block layout

Postcondition: Result overwrites existing out with transposed array

Definition at line 725 of file cumatrices.cu.

                                                                                                                {
   Complex *holder;
   cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(Complex));
   cudaMemcpy(holder,out,fast_in*fast_out*sizeof(Complex),cudaMemcpyDefault);
   Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
   cudaFree(holder);
}

References Complex, dimBlock, dimGrid, and Kernels::Transpose().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuUpDownPart() [1/2]

void cuUpDownPart	(	const unsigned int	na,
		Complex *	X0,
		Complex *	R1,
		dim3	dimBlock,
		dim3	dimGrid )

Up/Down partitioning of the pseudofermion field.

Parameters

[in]	na	Flavour index
[out]	X0	Partitioned field
[in]	R1	Full pseudofermion field
[in]	dimGrid,dimBlock	CUDA grid/block size

Postcondition: Result written into X0

Definition at line 291 of file cusu2hmc.cu.

                                                                                              {
   Kernels::cuUpDownPart<<<dimGrid,dimBlock>>>(na,X0,R1);   
}

References Complex, Kernels::cuUpDownPart(), dimBlock, and dimGrid.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ cuUpDownPart() [2/2]

__global__ void Kernels::cuUpDownPart	(	const unsigned int	na,
		Complex *	X0,
		Complex *	R1 )

Up/Down partitioning of the pseudofermion field.

Parameters

[in]	na	Flavour index
[out]	X0	Partitioned field
[in]	R1	Full pseudofermion field

Postcondition: Result written to X0

Definition at line 101 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;
   //Up/down partitioning (using only pseudofermions of flavour 1)
   for(unsigned int i = gthreadId; i<kvol;i+=gsize*bsize)
      for(unsigned short idirac = 0; idirac < ndirac; idirac++){
      //R1 has ngorkov spinors, but we only want the first four.
         X0[i+kvol*(0+nc*(idirac+ndirac*na))]=R1[i+kvol*(0+nc*idirac)];
         X0[i+kvol*(1+nc*(idirac+ndirac*na))]=R1[i+kvol*(1+nc*idirac)];
      }
}

References Complex, kvol, nc, and ndirac.

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

[in]	na	flavour index
[out]	smallPhi	The partitioned output
[in]	Phi	The pseudofermion field

Returns: Zero on success, integer error code otherwise

Postcondition: Result written into smallPhi

Definition at line 311 of file su2hmc.c.

{
   const char funcname[] = "Fill_Small_Phi";
   //BIG and small phi index
#ifdef USE_GPU
   cuFill_Small_Phi(na,smallPhi,Phi,dimBlock,dimGrid);
#else
#pragma omp parallel for simd aligned(smallPhi,Phi:AVX) collapse(3)
   for(unsigned int i = 0; i<kvol;i++)
      for(unsigned short idirac = 0; idirac<ndirac; idirac++)
         for(unsigned short ic= 0; ic<nc; ic++)
            //   PHI_index=i*16+j*2+k;
            smallPhi[i+kvol*(ic+nc*idirac)]=Phi[i+kvol*(ic+nc*(idirac+ngorkov*na))];
#endif
   return 0;
}

References Complex, cuFill_Small_Phi(), dimBlock, dimGrid, kvol, nc, ndirac, and ngorkov.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Gather()

template<typename T>

__global__ void Kernels::Gather	(	T *	x,
		T *	y,
		const unsigned int	n,
		unsigned int *	table,
		const unsigned short	mu )

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

Parameters

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

Definition at line 116 of file cuforce.cu.

      {
         //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
         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 int kvbmu=kvolHalo*mu;
         for(unsigned int i = gthreadId; i<kvol;i+=gsize*bsize)
            x[i]=y[table[i+kvbmu]+kvbmu];
      }

References kvol, and kvolHalo.

Here is the caller graph for this function:

◆ Init()

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

Initialises the system.

Parameters

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

Returns: Zero on success, integer error code otherwise

Postcondition: Contents of all out arguments overwritten

Definition at line 10 of file su2hmc.c.

                                            {
   const char funcname[] = "Init";
 
#ifdef _OPENMP
   omp_set_num_threads(nthreads);
#ifdef __USE_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, kvolHalo);
   Check_addr(id, ksize, ksizet, 0, kvolHalo);
#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(unsigned int i = 0; i<kvol; i++){
      dk[0][i]=akappa*chem1; dk[1][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(unsigned 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 USE_GPU
   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(unsigned 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
   unsigned short __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 USE_GPU
   cudaMemcpy(gamin,gamin_t,4*4*sizeof(short),cudaMemcpyHostToDevice);
#else
   memcpy(gamin,gamin_t,4*4*sizeof(short));
#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(unsigned short i=0;i<5;i++)
      for(unsigned short j=0;j<4;j++)
         gamval_t[i][j]*=akappa;
#endif
 
 
#ifdef USE_GPU
   cudaMemcpy(gamval,gamval_t,5*4*sizeof(Complex),cudaMemcpyHostToDevice);
   cuComplex_convert(gamval_f,gamval,20,true,dimBlockOne,dimGridOne);   
#else
   memcpy(gamval,gamval_t,5*4*sizeof(Complex));
   for(unsigned short i=0;i<5*4;i++)
      gamval_f[i]=(Complex_f)gamval[i];
#endif
 
   if(iread){
      if(!rank) printf("Calling Par_sread() for configuration: %i\n", iread);
      Par_sread(iread, beta, fmu, akappa, ajq,c_sw,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(unsigned int i=0; i<kvol;i++)
            for(unsigned short mu=0;mu<ndim;mu++){
               ut[0][i+kvolHalo*mu]=1+0*I;   ut[1][i+kvolHalo*mu]=0+0*I;
            }
      }
      else if(istart>0){
         //Ideally, we can use gsl_ranlux as the PRNG
#ifdef __RANLUX__
         for(unsigned int i=0; i<kvol;i++)
            for(unsigned short mu=0;mu<ndim;mu++){
               ut[0][i+kvolHalo*mu]=2*(gsl_rng_uniform(ranlux_instd)-0.5+I*(gsl_rng_uniform(ranlux_instd)-0.5));
               ut[1][i+kvolHalo*mu]=2*(gsl_rng_uniform(ranlux_instd)-0.5+I*(gsl_rng_uniform(ranlux_instd)-0.5));
            }
         //Last resort, Numerical Recipes' Ran2
#else
         for(unsigned int i=0; i<kvol;i++)
            for(unsigned short mu=0;mu<ndim;mu++){
               ut[0][i+kvolHalo*mu]=2*(ran2(&seed)-0.5+I*(ran2(&seed)-0.5));
               ut[1][i+kvolHalo*mu]=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 USE_GPU
      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 USE_GPU
#if (nproc>1) //Strided for multi-GPU
      for(unsigned short mu=0;mu<ndim;mu++){
         cudaMemcpy(u[0]+kvol*mu, ut[0]+kvolHalo*mu, kvol*sizeof(Complex),cudaMemcpyDefault);
         cudaMemcpy(u[1]+kvol*mu, ut[1]+kvolHalo*mu, kvol*sizeof(Complex),cudaMemcpyDefault);
      }
#else
      cudaMemcpy(u[0], ut[0], ndim*kvol*sizeof(Complex),cudaMemcpyDefault);
      cudaMemcpy(u[1], ut[1], ndim*kvol*sizeof(Complex),cudaMemcpyDefault);
#endif
#else
      for(unsigned short mu=0;mu<ndim;mu++){
         memcpy(u[0]+kvol*mu, ut[0]+kvolHalo*mu, kvol*sizeof(Complex));
         memcpy(u[1]+kvol*mu, ut[1]+kvolHalo*mu, kvol*sizeof(Complex));
      }
#endif
   }
#ifdef _DEBUG
   printf("Initialisation Complete\n");
#endif
   return 0;
}

References Addrc(), AVX, Check_addr(), Complex, Complex_f, cuComplex_convert(), cuReal_convert(), DHalo_swap_dir(), dimBlock, dimBlockOne, dimGrid, dimGridOne, halo, I, ksize, ksizet, kvol, kvol3, kvolHalo, ndim, NOINIT, npt, nthreads, Par_ranset(), Par_sread(), pcoord, rank, ranlux_instd, Reunitarise(), seed, and UP.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Init_CUDA()

void Init_CUDA	(	Complex *	u11t,
		Complex *	u12t,
		Complex	gamval[20],
		Complex_f	gamval_f[20],
		unsigned short	gamin[16],
		double *	dk4m,
		double *	dk4p,
		unsigned int *	iu,
		unsigned int *	id )

Initialise CUDA cuInit was taken already by CUDA (unsurprisingly).

Parameters

[in]	u11t,u12t	Trial gauge fields
[in]	gamval,gamval_f	Double/float precision gamma matrices rescaled by kappa
[in]	gamin	Gamma matrix indices
[in]	dk4m,dk4p	\(e^{-\mu}\) and \(e^\mu\)
[in]	iu,id	Up/lower halo indices

Todo: CUDA 13 changed how setting devices work, so it's not just an integer any more. Those lines are commented out for now. They are not critical. Only hints for unified memory management.

Definition at line 248 of file cusu2hmc.cu.

                                                       {
   const char *funcname = "Init_CUDA";
   int device=-1;
   cudaGetDevice(&device);
   //Initialise streams for concurrent kernels
   for(unsigned int i=0;i<(ndirac*nadj*ndim);i++)
      cudaStreamCreate(&streams[i]);
   //Set iu and id to mainly read in CUDA and prefetch them to the GPU
   //cudaMemPrefetchAsync(iu,ndim*kvol*sizeof(int),device,streams[0]);
   //cudaMemPrefetchAsync(id,ndim*kvol*sizeof(int),device,streams[1]);
   //cudaMemAdvise(iu,ndim*kvol*sizeof(int),cudaMemAdviseSetReadMostly,device);
   //cudaMemAdvise(id,ndim*kvol*sizeof(int),cudaMemAdviseSetReadMostly,device);
 
   //Gamma matrices and indices on the GPU
   //cudaMemcpy(gamin_d,gamin,4*4*sizeof(int),cudaMemcpyHostToDevice);
   //cudaMemcpy(gamval_d,gamval,5*4*sizeof(Complex),cudaMemcpyHostToDevice);
   //cudaMemcpy(gamval_f_d,gamval_f,5*4*sizeof(Complex_f),cudaMemcpyHostToDevice);
   //cudaMemAdvise(gamin,4*4*sizeof(int),cudaMemAdviseSetReadMostly,device);
   //cudaMemAdvise(gamval,5*4*sizeof(Complex),cudaMemAdviseSetReadMostly,device);
 
   //More prefetching and marking as read-only (mostly)
   //Prefetching Momentum Fields and Trial Fields to GPU
   //cudaMemAdvise(dk4p,(kvol+halo)*sizeof(double),cudaMemAdviseSetReadMostly,device);
   //cudaMemAdvise(dk4m,(kvol+halo)*sizeof(double),cudaMemAdviseSetReadMostly,device);
 
   //cudaMemPrefetchAsync(dk4p,(kvol+halo)*sizeof(double),device,streams[2]);
   //cudaMemPrefetchAsync(dk4m,(kvol+halo)*sizeof(double),device,streams[3]);
 
   //cudaMemPrefetchAsync(u11t, ndim*kvol*sizeof(Complex),device,streams[4]);
   //cudaMemPrefetchAsync(u12t, ndim*kvol*sizeof(Complex),device,streams[5]);
}

References Complex, Complex_f, nadj, ndim, ndirac, and streams.

Here is the caller graph for this function:

◆ Mixed_Sumto()

__global__ void Kernels::Mixed_Sumto	(	double *	d,
		float *	f,
		const unsigned int	n )

Sums a float array into a double array.

Parameters

[out]	d	The double array
[in]	f	The float array
[in]	n	The size of the arrays

Definition at line 520 of file cumatrices.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<n;i+=bsize*gsize)
         d[i]+=(double)f[i];
      return;
   }

Here is the caller graph for this function:

◆ P_allocate()

void P_allocate	(	Complex **	p,
		Complex **	r,
		Complex **	x1,
		Complex **	x2,
		Complex *	clover[2] )

Allocates double precision memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[out]	p	Holder for fermion field during inversion
[out]	r	Residue vector for conjugate gradient
[out]	x1,x2	\(M\) and \(M^\dagger M\)
[out]	clover	Clover Fields

Postcondition: Memory allocated

Definition at line 140 of file congrad.c.

                                                                                        {
#ifdef USE_GPU
#ifdef _DEBUG
   cudaMallocManaged((void **)&clover[0], 6*kvol*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&clover[1], 6*kvol*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)p, kfermHalo*sizeof(Complex), cudaMemAttachGlobal);
   cudaMallocManaged((void **)r, kferm*sizeof(Complex), cudaMemAttachGlobal);
   cudaMallocManaged((void **)x1, kfermHalo*sizeof(Complex), cudaMemAttachGlobal);
   cudaMallocManaged((void **)x2, kferm*sizeof(Complex), cudaMemAttachGlobal);
#else
   cudaMallocAsync((void **)&clover[0], 6*kvol*sizeof(Complex),streams[0]);
   cudaMallocAsync((void **)&clover[1], 6*kvol*sizeof(Complex),streams[1]);
   cudaMallocAsync((void **)p, kfermHalo*sizeof(Complex),streams[2]);
   cudaMallocAsync((void **)r, kferm*sizeof(Complex),streams[3]);
   cudaMallocAsync((void **)x1, kfermHalo*sizeof(Complex),streams[4]);
   cudaMallocAsync((void **)x2, kferm*sizeof(Complex),streams[5]);
#endif
   cudaDeviceSynchronise();
#else
   clover[0]=(Complex *)aligned_alloc(AVX,6*kvol*sizeof(Complex));
   clover[1]=(Complex *)aligned_alloc(AVX,6*kvol*sizeof(Complex));
   *p  = aligned_alloc(AVX, kfermHalo*sizeof(Complex));
   *r  = aligned_alloc(AVX, kferm*sizeof(Complex));
   *x1 = aligned_alloc(AVX, kfermHalo*sizeof(Complex));
   *x2 = aligned_alloc(AVX, kferm*sizeof(Complex));
#endif
}

References AVX, Complex, cudaDeviceSynchronise, kferm, kfermHalo, kvol, and streams.

Here is the caller graph for this function:

◆ P_allocate_f()

void P_allocate_f	(	Complex_f **	p_f,
		Complex_f **	r_f,
		Complex_f **	x1_f,
		Complex_f **	x2_f,
		Complex_f **	xi_f )

Allocates memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[out]	p_f	Holder for fermion field during inversion
[out]	r_f	Residue vector for conjugate gradient
[out]	x1_f,x2_f	\(M\) and \(M^\dagger M\)
[out]	xi_f	Accumulator array for conjugate gradient

Postcondition: Memory allocated

Definition at line 103 of file congrad.c.

                                                                                                        {
#ifdef USE_GPU
#ifdef _DEBUG
   cudaMallocManaged((void **)p_f, kfermHalo*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)r_f, kferm*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x1_f, kfermHalo*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x2_f, kferm*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)xi_f, kferm*sizeof(Complex_f),cudaMemAttachGlobal);
#else
   cudaMalloc((void **)p_f, kfermHalo*sizeof(Complex_f));
   cudaMalloc((void **)r_f, kferm*sizeof(Complex_f));
   cudaMalloc((void **)x1_f, kfermHalo*sizeof(Complex_f));
   cudaMalloc((void **)x2_f, kferm*sizeof(Complex_f));
   cudaMalloc((void **)xi_f, kferm*sizeof(Complex_f));
#endif
   cudaDeviceSynchronise();
#else
   *p_f  =  aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
   *r_f  =  aligned_alloc(AVX,kferm*sizeof(Complex_f));
   *x1_f =  aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
   *x2_f =  aligned_alloc(AVX,kferm*sizeof(Complex_f));
   *xi_f =  aligned_alloc(AVX,kferm*sizeof(Complex_f));
#endif
}

References AVX, Complex_f, cudaDeviceSynchronise, kferm, and kfermHalo.

Here is the caller graph for this function:

◆ P_free()

void P_free	(	Complex **	p,
		Complex **	r,
		Complex **	x1,
		Complex **	x2,
		Complex *	clover[2] )

Frees memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[in,out]	p	Holder for fermion field during inversion
[in,out]	x1,x2	\(M\) and \(M^\dagger M\)
[in,out]	r	Residue vector for conjugate gradient
[in,out]	clover	Clover fields

Postcondition: Memory freed

Definition at line 263 of file congrad.c.

                                                                                    {
#ifdef USE_GPU
#ifdef _DEBUG
   cudaFree(clover[0]); cudaFree(clover[1]);
   cudaFree(*p);cudaFree(*r);cudaFree(*x1); cudaFree(*x2);
#else
   cudaFreeAsync(clover[0],streams[0]); cudaFreeAsync(clover[1],streams[1]);
   cudaFreeAsync(*p,streams[2]);cudaFreeAsync(*r,streams[3]);cudaFreeAsync(*x1,streams[4]); cudaFreeAsync(*x2,streams[5]);
#endif
#else
   free(*p); free(*r); free(*x1); free(*x2);
   free(clover[0]); free(clover[1]);
#endif
}

References Complex, and streams.

Here is the caller graph for this function:

◆ P_free_f()

void P_free_f	(	Complex_f **	p_f,
		Complex_f **	r_f,
		Complex_f **	x1_f,
		Complex_f **	x2_f,
		Complex_f **	xi_f )

Frees memory needed for Congradp. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[in,out]	p_f	Holder for fermion field during inversion
[in,out]	r_f	Residue vector for conjugate gradient
[in,out]	x1_f,x2_f	\(M\) and \(M^\dagger M\)
[in,out]	xi_f	Conjugate gradient accumulator

Postcondition: Memory freed

Definition at line 243 of file congrad.c.

                                                                                                    {
#ifdef   USE_GPU
   cudaFree(*p_f); cudaFree(*r_f);cudaFree(*x1_f); cudaFree(*x2_f); cudaFree(*xi_f); 
#else
   free(*p_f); free(*r_f); free(*x1_f); free(*x2_f); free(*xi_f); 
#endif
}

References Complex_f.

Here is the caller graph for this function:

◆ Q_allocate()

void Q_allocate	(	Complex **	p,
		Complex **	x1,
		Complex **	x2,
		Complex *	clover[2] )

Allocates double precision memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[out]	p	Holder for fermion field during inversion
[out]	x1,x2	\(M\) and \(M^\dagger M\)
[out]	clover	Clover Fields

Postcondition: Memory allocated

Definition at line 64 of file congrad.c.

                                                                            {
   const char funcname[] = "Q_allocate";
#ifdef USE_GPU
#ifdef _DEBUG
   cudaMallocManaged((void **)&clover[0], 6*kvol*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&clover[1], 6*kvol*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)p, kferm2Halo*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x1, kferm2Halo*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x2, kferm2*sizeof(Complex),cudaMemAttachGlobal);
#else
   //First two have halo exchanges, so getting NCCL working is important
   cudaMallocAsync((void **)&clover[0], 6*kvol*sizeof(Complex),streams[0]);
   cudaMallocAsync((void **)&clover[1], 6*kvol*sizeof(Complex),streams[1]);
   cudaMallocAsync((void **)p, kferm2Halo*sizeof(Complex),streams[2]);
   cudaMallocAsync((void **)x1, kferm2Halo*sizeof(Complex),streams[3]);
   cudaMallocAsync((void **)x2, kferm2*sizeof(Complex),streams[4]);
#endif
#else
   clover[0]=(Complex *)aligned_alloc(AVX,6*kvol*sizeof(Complex));
   clover[1]=(Complex *)aligned_alloc(AVX,6*kvol*sizeof(Complex));
   *p=(Complex *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex));
   *x1=(Complex *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex));
   *x2=(Complex *)aligned_alloc(AVX,kferm2*sizeof(Complex));
#endif
   return;
}

References AVX, Complex, kferm2, kferm2Halo, kvol, and streams.

Here is the caller graph for this function:

◆ Q_allocate_f()

void Q_allocate_f	(	Complex_f **	p_f,
		Complex_f **	x1_f,
		Complex_f **	x2_f,
		Complex_f **	r_f,
		Complex_f **	X1_f )

Allocates memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[out]	p_f	Holder for fermion field during inversion
[out]	x1_f,x2_f	\(M\) and \(M^\dagger M\)
[out]	r_f	Residue vector for conjugate gradient
[out]	X1_f	Pseudofermion field

Postcondition: Memory allocated

Definition at line 23 of file congrad.c.

                                                                                                         {
   const char funcname[] = "Q_allocate";
#ifdef USE_GPU
#ifdef _DEBUG
   cudaMallocManaged((void **)p_f, kferm2Halo*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x1_f, kferm2Halo*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)x2_f, kferm2*sizeof(Complex_f),cudaMemAttachGlobal);
 
   cudaMallocManaged((void **)r_f, kferm2*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)X1_f, kferm2*sizeof(Complex_f),cudaMemAttachGlobal);
#else
   //First two have halo exchanges, so getting NCCL working is important
   cudaMallocAsync((void **)p_f, kferm2Halo*sizeof(Complex_f),streams[0]);
   cudaMallocAsync((void **)x1_f, kferm2Halo*sizeof(Complex_f),streams[1]);
   cudaMallocAsync((void **)x2_f, kferm2*sizeof(Complex_f),streams[2]);
 
   cudaMallocAsync((void **)r_f, kferm2*sizeof(Complex_f),streams[3]);
   cudaMallocAsync((void **)X1_f, kferm2*sizeof(Complex_f),streams[4]);
#endif
#else
   *p_f=(Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
   *x1_f=(Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
   *x2_f=(Complex_f *)aligned_alloc(AVX,kferm2*sizeof(Complex_f));
 
   *r_f=(Complex_f *)aligned_alloc(AVX,kferm2*sizeof(Complex_f));
   *X1_f=(Complex_f *)aligned_alloc(AVX,kferm2*sizeof(Complex_f));
#endif
   return;
}

References AVX, Complex_f, kferm2, kferm2Halo, and streams.

Here is the caller graph for this function:

◆ Q_free()

void Q_free	(	Complex **	p,
		Complex **	x1,
		Complex **	x2,
		Complex *	clover[2] )

Frees double precision memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[in,out]	p	Holder for fermion field during inversion
[in,out]	x1,x2	\(M\) and \(M^\dagger M\)
[in,out]	clover	clover fields

Postcondition: Memory freed

Definition at line 211 of file congrad.c.

                                                                        {
   const char funcname[] = "Qree";
#ifdef USE_GPU
#ifdef _DEBUG
   cudaDeviceSynchronise();
   cudaFree(*x1);cudaFree(*x2); cudaFree(*p);
   cudaFree(clover[0]); cudaFree(clover[1]);
#else
   //streams match the ones that allocated them.
   cudaFreeAsync(clover[0],streams[0]); cudaFreeAsync(clover[1],streams[1]);
   cudaFreeAsync(*p,streams[2]);cudaFreeAsync(*x1,streams[3]);cudaFreeAsync(*x2,streams[4]);
   cudaDeviceSynchronise();
#endif
#else
   free(*x1);free(*x2); free(*p); 
   free(clover[0]); free(clover[1]);
#endif
   return;
}

References Complex, cudaDeviceSynchronise, and streams.

Here is the caller graph for this function:

◆ Q_free_f()

void Q_free_f	(	Complex_f **	p_f,
		Complex_f **	x1_f,
		Complex_f **	x2_f,
		Complex_f **	r_f,
		Complex_f **	X1_f )

Frees memory needed for Congradq. Just to improve readability Note that since C does not modify it's arguments, you need to pass a pointer to the pointer you want to assign the memory to. This is similar behaviour to cudaMalloc.

Parameters

[in,out]	p_f	Holder for fermion field during inversion
[in,out]	x1_f,x2_f	\(M\) and \(M^\dagger M\)
[in,out]	r_f	Residue vector for conjugate gradient
[in,out]	X1_f	Pseudofermion field

Postcondition: Memory freed

Definition at line 181 of file congrad.c.

                                                                                                     {
   const char funcname[] = "Q_free";
#ifdef USE_GPU
#ifdef _DEBUG
   cudaDeviceSynchronise();
   cudaFree(*x1_f);cudaFree(*x2_f); cudaFree(*p_f);
   cudaFree(*r_f); cudaFree(*X1_f);
#else
   //streams match the ones that allocated them.
   cudaFreeAsync(*p_f,streams[0]);cudaFreeAsync(*x1_f,streams[1]);cudaFreeAsync(*x2_f,streams[2]);
   cudaDeviceSynchronise();
   cudaFreeAsync(*r_f,streams[3]); cudaFreeAsync(*X1_f,streams[4]);
#endif
#else
   free(*x1_f); free(*x2_f); free(*p_f); free(*r_f);free(*X1_f);
#endif
   return;
}

References Complex_f, cudaDeviceSynchronise, and streams.

Here is the caller graph for this function:

◆ Real_convert()

__global__ void Kernels::Real_convert	(	float *	a,
		double *	b,
		const unsigned int	len,
		const bool	dtof )

takes an array of real float and double precision numbers and converts the precision

Parameters

[in,out]	a	Float array
[in,out]	b	Double array
[in]	len	Number of elements to convert
[in]	dtof	If true, convert double to float. Otherwise convert float to double

Postcondition: Depending on the value of dtof, either the contents of a or b are overwritten with those of the other array in the opposite precision.

Definition at line 50 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;
 
   //True: Convert float to double
   if(dtof)
      for(unsigned int i = gthreadId; i<len;i+=gsize*bsize)
         a[i]=(float)b[i];
   //False: Convert double to float.
   else
      for(unsigned int i = gthreadId; i<len;i+=gsize*bsize)
         b[i]=(double)a[i];
}

Here is the caller graph for this function:

◆ reduce_sum()

template<typename T, unsigned int bsize>

__global__ void Kernels::reduce_sum	(	T *	g_in_data,
		T *	g_out_data,
		const unsigned int	n )

Performs a block reduction for sum.

Parameters

[in]	g_in_data	The input global data array
[out]	g_out_data	The output global data array
[in]	n	The size of the input array

Postcondition: sum saved in zeroth entry of g_out_data

Definition at line 543 of file cumatrices.cu.

                                                                                   {
         extern __shared__ T sdata[];  // stored in the shared memory
 
         // Each thread loading one element from global onto shared memory
         const unsigned short tid = threadIdx.x;
         unsigned int i = blockIdx.x*(bsize*2) + tid;
         const unsigned int gridSize = blockDim.x * 2 * gridDim.x;
         sdata[tid] = 0;
         while (i < n) {
            sdata[tid] += g_in_data[i];
            if (i + bsize < n) {
               sdata[tid] += g_in_data[i + bsize];
            }
            i += gridSize;
         }
         __syncthreads();
 
         // Perform reductions in steps, reducing thread synchronization
         // CUDA warps aren't guaranteed to be 32, and AMD vary. warpSize from the compiler fixes that.
         for(unsigned int s=bsize/2;s>=warpSize;s>>=1){
            if (tid < s)  
               sdata[tid] += sdata[tid + s];  __syncthreads();
         }
#if defined(__HIP_PLATFORM_AMD__) || defined(__HIP_PLATFORM_HCC__)
#define SU2_WARP_MASK 0xffffffffffffffffULL   // 64-bit lane mask (HIP/AMD)
#else
#define SU2_WARP_MASK 0xffffffffU             // 32-bit lane mask (CUDA)
#endif
         if (tid < warpSize)
            //Device::warpReduce_sum<T,bsize>(sdata, tid);
         {
            T val = sdata[tid];
            for (int offset = warpSize/2; offset > 0; offset >>= 1)
               val += __shfl_down_sync(SU2_WARP_MASK, val, offset);
            if (tid == 0) sdata[0] = val;
         }
 
         if (tid == 0){
            g_out_data[blockIdx.x] = sdata[0];
         }
      }

Here is the caller graph for this function:

◆ Reunitarise()

int Reunitarise ( Complex * ut[2] )

inline

Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.

If you're looking at the FORTRAN code be careful. There are two header files for the /trial/ header. One with u11 u12 (which was included here originally) and the other with u11t and u12t.

Parameters

[in,out] ut Trial fields to be reunitarised

Returns: Zero on success, integer error code otherwise

Postcondition: ut replaced with reunitarised gauge fields

Definition at line 342 of file su2hmc.c.

                                      {
   const char funcname[] = "Reunitarise";
#ifdef USE_GPU
   cuReunitarise(ut,dimGrid,dimBlock);
#else
#pragma omp parallel for simd
   for(unsigned short mu=0;mu<ndim;mu++)
      for(unsigned int i=0; i<kvol; i++){
         //Declaring anorm inside the loop will hopefully let the compiler know it
         //is safe to vectorise aggressively
         double anorm=sqrt(conj(ut[0][i+kvolHalo*mu])*ut[0][i+kvolHalo*mu]+conj(ut[1][i+kvolHalo*mu])*ut[1][i+kvolHalo*mu]);
         ut[0][i+kvolHalo*mu]/=anorm; ut[1][i+kvolHalo*mu]/=anorm;
      }
#endif
   return 0;
}

References Complex, conj(), cuReunitarise(), dimBlock, dimGrid, kvol, kvolHalo, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Transpose()

template<typename T>

__global__ void Kernels::Transpose	(	T *	out,
		const T *	in,
		const int	fast_in,
		const int	fast_out )

Swaps the order of the gauge field so that it is now SoA instead of AoS and it is nice and coalesced in memory.

Parameters

[out]	out	The flipped array
[in]	in	The original array
[in]	fast_out	The size of the slowest moving dimension. This is the lattice site when read in from disk
[in]	fast_in	The size of the fastest moving dimension. This is the direction index when read in from disk.

Definition at line 490 of file cumatrices.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;
 
         //The if/else here is only to ensure we maximise GPU bandwidth
         //Typically this is used to write back to the AoS/Coalseced format
         if(fast_out>fast_in){
            for(unsigned int x=gthreadId;x<fast_out;x+=gsize*bsize)
               for(unsigned int y=0; y<fast_in;y++)
                  out[y*fast_out+x]=in[x*fast_in+y];
         }
         //Typically this is used to write back to the SoA/saved config format
         else{
            for(unsigned int x=0; x<fast_out;x++)
               for(unsigned int y=gthreadId;y<fast_in;y+=gsize*bsize)
                  out[y*fast_out+x]=in[x*fast_in+y];
         }
      }

Here is the caller graph for this function:

◆ Transpose_c()

void Transpose_c	(	Complex_f *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 836 of file matrices.c.

                                                                              {
   const volatile char funcname[]="Transpose_c";
 
#ifdef USE_GPU
   cuTranspose_c(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   Complex_f *in = (Complex_f *)aligned_alloc(AVX,fast_in*fast_out*sizeof(Complex_f));
   memcpy(in,out,fast_in*fast_out*sizeof(Complex_f));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(int x=0;x<fast_out;x++)
         for(int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(int x=0; x<fast_out;x++)
         for(int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, Complex_f, cuTranspose_c(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ Transpose_d()

void Transpose_d	(	double *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 905 of file matrices.c.

                                                                           {
   const char funcname[]="Transpose_f";
 
#ifdef USE_GPU
   cuTranspose_d(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   double *in = (double *)aligned_alloc(AVX,fast_in*fast_out*sizeof(double));
   memcpy(in,out,fast_in*fast_out*sizeof(double));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(int x=0;x<fast_out;x++)
         for(int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(int x=0; x<fast_out;x++)
         for(int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, cuTranspose_d(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ Transpose_f()

void Transpose_f	(	float *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 882 of file matrices.c.

                                                                          {
   const char funcname[]="Transpose_f";
 
#ifdef USE_GPU
   cuTranspose_f(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   float *in = (float *)aligned_alloc(AVX,fast_in*fast_out*sizeof(float));
   memcpy(in,out,fast_in*fast_out*sizeof(float));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(int x=0;x<fast_out;x++)
         for(int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(int x=0; x<fast_out;x++)
         for(int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, cuTranspose_f(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ Transpose_I()

void Transpose_I	(	int *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 928 of file matrices.c.

                                                                        {
   const char funcname[]="Transpose_I";
 
#ifdef USE_GPU
   cuTranspose_I(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   int *in = (int *)aligned_alloc(AVX,fast_in*fast_out*sizeof(int));
   memcpy(in,out,fast_in*fast_out*sizeof(int));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(int x=0;x<fast_out;x++)
         for(int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(int x=0; x<fast_out;x++)
         for(int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, cuTranspose_I(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ Transpose_U()

void Transpose_U	(	unsigned int *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 951 of file matrices.c.

                                                                                 {
   const char funcname[]="Transpose_I";
 
#ifdef USE_GPU
   cuTranspose_U(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   unsigned int *in = (unsigned int *)aligned_alloc(AVX,fast_in*fast_out*sizeof(unsigned int));
   memcpy(in,out,fast_in*fast_out*sizeof(unsigned int));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(unsigned int x=0;x<fast_out;x++)
         for(unsigned int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(unsigned int x=0; x<fast_out;x++)
         for(unsigned int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, cuTranspose_U(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ Transpose_z()

void Transpose_z	(	Complex *	out,
		const int	fast_in,
		const int	fast_out )

inline

In place transpose used to convert from AoS to SoA memory layout.

Parameters

[in,out]	out	The array being transposed
[in]	fast_in	The old outermost/fastest index
[in]	fast_out	The new outermost/fastest index

Postcondition: Result overwrites existing out with transposed array

Definition at line 859 of file matrices.c.

                                                                            {
   const volatile char funcname[]="Transpose_c";
 
#ifdef USE_GPU
   cuTranspose_z(out,fast_in,fast_out,dimGrid,dimBlock);
#else
   Complex *in = (Complex *)aligned_alloc(AVX,fast_in*fast_out*sizeof(Complex));
   memcpy(in,out,fast_in*fast_out*sizeof(Complex));
   //Typically this is used to write back to the AoS/Coalseced format
   if(fast_out>fast_in){
      for(int x=0;x<fast_out;x++)
         for(int y=0; y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   //Typically this is used to write back to the SoA/saved config format
   else{
      for(int x=0; x<fast_out;x++)
         for(int y=0;y<fast_in;y++)
            out[y*fast_out+x]=in[x*fast_in+y];
   }
   free(in);
#endif
}

References AVX, Complex, cuTranspose_z(), dimBlock, and dimGrid.

Here is the call graph for this function:

◆ UpDownPart()

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

inline

Up/Down partitioning of the pseudofermion field.

Parameters

[in]	na	Flavour index
[out]	X0	Partitioned field
[in]	R1	Full pseudofermion field

Returns: Zero on success, integer error code otherwise

Postcondition: Result written into X0

Definition at line 327 of file su2hmc.c.

                                                                      {
#ifdef USE_GPU
   cuUpDownPart(na,X0,R1,dimBlock,dimGrid);
   cudaDeviceSynchronise();
#else
#pragma omp parallel for simd collapse(2) aligned(X0,R1:AVX)
   for(unsigned int i=0; i<kvol; i++)
      for(unsigned short idirac=0;idirac<ndirac;idirac++){
      //R1 has ngorkov spinors, but we only want the first four.
         X0[i+kvol*(0+nc*(idirac+ndirac*na))]=R1[i+kvol*(0+nc*idirac)];
         X0[i+kvol*(1+nc*(idirac+ndirac*na))]=R1[i+kvol*(1+nc*idirac)];
      }
#endif
   return 0;
}

References Complex, cudaDeviceSynchronise, cuUpDownPart(), dimBlock, dimGrid, kvol, nc, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ warpReduce_sum()

template<typename T, unsigned int bsize>

__device__ void Device::warpReduce_sum	(	volatile T *	sdata,
		const unsigned int	tid )

Performs a warp reduction for sum.

Parameters

[in,out]	sdata	The shared data array
[in]	tid	The thread ID

Postcondition: Sum written to zeroth index of sdata

Definition at line 21 of file cumatrices.cu.

                                                                               {
         if(bsize >= 64) sdata[tid] += sdata[tid + 32];
         if(bsize >= 32) sdata[tid] += sdata[tid + 16];
         if(bsize >= 16) sdata[tid] += sdata[tid + 8];
         if(bsize >= 8) sdata[tid] += sdata[tid + 4];
         if(bsize >= 4) sdata[tid] += sdata[tid + 2];
         if(bsize >= 2) sdata[tid] += sdata[tid + 1];
      }

◆ 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

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

Returns: Zero on success, integer error code otherwise

Postcondition: Contents of x replaced with output

Definition at line 301 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
#pragma omp parallel for simd aligned (x,y,table:AVX)
   for(unsigned int i=0; i<n; i++)
      x[i]=y[table[i+kvol*mu]+kvolHalo*mu];
   return 0;
}

References Complex, kvol, and kvolHalo.

su2hmc

Topics

Macros

Functions

Detailed Description

Macro Definition Documentation

◆ MAX

◆ MIN

Function Documentation

◆ blockInit()

◆ C_gather()

◆ ComplexConvert()

◆ conj()

◆ cuComplex_convert()

◆ cuFill_Small_Phi() [1/2]

◆ cuFill_Small_Phi() [2/2]

◆ cuMixed_Sumto()

◆ cuReal_convert()

◆ cureduce_sum_d()

◆ cureduce_sum_f()

◆ cuReunitarise() [1/2]

◆ cuReunitarise() [2/2]

◆ cuTranspose_c()

◆ cuTranspose_d()

◆ cuTranspose_f()

◆ cuTranspose_I()

◆ cuTranspose_U()

◆ cuTranspose_z()

◆ cuUpDownPart() [1/2]

◆ cuUpDownPart() [2/2]

◆ Fill_Small_Phi()

◆ Gather()

◆ Init()

◆ Init_CUDA()

◆ Mixed_Sumto()

◆ P_allocate()

◆ P_allocate_f()

◆ P_free()

◆ P_free_f()

◆ Q_allocate()

◆ Q_allocate_f()

◆ Q_free()

◆ Q_free_f()

◆ Real_convert()

◆ reduce_sum()

◆ Reunitarise()

◆ Transpose()

◆ Transpose_c()

◆ Transpose_d()

◆ Transpose_f()

◆ Transpose_I()

◆ Transpose_U()

◆ Transpose_z()

◆ UpDownPart()

◆ warpReduce_sum()

◆ Z_gather()