Clover Force

Collaboration diagram for Clover Force:

Classes
struct	Bilinear_a
	Structure of arrays for Hermitian bilinear \(X_{\mu\nu}\) in memory. More...
struct	Bilinear
	Hermitian bilinear \(X_{\mu\nu}\) on the local stack. More...

Functions
void	CalcXmunu (Bilinear_a Xmunu, Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, const unsigned short mu, const unsigned short nu)
	Gets \(X_{\mu\nu}\) for the clover force.
void	Clov_Force (double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, unsigned int *iu, unsigned int *id, const float akappa)
	Gets the clover contribution to the force.
void	cuCalcXmunu (Bilinear_a Xmunu, Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, const unsigned short mu, const unsigned short nu)
	CUDA wrapper for CalcXmunu. Only called during testing to be honest.
int	cuClov_Force (double *dSdpi, Complex_f *ut[nc], Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, const unsigned int *iu, const unsigned int *id, const float akappa)
	CUDA wrapper for Clover_Force.
static void	GetBilinear (Complex_f Z[nc *nc], Bilinear_a Xmn, unsigned int ind)
	Loads the compacted bilinear form into a \(2\times2\) complex valued matrix.
static void	GLeft (Complex_f out[4], const Complex_f G[2], const Complex_f X[4])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the left.
static void	GRight (Complex_f out[4], const Complex_f G[2], const Complex_f X[4])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the right.
static void	GSandwich (Complex_f out[4], Complex_f tmp[4], const Complex_f Gl[2], const Complex_f X[4], const Complex_f Gr[2])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the left and the right.
template<typename T>
__device__ void	Device::cuByGenLeft (T a[nc], const unsigned short gen)
	Multiply leaf (or part of one) by generator from left.
template<typename T>
__device__ void	Device::cuByGenRight (T a[nc], const unsigned short gen)
	Multiply leaf (or part of one) by generator from right.
template<typename T>
__device__ void	Device::cuGLeft (T out[4], const T G[2], const T X[4])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the left.
template<typename T>
__device__ void	Device::cuGRight (T out[4], const T G[2], const T X[4])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the right.
template<typename T>
__device__ void	Device::cuGSandwich (T out[4], T tmp[4], const T Gl[2], const T X[4], const T Gr[2])
	Multiplies \( X_{\mu\nu}\) by a gauge field from the left and the right.
__device__ void	Device::GetBilinear (Complex_f Z[nc *nc], Bilinear_a Xmn, unsigned int ind)
	Loads the compacted bilinear form into a \(2\times2\) complex valued matrix.
template<typename T>
__global__ void	Kernels::cuCalcXmunu (Bilinear_a Xmunu, const complex< T > *X1, const complex< T > *X2, const complex< T > *sigval, const unsigned short *sigin, const unsigned short clov)
	Gets \(X_{\mu\nu}\) for the clover force.
template<typename T>
__global__ void	Kernels::Clov_Force (double *dSdpi, const complex< T > *u11t, const complex< T > *u12t, Bilinear_a Xmn, const complex< T > *sigval, const unsigned short *sigin, const unsigned int *iu, const unsigned int *id, const float akappa, const unsigned short mu, const unsigned short nu)
	cuGets the clover contribution to the force

Detailed Description

Function Documentation

CalcXmunu()

void CalcXmunu	(	Bilinear_a	Xmunu,
		Complex_f *	X1,
		Complex_f *	X2,
		const Complex_f *	sigval,
		const unsigned short *	sigin,
		const unsigned short	mu,
		const unsigned short	nu )

Gets \(X_{\mu\nu}\) for the clover force.

Parameters

[out]	Xmunu	All Xmunu values
[in]	X1	Congrad output \(\left(M^\dagger M\right)\Phi\)
[in]	X2	\(M\left(M^\dagger M\right)^{-1}\Phi\)
[in]	sigval	\(\sigma_{\mu\nu}\) scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	Dirac index of \(\sigma_{\mu\nu}\)
[in]	mu,nu	Lattice directions

Postcondition

Bilinears written to Xmunu

Definition at line 437 of file clover.c.

                                                       {
   const char funcname[] = "Xmunu";
#ifdef USE_GPU
   cuCalcXmunu(Xmunu,X1,X2,sigval,sigin,mu,nu);
#else
   unsigned short clov;
   //Get sign and index of @f$\sigma_{\mu\nu}@f$ correct
   clov = (mu==0) ? nu-1 : mu+nu;
#pragma omp parallel for simd aligned(X1,X2:AVX)
   for(unsigned int i=0;i<kvol;i++){
      //Buffer. Four registers
      Bilinear Xmn; 
      Xmn.diag[0]=Xmn.diag[1]=0;Xmn.offd=0;
      for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
         const unsigned short sind = sigin[clov*ndirac+(idirac>>1)]<<1;
         const Complex_f sig = sigval[clov*ndirac+(idirac>>1)];
#pragma unroll
         for(unsigned short c1=0;c1<nc;c1++){
            //Spinors (rows) So we only load from memory once.
            const Complex_f X1s = X1[i+kvolHalo*(sind+c1)];
            const Complex_f X2s = X2[i+kvolHalo*(sind+c1)];
#pragma unroll
            for(unsigned short c2=0;c2<nc;c2++){
               //The second off diagonal term is the conjugate of the first
               if(c1==1&&c2==0)
                  continue;
               //Conjugated spinor (columns).
               const Complex_f X1c = conjf(X1[i+kvolHalo*(idirac+c2)]);
               const Complex_f X2c = conjf(X2[i+kvolHalo*(idirac+c2)]);
               if(c1==c2)
                  Xmn.diag[c1]+=creal(sig*(X2s*X1c+X1s*X2c));
               else
                  Xmn.offd+=sig*(X2s*X1c+X1s*X2c);
 
            }
         }
      }
      //And write back to global memory.
#pragma unroll
      //Write the diagonals
      for(unsigned short c=0;c<nc;c++)
         Xmunu.diag[i+kvolHalo*c]=Xmn.diag[c];
      //And the off diagonal terms
      Xmunu.offd[i]=Xmn.offd;
   }
