cuclover.cu

Go to the documentation of this file.

#include <clover.h>
 
//CUDA Device code
namespace Device{
   template <typename T>
      __device__ void cuByGenLeft(T a[nc],const unsigned short gen){
         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;
      }

   template <typename T>
      __device__ void cuByGenRight(T a[nc],const unsigned short gen){
         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;
      }

   template <typename T>
      __device__ void Half_Leaf(complex<T> Leaves[nc], complex<T> *u11t, complex<T> *u12t, complex<T> a[nc], unsigned int *iu,\
            unsigned int *id, const unsigned int i, const unsigned short mu, const unsigned short nu, const unsigned short leaf){
 
         unsigned int ind; unsigned int double_ind=id[nu*kvol+ind];
         switch(leaf){
            case(0):
               a[0]=u11t[i+kvolHalo*mu]; a[1]=u12t[i+kvolHalo*mu];
               ind = iu[mu*kvol+i]; 

               Leaves[0]=a[0]*u11t[ind+kvolHalo*nu]-a[1]*conj(u12t[ind+kvolHalo*nu]);
               Leaves[1]=a[0]*u12t[ind+kvolHalo*nu]+a[1]*conj(u11t[ind+kvolHalo*nu]);
               break;
            case(1):
                  ind = id[mu*kvol+i];
                  a[0]=u11t[i+kvolHalo*nu]; a[1]=u12t[i+kvolHalo*nu];
                  //Awkward index...
                  double_ind=iu[nu*kvol+ind];
                  Leaves[0]=a[0]*conj(u11t[double_ind+kvolHalo*mu])+a[1]*conj(u12t[double_ind+kvolHalo*mu]);
                  Leaves[1]=-a[0]*u12t[double_ind+kvolHalo*mu]+a[1]*u11t[double_ind+kvolHalo*mu];
               break;
            case(2):
               ind = id[nu*kvol+i];
               //Daggered. So Conj what goes into a[0] and negate what goes into a[1]
               a[0]=conj(u11t[ind+kvolHalo*nu]); a[1]=-u12t[ind+kvolHalo*nu];

               Leaves[0]=a[0]*u11t[ind+kvolHalo*mu]-a[1]*conj(u12t[ind+kvolHalo*mu]);
               //Don't forget negatiion of second term was handled earlier!
               Leaves[1]=a[0]*u12t[ind+kvolHalo*mu]+a[1]*conj(u11t[ind+kvolHalo*mu]);
               break;
            case(3):
                  ind  =  id[i+kvol*mu];
                  //Daggered. So Conj what goes into a[0] and negate what goes into a[1]
                  a[0]=conj(u11t[ind+kvolHalo*mu]); a[1]=-u12t[ind+kvolHalo*mu];
                  //Another awkward index
                  double_ind=id[nu*kvol+ind];

                  Leaves[0]=a[0]*conj(u11t[double_ind+kvolHalo*nu])+a[1]*conj(u12t[double_ind+kvolHalo*nu]);
                  Leaves[1]=-a[0]*u12t[double_ind+kvolHalo*nu]+a[1]*u11t[double_ind+kvolHalo*nu];
               break;
         }
         return;
      }