#endif
   return;
}

References Complex_f, creal, cuCalcXmunu(), Bilinear::diag, Bilinear_a::diag, kvol, kvolHalo, nc, ndirac, Bilinear::offd, and Bilinear_a::offd.

Here is the call graph for this function:

Here is the caller graph for this function:

Clov_Force() [1/2]

void Clov_Force	(	double *	dSdpi,
		Complex_f *	ut[2],
		Complex_f *	X1,
		Complex_f *	X2,
		const Complex_f *	sigval,
		const unsigned short *	sigin,
		unsigned int *	iu,
		unsigned int *	id,
		const float	akappa )

Gets the clover contribution to the force.

Parameters

[in,out]	dSdpi	Force
[in]	ut	Gauge fields
[in]	X1	Congrad output \(\left(M^\dagger M\right)\Phi\)
[in]	X2	\(M\left(M^\dagger M\right)^{-1}\Phi\)
[in]	sigval	\(\sigma_{\mu\nu}\) scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	Dirac index of \(\sigma_{\mu\nu}\)
[in]	iu,id	Neighbouring sites
[in]	akappa	Hopping parameter

Postcondition

Force contribution added to dSdpi

Definition at line 534 of file clover.c.

                                                                                          {
   const char funcname[] = "Clov_Force";
#ifdef USE_GPU
   cuClov_Force(dSdpi,ut,X1,X2,sigval,sigin,iu,id,akappa);
#else
   //Allocate the @f$X_{\mu\nu}@f$ array
   unsigned short nclov=6; unsigned short clov=0;
   Bilinear_a Xmn[nclov];
   //And get the @f$X_{\mu\nu}@f$ values
   //Loop over @f$\mu@f$ and @f$\nu@f$. Symmetry means we actually only need half the terms
   for(unsigned short mu=0;mu<ndim-1;mu++)
      for(unsigned short nu=mu+1;nu<ndim;nu++)
         if(mu!=nu){
            clov = (mu==0) ? nu-1 : mu+nu;
            Xmn[clov].diag=(float *)aligned_alloc(AVX,2*kvolHalo*sizeof(float));
            Xmn[clov].offd=(Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f));
            CalcXmunu(Xmn[clov],X1,X2,sigval,sigin,mu,nu);
#if(nproc>1)
            SHalo_swap_all(Xmn[clov].diag,2);
            CHalo_swap_all(Xmn[clov].offd,1);
#endif
         }
   for(unsigned short mu=0;mu<ndim;mu++)
      for(unsigned short nu=0;nu<ndim;nu++)
         if(mu!=nu){
            if(mu<nu)
               clov = (mu==0) ? nu-1 : mu+nu;
            else
               clov = (nu==0) ? mu-1 : mu+nu;
#pragma omp parallel for
            for(unsigned int i=0;i<kvol;i++){
               //This is where it gets messy. Using HiRep/OpenQCD labelling for different intermediate values
               //But recycling to reduce register pressure on GPU
               //First up, W0, W1 and W6 match their Documentation values
               Complex_f W0[2], W1[2], W6[2];   
               //Get the correct site. 
               unsigned int ind = id[i+kvol*nu];
               //Gauge field @f$U_\nu\left(i-\hat{\nu}\right)
               W1[0]=ut[0][ind+kvolHalo*nu]; W1[1]=ut[1][ind+kvolHalo*nu];
 
               //@f$Z_2=X_{\mu\nu}\left(x-\hat{\nu}\right)@f$
               Complex_f Z[nc*nc];
               GetBilinear(Z,Xmn[clov],ind);
 
               //W0 is @f$U^\dagger_\mu\left(x-\hat{nu}\right)@f$
               W0[0]=conjf(ut[0][ind+kvolHalo*mu]); W0[1]=-ut[1][ind+kvolHalo*mu];
 
               //Need a temporary Z buffers for the intermediate result
               Complex_f Zbuff1[nc*nc], Zbuff2[nc*nc];
               GSandwich(Zbuff1,Zbuff2,W0,Z,W1);
 
               //@f$W_6=W_0 W_1@f$
               W6[0]=W0[0]*W1[0]-W0[1]*conjf(W1[1]); W6[1]=W0[0]*W1[1]+W0[1]*conjf(W1[0]);
 
               //Z_3 is the @f$X_{\mu\nu}\left(x+\hat{\mu}-\hat{\nu}\right)@f$. Store in Z
               ind=iu[ind+kvol*mu];
               GetBilinear(Z,Xmn[clov],ind);
 
               //Need a second Zbuffer for another intermediate result.
               GRight(Zbuff2,W6,Z);
               //Sum the two results into Zbuff1. Then scale by -W5
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++)
                  Zbuff1[c]+=Zbuff2[c];
               //W5 is @f$U^\dagger_\nu\left(x+\hat{\mu}-\hat{\nu}\right)@f$
               Complex_f W5[2];
               W5[0]=conjf(ut[0][ind+kvolHalo*nu]); W5[1]=-ut[1][ind+kvolHalo*nu];
               //Now multiply by @f$W_5@f$ from the left into Zbuff2
               GLeft(Zbuff2,W5,Zbuff1);
 
               //Intermediate results from the four parts of the sum.
               Complex_f F_int[4];
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++)
                  //Negative as it is @f$-W_5@f$
                  F_int[c]=-Zbuff2[c];
 
               //Now we repeat for the last term in the sum. Recycling along the way.
               //First store @f$W_2=U_\nu\left(x+\hat{\mu}\right)@f$ into W0.
               ind=iu[i+kvol*mu];
               W0[0]=ut[0][ind+kvolHalo*nu]; W0[1]=ut[1][ind+kvolHalo*nu];
               //@f$W_3=U^\dagger_\mu\left(x+\hat{\nu}\right)@f$. Storing it in W1
               ind=iu[i+kvol*nu];
               W1[0]=conjf(ut[0][ind+kvolHalo*mu]); W1[1]=-ut[1][ind+kvolHalo*mu];
               //@f$Z_4=X_{\mu\nu}\left(x+\hat{\mu}+\hat{\nu}\right)@f$. Storing in Z
               ind=iu[ind+kvol*mu];
               GetBilinear(Z,Xmn[clov],ind);
               //Calculate and write into Zbuff1
               GSandwich(Zbuff1,Zbuff2,W0,Z,W1);
 
               //@f$W_7=W_0 W_1@f$
               Complex_f W7[2];
               W7[0]=W0[0]*W1[0]-W0[1]*conjf(W1[1]); W7[1]=W0[0]*W1[1]+W0[1]*conjf(W1[0]);
               //@f$Z_5=X_{\mu\nu}\left(x+\hat{\nu}\right)@f$
               ind=iu[i+kvol*nu]; 
               GetBilinear(Z,Xmn[clov],ind);
               //And calculate the second term
               GLeft(Zbuff2,W7,Z);
               //Sum the two results into Zbuff1.
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++)
                  Zbuff1[c]+=Zbuff2[c];
 
               //W4 is @f$U^\dagger_\nu\left(x\right)@f$
               Complex_f W4[2];
               W4[0]=conjf(ut[0][i+kvolHalo*nu]); W4[1]=-ut[1][i+kvolHalo*nu];
               //Now multiply by @f$W_4@f$ from the right into Zbuff2
               GRight(Zbuff2,W4,Zbuff1);
 
               //Intermediate results from the four parts of the sum.
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++)
                  F_int[c]+=Zbuff2[c];
               //The last thing we need is @f$W_8=W_7W_4-W_5W_6@f$. Do it in parts and store intermediates in W0 and W1
               W0[0]=W7[0]*W4[0]-W7[1]*conjf(W4[1]); W0[1]=W7[0]*W4[1]+W7[1]*conjf(W4[0]);
               W1[0]=W5[0]*W6[0]-W5[1]*conjf(W6[1]); W1[1]=W5[0]*W6[1]+W5[1]*conjf(W6[0]);
               //Store W8 in W0
               W0[0]-=W1[0]; W0[1]-=W1[1];
 
               //Now load @f$Z_0=X_{\mu\nu}(x)@f$
               GetBilinear(Z,Xmn[clov],i);
               GLeft(Zbuff1,W0,Z);
               //And sum intermediate
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++)
                  F_int[c]+=Zbuff1[c];
 
               //Now load @f$Z_1=X_{\mu\nu}\left(x+\hat{mu})@f$
               ind=iu[i+kvol*mu];
               GetBilinear(Z,Xmn[clov],ind);
               GRight(Zbuff1,W0,Z);
               //And sum intermediate
#pragma unroll
               for(unsigned short c=0;c<nc*nc;c++){
                  F_int[c]+=Zbuff1[c];
                  //See if this works...
                  F_int[c]*=-I;
               }
 
               //Excellent. Now we just need to multiply by the derivative term
               W0[0]=ut[0][i+kvolHalo*mu]; W0[1]=ut[1][i+kvolHalo*mu];
               for(unsigned short gen=0;gen<nadj;gen++){
                  W1[0]=W0[0]; W1[1]=W0[1];
                  ByGenLeft(W1,gen);
                  GLeft(Zbuff1,W1,F_int);
                  //Sum of the real part of the trace.
                  float dSdpis=crealf(Zbuff1[0])+crealf(Zbuff1[3]);
                  if(mu<nu)
                     dSdpi[i+kvol*(gen*ndim+mu)] -=akappa*dSdpis/4.0f;
                  else
                     dSdpi[i+kvol*(gen*ndim+mu)] +=akappa*dSdpis/4.0f;
               }
            }
         }
   for(clov=0;clov<nclov;clov++){
      free(Xmn[clov].diag); free(Xmn[clov].offd);
   }
#endif
   return;
}

References AVX, ByGenLeft(), CalcXmunu(), CHalo_swap_all(), Complex_f, cuClov_Force(), Bilinear_a::diag, GetBilinear(), GLeft(), GRight(), GSandwich(), I, kvol, kvolHalo, nadj, nc, ndim, Bilinear_a::offd, and SHalo_swap_all().