   template <typename T>
      __device__ void Leaf(complex<T> *u11t, complex<T> *u12t, complex<T> Leaves[nc],\
            unsigned int *iu, unsigned int *id, unsigned int i,const unsigned short mu,\
            const unsigned short nu,const unsigned short leaf){
         complex<T> a[nc];
         Half_Leaf(Leaves,u11t,u12t,a,iu,id,i,mu,nu,leaf);
         unsigned int ind,double_ind;
         switch(leaf){
            case(0)://Need braces for strict C++ standard compliance as a variable is declared in the case
               {
                  ind = iu[nu*kvol+i]; 
                  a[0]=Leaves[0]*conj(u11t[ind+kvolHalo*mu])+Leaves[1]*conj(u12t[ind+kvolHalo*mu]);
                  a[1]=-Leaves[0]*u12t[ind+kvolHalo*mu]+Leaves[1]*u11t[ind+kvolHalo*mu];

                  Leaves[0]=a[0]*conj(u11t[i+kvolHalo*nu])+a[1]*conj(u12t[i+kvolHalo*nu]);
                  Leaves[1]=-a[0]*u12t[i+kvolHalo*nu]+a[1]*u11t[i+kvolHalo*nu];
 
               }
               break;
            case(1):
               ind = id[mu*kvol+i];

               a[0]=Leaves[0]*conj(u11t[ind+kvolHalo*nu])+Leaves[1]*conj(u12t[ind+kvolHalo*nu]);
               a[1]=-Leaves[0]*u12t[ind+kvolHalo*nu]+Leaves[1]*u11t[ind+kvolHalo*nu];

               Leaves[0]=a[0]*u11t[ind+kvolHalo*mu]-a[1]*conj(u12t[ind+kvolHalo*mu]);
               Leaves[1]=a[0]*u12t[ind+kvolHalo*mu]+a[1]*conj(u11t[ind+kvolHalo*mu]);
               break;
            case(2):
                  ind = id[nu*kvol+i]; 
                  double_ind=iu[mu*kvol+ind];
                  a[0]=Leaves[0]*u11t[double_ind+kvolHalo*nu]-Leaves[1]*conj(u12t[double_ind+kvolHalo*nu]);
                  a[1]=Leaves[0]*u12t[double_ind+kvolHalo*nu]+Leaves[1]*conj(u11t[double_ind+kvolHalo*nu]);

                  Leaves[0]=a[0]*conj(u11t[i+kvolHalo*mu])+a[1]*conj(u12t[i+kvolHalo*mu]);
                  Leaves[1]=-a[0]*u12t[i+kvolHalo*mu]+a[1]*u11t[i+kvolHalo*mu];
 
               break;
            case(3)://Need braces for strict C++ standard compliance as a variable is declared in the case
                  ind = id[nu*kvol+i]; 
                  double_ind=id[mu*kvol+ind];

                  a[0]=Leaves[0]*u11t[double_ind+kvolHalo*mu]-Leaves[1]*conj(u12t[double_ind+kvolHalo*mu]);
                  a[1]=Leaves[0]*u12t[double_ind+kvolHalo*mu]+Leaves[1]*conj(u11t[double_ind+kvolHalo*mu]);

                  Leaves[0]=a[0]*u11t[ind+kvolHalo*nu]-a[1]*conj(u12t[ind+kvolHalo*nu]);
                  Leaves[1]=a[0]*u12t[ind+kvolHalo*nu]+a[1]*conj(u11t[ind+kvolHalo*nu]);
               break;
         }
         return;
      }

   template <typename T>
      __device__ void cuGLeft(T out[4],const T G[2], const T X[4]){
         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;
      }

   template <typename T>
      __device__ void cuGRight(T out[4],const T G[2], const T X[4]){
         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;
      }

   template <typename T>
      __device__ void cuGSandwich(T out[4],T tmp[4], const T Gl[2], const T X[4],const T Gr[2]){
         cuGRight(tmp,Gr,X);
         cuGLeft(out,Gl,tmp);
         return;
      }

   __device__ void GetBilinear(Complex_f Z[nc*nc], Bilinear_a Xmn,unsigned int ind){
      Z[0]=Xmn.diag[ind]; Z[3]=Xmn.diag[ind+kvolHalo];
      Z[1]=Xmn.offd[ind]; Z[2]=conj(Z[1]);
      return;
   }
}