Here is the call graph for this function:

Here is the caller graph for this function:

Clov_Force() [2/2]

template<typename T>

__global__ void Kernels::Clov_Force	(	double *	dSdpi,
		const complex< T > *	u11t,
		const complex< T > *	u12t,
		Bilinear_a	Xmn,
		const complex< T > *	sigval,
		const unsigned short *	sigin,
		const unsigned int *	iu,
		const unsigned int *	id,
		const float	akappa,
		const unsigned short	mu,
		const unsigned short	nu )

cuGets the clover contribution to the force

Parameters

[out]	dSdpi	Force.
[in]	u11t,u12t	Gauge fields
[in]	Xmn	Pseudofermion bilinear contraction
[in]	sigval	\(\sigma_{\mu\nu}\) scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	Dirac index of \(\sigma_{\mu\nu}\)
[in]	iu,id	Neighbouring sites
[in]	akappa	Hopping Parameter
[in]	mu,nu	Clover direction

Postcondition

Clover force is added to dSdpi

Definition at line 427 of file cuclover.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;
         //And get the @f$X_{\mu\nu}@f$ values
         //Loop over @f$\mu@f$ and @f$\nu@f$,
         for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
            //Buffer for intermediate force calculation. One for each generator.
            //This is where it gets messy. Using HiRep/OpenQCD labelling for different intermediate values
            //But recycling to reduce register pressure on GPU
            //First up, W0, W1 and W6 match their Documentation values
            complex<T> W0[2]; complex<T> W1[2]; complex<T> W6[2]; 
            //Get the correct site. Originally ind and  stood for up and down. Then I realised only one was needed
            //at a time and am too lazy to change it everywhere.
            unsigned int ind = id[i+kvol*nu];
            //Gauge field @f$U_\nu\left(i-\hat{\nu}\right)
            W1[0]=u11t[ind+kvolHalo*nu]; W1[1]=u12t[ind+kvolHalo*nu];
 
            //@f$Z_2=X_{\mu\nu}\left(i-\hat{\nu}\right)@f$
            complex<T> Z[nc*nc];
            Device::GetBilinear(Z,Xmn,ind);
 
            //W0 is @f$U^\dagger_\mu@f(x-\hat{nu}\right)@f$
            W0[0]=conj(u11t[ind+kvolHalo*mu]); W0[1]=-u12t[ind+kvolHalo*mu];
 
            //Need a temporary Z buffers for the intermediate result
            complex<T> Zbuff1[nc*nc];complex<T> Zbuff2[nc*nc];
            Device::cuGSandwich(Zbuff1,Zbuff2,W0,Z,W1);
 
            //@f$W_6=W_0 W_1@f$
            W6[0]=W0[0]*W1[0]-W0[1]*conj(W1[1]); W6[1]=W0[0]*W1[1]+W0[1]*conj(W1[0]);
 
            //Z3 is the @f$X_{\mu\nu}\left(x+\hat{\mu}-\hat{\nu}\right)@f$. Store in Z
            ind=iu[ind+kvol*mu];
            Device::GetBilinear(Z,Xmn,ind);
 
            //Need a second Zbuffer for another intermediate result.
            Device::cuGRight(Zbuff2,W6,Z);
            //Sum the two results into Zbuff1. Then scale by -W5
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++)
               Zbuff1[c]+=Zbuff2[c];
            //W5 is @f$U^\dagger_\nu\left(x+\hat{\mu}-\hat{\nu}\right)@f$
            complex<T> W5[2];
            W5[0]=conj(u11t[ind+kvolHalo*nu]); W5[1]=-u12t[ind+kvolHalo*nu];
            //Now multiply by @f$W_5@f$ from the left into Zbuff2
            Device::cuGLeft(Zbuff2,W5,Zbuff1);
 
            //Intermediate results from the four parts of the sum.
            complex<T> F_int[4];
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++)
               //Negative as it is @f$-W_5@f$
               F_int[c]=-Zbuff2[c];
 
            //Now we repeat for the last term in the sum. Recycling along the way.
            //First store @f$W_2=U_\nu\left(x+\hat{\mu}\right)@f$ into W0.
            ind=iu[i+kvol*mu];
            W0[0]=u11t[ind+kvolHalo*nu]; W0[1]=u12t[ind+kvolHalo*nu];
            //@f$W_3=U^\dagger_\mu\left(x+\hat{\nu}\right). Storing it in W1
            ind=iu[i+kvol*nu];
            W1[0]=conj(u11t[ind+kvolHalo*mu]); W1[1]=-u12t[ind+kvolHalo*mu];
            //@f$Z_4=X_{\mu\nu}\left(x+\hat{\mu}+\hat{\nu}\right)@f$. Storing in Z
            ind=iu[ind+kvol*mu];
            Device::GetBilinear(Z,Xmn,ind);
            //Calculate and write into Zbuff1
            Device::cuGSandwich(Zbuff1,Zbuff2,W0,Z,W1);
 
            //@f$W_7=W_0 W_1@f$
            complex<T> W7[2];
            W7[0]=W0[0]*W1[0]-W0[1]*conj(W1[1]); W7[1]=W0[0]*W1[1]+W0[1]*conj(W1[0]);
            //@f$Z_5=X_{\mu\nu}\left(x+\hat{\nu}\right)@f$
            ind=iu[i+kvol*nu]; 
            Device::GetBilinear(Z,Xmn,ind);
            //And calculate the second term
            Device::cuGLeft(Zbuff2,W7,Z);
            //Sum the two results into Zbuff1.
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++)
               Zbuff1[c]+=Zbuff2[c];
            //W4 is @f$U^\dagger_\nu\left(x\right)@f$
            complex<T> W4[2];
            W4[0]=conj(u11t[i+kvolHalo*nu]); W4[1]=-u12t[i+kvolHalo*nu];
            //Now multiply by @f$W_4@f$ from the right into Zbuff2
            Device::cuGRight(Zbuff2,W4,Zbuff1);
 
            //Intermediate results from the four parts of the sum.
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++)
               F_int[c]+=Zbuff2[c];
            //The last thing we need is @f$W_8=W_7W_4-W_5W_6@f$. Do it in parts and store intermediates in W0 and W1
            W0[0]=W7[0]*W4[0]-W7[1]*conj(W4[1]); W0[1]=W7[0]*W4[1]+W7[1]*conj(W4[0]);
            W1[0]=W5[0]*W6[0]-W5[1]*conj(W6[1]); W1[1]=W5[0]*W6[1]+W5[1]*conj(W6[0]);
            //Store W8 in W0
            W0[0]-=W1[0]; W0[1]-=W1[1];
 
            //Now load @f$@Z_0=X_{\mu\nu}(x)@f$
            Device::GetBilinear(Z,Xmn,i);
            Device::cuGLeft(Zbuff1,W0,Z);
            //And sum intermediate
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++)
               F_int[c]+=Zbuff1[c];
 
            //Now load @f$@Z_1=X_{\mu\nu}(x)@f$
            ind=iu[i+kvol*mu];
            Device::GetBilinear(Z,Xmn,ind);
            Device::cuGRight(Zbuff1,W0,Z);
            //And sum intermediate
#pragma unroll
            for(unsigned short c=0;c<nc*nc;c++){
               F_int[c]+=Zbuff1[c];
               F_int[c]*=-I_f;
            }
 
            //Excellent. Now we just need to multiply by the derivative term
            W0[0]=u11t[i+kvolHalo*mu]; W0[1]=u12t[i+kvolHalo*mu];
            for(unsigned short gen=0;gen<nadj;gen++){
               W1[0]=W0[0]; W1[1]=W0[1];
               cuByGenLeft(W1,gen);
               Device::cuGLeft(Zbuff1,W1,F_int);
               //Sum of the real part of the trace.
               float dSdpis=creal(Zbuff1[0])+creal(Zbuff1[3]);
               //tmp lets us control the number of registers explictly
               if(mu<nu)
                  dSdpi[i+kvol*(gen*ndim+mu)] -=akappa*dSdpis/4.0f;
               else
                  dSdpi[i+kvol*(gen*ndim+mu)] +=akappa*dSdpis/4.0f;
 
            }
         }
         return;
      }

References conj(), creal, Device::cuByGenLeft(), Device::cuGLeft(), Device::cuGRight(), Device::cuGSandwich(), Device::GetBilinear(), I_f, kvol, kvolHalo, nadj, nc, and ndim.

Here is the call graph for this function:

Here is the caller graph for this function:

cuByGenLeft()

template<typename T>

__device__ void Device::cuByGenLeft	(	T	a[nc],
		const unsigned short	gen )

Multiply leaf (or part of one) by generator from left.

The leaves contributing to each force term need to be scaled by the generator, but the generator appears at different points in each leaf. This routine multiples by the generator from the left side.

Parameters

[in,out]	a	The leaf or partial leaf. Result saved in place
[in]	gen	What generator are we multiplying by?

\(i\sigma_x\)

\(i\sigma_y\)

\(i\sigma_z\)

Definition at line 24 of file cuclover.cu.

                                                                   {
         T tmp = a[0];
         switch(gen){
            case(0):
               a[0] = T(-a[1].imag(), -a[1].real());
               a[1] = T( tmp.imag(),  tmp.real());
               break;
            case(1):
               a[0] = -conj(a[1]);
               a[1] = conj(tmp);
               break;
            case(2):
               a[0] = T(-a[0].imag(), a[0].real());
               a[1] = T(-a[1].imag(), a[1].real());
               break;
         }
         return;
      }