namespace Kernels{
   using namespace Device;
   template <typename T>
      __global__ void Half_Leaves(complex<T> *hLeaves0,complex<T> *hLeaves1,complex<T> *u11t,complex<T> *u12t,\
            unsigned int *iu,unsigned int *id,const unsigned short mu,const unsigned short nu){
         const volatile int gsize = gridDim.x*gridDim.y*gridDim.z;
         const volatile int bsize = blockDim.x*blockDim.y*blockDim.z;
         const volatile int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
         const volatile int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
         const volatile int gthreadId= blockId * bsize+bthreadId;
 
         complex<T> Leaves[nc], a[nc];
         for(unsigned short leaf=0;leaf<ndim;leaf++)
            for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
               Half_Leaf(Leaves,u11t,u12t,a,iu,id,i,mu,nu,leaf);
               hLeaves0[i+kvol*leaf]=Leaves[0]; hLeaves1[i+kvol*leaf]=Leaves[1];
            }
         return;
      }

   template <typename T>
      __global__  void Full_Clover(complex<T> *clover1, complex<T> *clover2,\
            complex<T> *u11t, complex<T> *u12t, unsigned int *iu, unsigned int *id, int mu, int nu){
         const volatile int gsize = gridDim.x*gridDim.y*gridDim.z;
         const volatile int bsize = blockDim.x*blockDim.y*blockDim.z;
         const volatile int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
         const volatile int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
         const volatile int gthreadId= blockId * bsize+bthreadId;
 
         complex<T> Leaves[2];
         for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
            clover1[i]=0;clover2[i]=0;
            for(unsigned short leaf=0;leaf<ndim;leaf++)
            {
               //Pointer arithemetic on the leaves.
               Leaf(u11t,u12t,Leaves,iu,id,i,mu,nu,leaf);
               clover1[i]+=Leaves[0]; clover2[i]+=Leaves[1];
            }

            clover1[i]=clover1[i].imag();    clover1[i]*=(1.0f/4.0f);
            //    clover1[i]=clover1[1].imag()/4.0f;

            clover2[i]+=clover2[i];             clover2[i]*=(-I_f/8.0f);
         }
         return;
      }

   template <typename T>
      __global__ void cuCalcXmunu(Bilinear_a Xmunu, const complex<T> *X1, const complex<T> *X2,
            const complex<T> *sigval, const unsigned short *sigin,const unsigned short clov){
         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;
      }

   template <typename T>
      __global__ void 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){
         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;
      }
 
   //Clover multiplication
   template <typename T>
      __global__ void ByClover(complex<T> *phi, complex<T> *r, complex<T> *clover1, complex<T> *clover2, complex<T> *sigval, const float akappa, unsigned short *sigin, bool dag){
         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+=bsize*gsize){
            //Prefetched r and Phi array
            complex<T> phi_s[ngorkov][nc];
#pragma unroll
            for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++)
               for(unsigned short c=0; c<nc; c++){
                  phi_s[igorkov][c]=0;
               }
            complex<T> r_s[nc];
            complex<T> clov_s[nc];
#pragma unroll
            for(unsigned short clov=0;clov<6;clov++){
               clov_s[0]=clover1[clov*kvol+i]; clov_s[1]=clover2[clov*kvol+i];
               for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
                  //Mod 4 done bitwise. In general n mod 2^m = n & (2^m-1)
                  const unsigned short idirac = igorkov&3;
                  const unsigned short sind = (igorkov<4) ? sigin[clov*ndirac+idirac] : sigin[clov*ndirac+idirac]+4;
#pragma unroll
                  for(unsigned short c=0; c<nc; c++)
                     r_s[c]= r[i+kvolHalo*(sind*nc+c)];
                  phi_s[igorkov][0]+=sigval[clov*ndirac+idirac]*(creal(clov_s[0])*r_s[0]+clov_s[1]*r_s[1]);
                  //Clover is in the Lie Algebra, not Lie group. So signs are correct here.
                  phi_s[igorkov][1]+=sigval[clov*ndirac+idirac]*(conj(clov_s[1])*r_s[0]-creal(clov_s[0])*r_s[1]);
               }
            }
#pragma unroll
            for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++)
               for(unsigned short c=0; c<nc; c++){
                  //dag is just to do with the output layout and if it has a halo
                  if(dag)
                     phi[i+kvol*(nc*igorkov+c)]+=akappa*phi_s[igorkov][c];
                  else
                     phi[i+kvolHalo*(nc*igorkov+c)]+=akappa*phi_s[igorkov][c];
               }
         }
         return;
      }

   template <typename T>
      __global__ void HbyClover(complex<T> *phi, complex<T> *r, complex<T> *clover1, complex<T> *clover2,complex<T> *sigval, const float akappa, unsigned short *sigin,bool dag){
         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+=bsize*gsize){
            //Prefetched r and Phi array
            complex<T> phi_s[ndirac*nc];
#pragma unroll
            for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc)
               for(unsigned short c=0; c<nc; c++){
                  phi_s[idirac+c]=0;
               }
            complex<T> r_s[nc]; complex<T> clov_s[nc];
#pragma unroll
            for(unsigned short clov=0;clov<6;clov++){
               clov_s[0]=clover1[clov*kvol+i]; clov_s[1]=clover2[clov*kvol+i];
               for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
                  const unsigned short sind = sigin[clov*ndirac+(idirac>>1)] << (nc-1);
#pragma unroll
                  for(unsigned short c=0; c<nc; c++){
                     r_s[c]= r[i+kvolHalo*(sind+c)];
                  }
                  const complex<T> sig=sigval[clov*ndirac+(idirac>>1)];
                  phi_s[idirac+0]+=sig*(creal(clov_s[0])*r_s[0]+clov_s[1]*r_s[1]);
                  //Clover is in the Lie Algebra, not Lie group. So signs are correct here.
                  phi_s[idirac+1]+=sig*(conj(clov_s[1])*r_s[0]-creal(clov_s[0])*r_s[1]);
               }
            }
#pragma unroll
            for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc)
               for(unsigned short c=0; c<nc; c++)
                  //dag is just to do with the output layout and if it has a halo
                  if(dag)
                     phi[i+kvol*(c+idirac)]+=akappa*phi_s[idirac+c];
                  else
                     phi[i+kvolHalo*(c+idirac)]+=akappa*phi_s[idirac+c];
         }
         return;
      }
}
using namespace Kernels;
 
//Calling Wrappers
//This gets called by C so cannot be templated...
int cuClover(Complex_f *clover[nc],Complex_f *ut[nc], unsigned int *iu, unsigned int *id){
   const char funcname[]="cuClover";
#ifdef _DEBUG
   cudaMallocManaged((void **)&clover[0],6*kvol*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&clover[1],6*kvol*sizeof(Complex_f),cudaMemAttachGlobal);
#else
   cudaMallocAsync((void **)&clover[0],6*kvol*sizeof(Complex_f),streams[0]);
   cudaMallocAsync((void **)&clover[1],6*kvol*sizeof(Complex_f),streams[1]);
#endif
   for(unsigned short mu=0;mu<ndim-1;mu++)
      for(unsigned short nu=mu+1;nu<ndim;nu++)
         if(mu!=nu){
            //Clover index
            unsigned short clov = (mu==0) ? nu-1 :mu+nu;
            //Allocate clover memory
            //Note that the clover is completely local, so doesn't need a halo for MPI
            Full_Clover<<<dimGrid,dimBlock,0,streams[clov]>>>(clover[0]+clov*kvol,clover[1]+clov*kvol,\
                  ut[0],ut[1],iu,id,mu,nu);
         }
   cudaDeviceSynchronise();
   return 0;
}
void cuByClover(Complex *phi, Complex *r, Complex *clover[nc],Complex *sigval,const float akappa, unsigned short *sigin, bool dag){
   Kernels::ByClover<<<dimGrid,dimBlock>>>(phi,r,clover[0],clover[1],sigval,akappa,sigin,dag);
}
void cuHbyClover(Complex *phi, Complex *r, Complex *clover[nc],Complex *sigval, const float akappa, unsigned short *sigin, bool dag){
   Kernels::HbyClover<<<dimGrid,dimBlock>>>(phi,r,clover[0],clover[1],sigval,akappa,sigin,dag);
}  
void cuByClover_f(Complex_f *phi, Complex_f *r, Complex_f *clover[nc],Complex_f *sigval, const float akappa, unsigned short *sigin,bool dag){
   Kernels::ByClover<<<dimGrid,dimBlock>>>(phi,r,clover[0],clover[1],sigval,akappa,sigin,dag);
}
void cuHbyClover_f(Complex_f *phi, Complex_f *r, Complex_f *clover[nc],Complex_f *sigval, const float akappa, unsigned short *sigin,bool dag){
   Kernels::HbyClover<<<dimGrid,dimBlock>>>(phi,r,clover[0],clover[1],sigval,akappa,sigin,dag);
}
 
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){
   //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;
}
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){
   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;
}