References conj(), and nc.

Here is the call graph for this function:

Here is the caller graph for this function:

cuByGenRight()

template<typename T>

__device__ void Device::cuByGenRight	(	T	a[nc],
		const unsigned short	gen )

Multiply leaf (or part of one) by generator from right.

The leaves contributing to each force term need to be scaled by the generator, but the generator appears at different points in each leaf. This routine multiples by the generator from the right side.

Parameters

[in,out]	a	The leaf or partial leaf. Result saved in place
[in]	gen	What generator are we multiplying by?

\(i\sigma_x\)

\(i\sigma_y\)

\(i\sigma_z\)

Definition at line 56 of file cuclover.cu.

                                                                    {
         T tmp = a[0];
         switch(gen){
            case(0):
               a[0] = T(-a[1].imag(), a[1].real());
               a[1] = T(-tmp.imag(),  tmp.real());
               break;
            case(1):
               a[0]=-a[1]; a[1]=tmp;
               break;
            case(2):
               a[0] = T(-a[0].imag(),  a[0].real());
               a[1] = T( a[1].imag(), -a[1].real());
               break;
         }
         return;
      }

References nc.

cuCalcXmunu() [1/2]

void cuCalcXmunu	(	Bilinear_a	Xmunu,
		Complex_f *	X1,
		Complex_f *	X2,
		const Complex_f *	sigval,
		const unsigned short *	sigin,
		const unsigned short	mu,
		const unsigned short	nu )

CUDA wrapper for CalcXmunu. Only called during testing to be honest.

Parameters

[out]	Xmunu	All Xmunu values
[in]	X1	Congrad output \(\left(M^\dagger M\right)\Phi\)
[in]	X2	\(M\left(M^\dagger M\right)^{-1}\Phi\)
[in]	sigval	\(\sigma_{\mu\nu}\) scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	Dirac index of \(\sigma_{\mu\nu}\)
[in]	mu,nu	Lattice directions

Postcondition

Bilinears written to Xmunu

Definition at line 725 of file cuclover.cu.

                                                       {
   //Get sign and index of @f$\sigma_{\mu\nu}@f correct
   unsigned short clov = (mu==0) ? nu-1 : mu+nu;
   Kernels::cuCalcXmunu<<<dimGrid,dimBlock,0,streams[clov]>>>(Xmunu,X1,X2,sigval,sigin,clov);
   return;
}

References Complex_f, and Kernels::cuCalcXmunu().

Here is the call graph for this function:

Here is the caller graph for this function:

cuCalcXmunu() [2/2]

template<typename T>

__global__ void Kernels::cuCalcXmunu	(	Bilinear_a	Xmunu,
		const complex< T > *	X1,
		const complex< T > *	X2,
		const complex< T > *	sigval,
		const unsigned short *	sigin,
		const unsigned short	clov )

Gets \(X_{\mu\nu}\) for the clover force.

Parameters

[out]	Xmunu	All Xmunu values
[in]	X1	Congrad output \(\left(M^\dagger M\right)\Phi\)
[in]	X2	\(M\left(M^\dagger M\right)^{-1}\Phi\)
[in]	sigval	\(\sigma_{\mu\nu}\) scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	Dirac index of \(\sigma_{\mu\nu}\)
[in]	clov	Index of clover being used

Postcondition

Contents of Xmunu overwritten

Definition at line 369 of file cuclover.cu.

                                                                                            {
         const char funcname[] = "Xmunu";
         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;
         //Get sign and index of @f$\sigma_{\mu\nu}@f correct
         for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
            //Buffer. Four registers
            Bilinear Xmn; 
            Xmn.diag[0]=Xmn.diag[1]=0;Xmn.offd=0;
            for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
               const unsigned short sind = sigin[clov*ndirac+(idirac>>1)]<<1;
               const complex<T> sig = sigval[clov*ndirac+(idirac>>1)];
               for(unsigned short c1=0;c1<nc;c1++){
                  //Spinors (rows) So we only load from memory once.
                  const complex<T> X1s = X1[i+kvolHalo*(sind+c1)];
                  const complex<T> X2s = X2[i+kvolHalo*(sind+c1)];
                  for(unsigned short c2=0;c2<nc;c2++){
                     //The second off diagonal term is the conjugate of the first
                     if(c1==1&&c2==0)
                        continue;
                     //Conjugated spinor (columns).
                     const complex<T> X1c = conj(X1[i+kvolHalo*(idirac+c2)]);
                     const complex<T> X2c = conj(X2[i+kvolHalo*(idirac+c2)]);
                     if(c1==c2)
                        Xmn.diag[c1]+=(T)(sig*(X2s*X1c+X1s*X2c)).real();
                     else
                        Xmn.offd+=sig*(X2s*X1c+X1s*X2c);
                  }
               }
            }
            //Write the diagonals
            for(unsigned short c=0;c<nc;c++)
               Xmunu.diag[i+kvolHalo*c]=Xmn.diag[c];
            //And the off diagonal terms
            Xmunu.offd[i]=Xmn.offd;
         }
         return;
      }

References conj(), Bilinear::diag, Bilinear_a::diag, kvol, kvolHalo, nc, ndirac, Bilinear::offd, and Bilinear_a::offd.

Here is the call graph for this function:

Here is the caller graph for this function:

cuClov_Force()

int cuClov_Force	(	double *	dSdpi,
		Complex_f *	ut[nc],
		Complex_f *	X1,
		Complex_f *	X2,
		const Complex_f *	sigval,
		const unsigned short *	sigin,
		const unsigned int *	iu,
		const unsigned int *	id,
		const float	akappa )

CUDA wrapper for Clover_Force.

Parameters

[in,out]	dSdpi	Force
[in]	ut	Gauge fields
[in]	X1	\(\left(M^\dagger M\right)^{-1} \Psi\)
[in]	X2	\(M\left(M^\dagger M\right)^{-1} \Psi\)
[in]	sigval	\( \sigma_{\mu\nu}\) entries scaled by \(c_sw\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	iu,id	Up/down indices
[in]	akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Postcondition

Force contribution added to dSdpi

Definition at line 732 of file cuclover.cu.

                                                                                                      {
   const char funcname[]="Clov_Force";
 
   const unsigned short nclov = (ndim-1)*(ndim-2);
   //Too many pointers here but not bothered doing it correctly. Overhead is basically zero.
   Bilinear_a Xmn[nclov];
   //Allocate half-leaf memory. We will have one stream for each direction
   for(unsigned short mu=0;mu<ndim-1;mu++)
      for(unsigned short nu=mu;nu<ndim;nu++)
         if(mu!=nu){
            //Get sign and index of @f$\sigma_{\mu\nu}@f correct
            unsigned short clov = (mu==0) ? nu-1 : mu+nu;
            //Allocate and evaluate @f$X_{\mu\nu}@f$ terms
            cudaMallocAsync((void **)&Xmn[clov].diag,nc*kvol*sizeof(float),streams[clov]);
            cudaMallocAsync((void **)&Xmn[clov].offd,nc*kvol*sizeof(complex<float>),streams[clov]);
            Kernels::cuCalcXmunu<<<dimGrid,dimBlock,0,streams[clov]>>>(Xmn[clov],X1,X2,sigval,sigin,clov);
         }
   cudaDeviceSynchronise();
 
   for(unsigned short mu=0;mu<ndim;mu++)
      for(unsigned short nu=0;nu<ndim;nu++)
         if(mu!=nu){
            unsigned short clov=0;
            if(mu<nu)
               clov = (mu==0) ? nu-1 : mu+nu;
            else
               clov = (nu==0) ? mu-1 : mu+nu;
            //Compute force for @f$\mu\nu@f$ and @f$\nu\mu@f$
            Kernels::Clov_Force<<<dimGrid,dimBlock,0,streams[mu]>>>(dSdpi,ut[0],ut[1],Xmn[clov],\
                  sigval,sigin,iu,id,akappa,mu,nu);
 
         }
   cudaDeviceSynchronise();
   for(unsigned short clov=0;clov<nclov;clov++){
      cudaFreeAsync(Xmn[clov].diag,NULL); cudaFreeAsync(Xmn[clov].offd,NULL);
   }
   cudaDeviceSynchronise();
   return 0;
}

References Kernels::Clov_Force(), Complex_f, Kernels::cuCalcXmunu(), cudaDeviceSynchronise, kvol, nc, ndim, and streams.

Here is the call graph for this function:

Here is the caller graph for this function:

cuGLeft()

template<typename T>

__device__ void Device::cuGLeft	(	T	out[4],
		const T	G[2],
		const T	X[4] )

Multiplies \( X_{\mu\nu}\) by a gauge field from the left.

Parameters

[out]	out	Result
[in]	X	\(X_{\mu\nu}(x)\)
[in]	G	Gauge field

Definition at line 224 of file cuclover.cu.

                                                                  {
         out[0]=G[0]*X[0]+G[1]*X[2];
         out[1]=G[0]*X[1]+G[1]*X[3];
         out[2]=-conj(G[1])*X[0]+conj(G[0])*X[2];
         out[3]=-conj(G[1])*X[1]+conj(G[0])*X[3];
         return;
      }

References conj().

Here is the call graph for this function:

Here is the caller graph for this function:

cuGRight()

template<typename T>

__device__ void Device::cuGRight	(	T	out[4],
		const T	G[2],
		const T	X[4] )

Multiplies \( X_{\mu\nu}\) by a gauge field from the right.

Parameters

[out]	out	Result
[in]	X	\(X_{\mu\nu}(x)\)
[in]	G	Gauge field

Definition at line 240 of file cuclover.cu.

                                                                   {
         out[0]=G[0]*X[0]-conj(G[1])*X[1];
         out[1]=G[1]*X[0]+conj(G[0])*X[1];
         out[2]=G[0]*X[2]-conj(G[1])*X[3];
         out[3]=G[1]*X[2]+conj(G[0])*X[3];
         return;
      }

References conj().

Here is the call graph for this function:

Here is the caller graph for this function:

cuGSandwich()

template<typename T>

__device__ void Device::cuGSandwich	(	T	out[4],
		T	tmp[4],
		const T	Gl[2],
		const T	X[4],
		const T	Gr[2] )

Multiplies \( X_{\mu\nu}\) by a gauge field from the left and the right.

Parameters

[out]	out	Result
[in]	tmp	Buffer for intermediate result. Passing as an argument to reduce register pressure.
[in]	X	\(X_{\mu\nu}(x)\)
[in]	Gl,Gr	Left/Right Gauge fields

Definition at line 257 of file cuclover.cu.

                                                                                               {
         cuGRight(tmp,Gr,X);
         cuGLeft(out,Gl,tmp);
         return;
      }

References cuGLeft(), and cuGRight().

Here is the call graph for this function:

Here is the caller graph for this function:

GetBilinear() [1/2]

__device__ void Device::GetBilinear	(	Complex_f	Z[nc *nc],
		Bilinear_a	Xmn,
		unsigned int	ind )

Loads the compacted bilinear form into a \(2\times2\) complex valued matrix.

Parameters

[out]	Z	Complex matrix on stack
[in]	Xmn	Hermitian form of \(X_{\mu\nu}\) in memory
[in]	ind	Index of \(X_{\mu\nu}\) we're interested in

Postcondition

Contents of Z are replaced by bilinear from memory

Definition at line 272 of file cuclover.cu.

                                                                                   {
      Z[0]=Xmn.diag[ind]; Z[3]=Xmn.diag[ind+kvolHalo];
      Z[1]=Xmn.offd[ind]; Z[2]=conj(Z[1]);
      return;
   }

References Complex_f, conj(), Bilinear_a::diag, kvolHalo, nc, and Bilinear_a::offd.

Here is the call graph for this function:

Here is the caller graph for this function:

GetBilinear() [2/2]

void GetBilinear ( Complex_f Z[nc *nc], Bilinear_a Xmn, unsigned int ind )

inlinestatic

Loads the compacted bilinear form into a \(2\times2\) complex valued matrix.

Parameters

[out]	Z	Complex matrix on stack
[in]	Xmn	Hermitian form of \(X_{\mu\nu}\) in memory
[in]	ind	Index of \(X_{\mu\nu}\) we're interested in

Postcondition

Contents of Z are replaced by bilinear from memory

Definition at line 432 of file clover.c.

                                                                                   {
   Z[0]=Xmn.diag[ind]; Z[3]=Xmn.diag[ind+kvolHalo];
   Z[1]=Xmn.offd[ind]; Z[2]=conjf(Z[1]);
   return;
}

References Complex_f, Bilinear_a::diag, kvolHalo, nc, and Bilinear_a::offd.

Here is the caller graph for this function:

GLeft()

void GLeft ( Complex_f out[4], const Complex_f G[2], const Complex_f X[4] )

inlinestatic

Multiplies \( X_{\mu\nu}\) by a gauge field from the left.

Parameters

[out]	out	Result
[in]	X	\(X_{\mu\nu}(x)\)
[in]	G	Gauge field

Definition at line 495 of file clover.c.

                                                                                     {
   out[0]=G[0]*X[0]+G[1]*X[2];
   out[1]=G[0]*X[1]+G[1]*X[3];
   out[2]=-conj(G[1])*X[0]+conj(G[0])*X[2];
   out[3]=-conj(G[1])*X[1]+conj(G[0])*X[3];
   return;
}

References Complex_f, and conj().

Here is the call graph for this function:

Here is the caller graph for this function:

GRight()

void GRight ( Complex_f out[4], const Complex_f G[2], const Complex_f X[4] )

inlinestatic

Multiplies \( X_{\mu\nu}\) by a gauge field from the right.

Parameters

[out]	out	Result
[in]	X	\(X_{\mu\nu}(x)\)
[in]	G	Gauge field

Definition at line 510 of file clover.c.

                                                                                      {
   out[0]=G[0]*X[0]-conj(G[1])*X[1];
   out[1]=G[1]*X[0]+conj(G[0])*X[1];
   out[2]=G[0]*X[2]-conj(G[1])*X[3];
   out[3]=G[1]*X[2]+conj(G[0])*X[3];
   return;
}

References Complex_f, and conj().

Here is the call graph for this function:

Here is the caller graph for this function:

GSandwich()

void GSandwich ( Complex_f out[4], Complex_f tmp[4], const Complex_f Gl[2], const Complex_f X[4], const Complex_f Gr[2] )

inlinestatic

Multiplies \( X_{\mu\nu}\) by a gauge field from the left and the right.

Parameters

[out]	out	Result
[in,out]	tmp	Buffer for intermediate result. Passing as an argument to reduce register pressure.
[in]	X	\(X_{\mu\nu}(x)\)
[in]	Gl,Gr	Left/Right Gauge fields

Postcondition

Contents of tmp will be changed during the call

Definition at line 528 of file clover.c.

                                                                                                                                  {
   GRight(tmp,Gr,X);
   GLeft(out,Gl,tmp);
   return;
}

References Complex_f, GLeft(), and GRight().

Here is the call graph for this function:

Here is the caller graph for this function: