su2hmc
Loading...
Searching...
No Matches
Fermion matrix products
Collaboration diagram for Fermion matrix products:

Topics

 Clover Multiplication routines

Functions

int Dslash (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa)
 Evaluates \(\Phi=M r\) in double precision.
int Dslashd (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa)
 Evaluates \(\Phi=M^\dagger r\) in double precision.
int Hdslash (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa)
 Evaluates \(\Phi=M r\) in double precision.
int Hdslashd (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa)
 Evaluates \(\Phi=M^\dagger r\) in double precision.
int Dslash_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa)
 Evaluates \(\Phi=M r\) in single precision.
int Dslashd_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa)
 Evaluates \(\Phi=M^\dagger r\) in single precision.
int Hdslash_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa)
 Evaluates \(\Phi=M r\) in single precision.
int Hdslashd_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa)
 Evaluates \(\Phi=M^\dagger r\) in single precision.
void cuDslash (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M r\) in double precision.
void cuDslashd (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.
void cuHdslash (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.
void cuHdslashd (Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.
void cuDslash_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.
void cuDslashd_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.
void cuHdslash_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M r\) in single precision.
void cuHdslashd_f (Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
 GPU calling wrapper for \(\Phi=M^\dagger r\) in single precision.
template<typename T>
__global__ void Kernels::cuDslash (complex< T > *phi, complex< T > *r, complex< T > *u11t, complex< T > *u12t, const unsigned int *iu, const unsigned int *id, complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa)
 Evaluates \(\Phi=M r\).
template<typename T>
__global__ void Kernels::cuDslashd (complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, const unsigned int *iu, const unsigned int *id, complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa)
 Evaluates \(\Phi=M^\dagger r\).
template<typename T>
__global__ void Kernels::cuHdslash (complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, unsigned int *iu, unsigned int *id, __constant__ complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const __grid_constant__ float akappa)
 Evaluates \(\Phi=Mr\) using up/down partitioning.
template<typename T>
__global__ void Kernels::cuHdslashd (complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, unsigned int *iu, unsigned int *id, __constant__ complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const __grid_constant__ float akappa)
 Evaluates \(\Phi=M^\dagger r\) using up/down partitioning.

Detailed Description

Function Documentation

◆ cuDslash() [1/2]

void cuDslash ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
Complex_f jqq,
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 609 of file cumatrices.cu.

611 {
612 const char funcname[] = "Dslash";
613 int cuCpyStat=0;
614 for(unsigned short j=0;j<nc*ngorkov;j++)
615 if((cuCpyStat=cudaMemcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex),cudaMemcpyDefault))){
616 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
617 CPYERROR,funcname,cuCpyStat);
618 exit(cuCpyStat);
619 }
620 Kernels::cuDslash<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],jqq,akappa);
621 return;
622}
CPYERROR
#define CPYERROR
Copy failed.
Definition errorcodes.h:66
Kernels::cuDslash
__global__ void cuDslash(complex< T > *phi, complex< T > *r, complex< T > *u11t, complex< T > *u12t, const unsigned int *iu, const unsigned int *id, complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa)
Evaluates .
Definition cumatrices.cu:48
nc
#define nc
Colours.
Definition sizes.h:182
ngorkov
#define ngorkov
Gor'kov indices.
Definition sizes.h:190
kvol
#define kvol
Sublattice volume.
Definition sizes.h:163
Complex
#define Complex
Double precision complex number.
Definition sizes.h:64
kvolHalo
#define kvolHalo
Subvolume + halo size.
Definition sizes.h:234

References Complex, Complex_f, CPYERROR, Kernels::cuDslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuDslash() [2/2]

template<typename T>
__global__ void Kernels::cuDslash ( complex< T > * phi,
complex< T > * r,
complex< T > * u11t,
complex< T > * u12t,
const unsigned int * iu,
const unsigned int * id,
complex< T > gamval[20],
const unsigned short gamin[16],
const T * dk4m,
const T * dk4p,
const Complex_f jqq,
const float akappa )

Evaluates \(\Phi=M r\).

Parameters
[in,out]phiThe product
[in]rThe array being acted on by M
[in]u11t,u12tGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk4m,dk4p\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result added to phi

Definition at line 48 of file cumatrices.cu.

49 {
50 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
51 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
52 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
53 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
54 const unsigned int gthreadId= blockId * bsize+bthreadId;
55
56 for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
57 complex<T> ru[nc]; complex<T> rd[nc];
58 complex<T> rgu[nc]; complex<T> rgd[nc];
59 complex<T> phi_s[ngorkov*nc];
60 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
61 unsigned short igork = ((idirac>>1)+4)<<1;
62 unsigned int ind_d =4*ndirac+(idirac>>1);
63 complex<T> a_1=conj(jqq)*gamval[ind_d];
64 //We subtract a_2, hence the minus
65 complex<T> a_2=-jqq*gamval[ind_d];
66 ind_d=i+kvolHalo*(idirac); unsigned int ind_g=i+kvolHalo*(igork);
67 phi_s[idirac]=phi[ind_d]+a_1*r[ind_g];
68 phi_s[igork]=phi[ind_g]+a_2*r[ind_d];
69 ind_d+=kvolHalo; ind_g+=kvolHalo;
70 phi_s[idirac+1]=phi[ind_d]+a_1*r[ind_g];
71 phi_s[igork+1]=phi[ind_g]+a_2*r[ind_d];
72 }
73 complex<T> u11s; complex<T> u12s;
74 complex<T> u11sd; complex<T> u12sd;
75 unsigned int ind;
76 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
77#ifndef NO_SPACE
78 for(unsigned short mu = 0; mu <3; mu++){
79 ind = i+kvol*mu;
80 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
81 ind = i+kvolHalo*mu;
82 u11s=u11t[ind]; u12s=u12t[ind];
83 ind = did+kvolHalo*mu;
84 u11sd=u11t[ind]; u12sd=u12t[ind];
85 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
86 unsigned short idirac=igorkov&3;
87 unsigned short gind=mu*ndirac+idirac;
88 const complex<T> gam=gamval[gind];
89 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing in the dirac term.
90 unsigned short igork1 = (igorkov<4) ? gamin[gind] : gamin[gind]+4;
91 for(unsigned short c=0;c<nc;c++){
92 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
93 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
94 }
95 //Wilson + Dirac term in that order. Definitely easier
96 phi_s[igorkov*nc]+=-akappa*(u11s*ru[0]+ u12s*ru[1]+\
97 conj(u11sd)*rd[0]- u12sd*rd[1]);
98 //Dirac term
99 phi_s[igorkov*nc]+=gam*(u11s*rgu[0]+ u12s*rgu[1]-\
100 conj(u11sd)*rgd[0]+ u12sd*rgd[1]);
101
102 phi_s[igorkov*nc+1]+=-akappa*(-conj(u12s)*ru[0]+ conj(u11s)*ru[1]+\
103 conj(u12sd)*rd[0]+ u11sd*rd[1]);
104 //Dirac term
105 phi_s[igorkov*nc+1]+=gam*(-conj(u12s)*rgu[0]+ conj(u11s)*rgu[1]-\
106 conj(u12sd)*rgd[0]- u11sd*rgd[1]);
107 }
108 }
109 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
110 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
111 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
112#endif
113#ifndef NO_TIME
114 ind=i+kvolHalo*3;
115 u11s=u11t[ind]; u12s=u12t[ind];
116 const T dk4ms=dk4m[i]; const T dk4ps=dk4p[i];
117 ind=i+kvol*3;
118 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
119 ind=did+kvolHalo*3;
120 u11sd=u11t[ind]; u12sd=u12t[ind];
121 const T dk4msd=dk4m[did]; const T dk4psd=dk4p[did];
122 for(unsigned short igorkov=0;igorkov<ndirac;igorkov++){
123 unsigned short igork1 = gamin[3*ndirac+igorkov];
124 for(unsigned short c=0;c<nc;c++){
125 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
126 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
127 }
128 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
129 phi_s[igorkov*nc]+=
130 -dk4ps*(u11s*(ru[0]-rgu[0]) +u12s*(ru[1]-rgu[1]))
131 -dk4msd*(conj(u11sd)*(rd[0]+rgd[0]) -u12sd *(rd[1]+rgd[1]));
132 phi[i+kvolHalo*(igorkov*nc)]=phi_s[igorkov*nc];
133
134 phi_s[igorkov*nc+1]+=
135 -dk4ps*(-conj(u12s)*(ru[0]-rgu[0]) +conj(u11s)*(ru[1]-rgu[1]))
136 -dk4msd*(conj(u12sd)*(rd[0]+rgd[0]) +u11sd *(rd[1]+rgd[1]));
137 phi[i+kvolHalo*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
138 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
139 //the FORTRAN code did it.
140 igork1 += 4;
141 //And the gorkov terms. Note that dk4p and dk4m swap positions compared to the above
142 for(unsigned short c=0;c<nc;c++){
143 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
144 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
145 }
146 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
147 phi_s[igorkovPP*nc]+=-dk4ms*(u11s*(ru[0]-rgu[0])+ u12s*(ru[1]-rgu[1]))-
148 dk4psd*(conj(u11sd)*(rd[0]+rgd[0])- u12sd*(rd[1]+rgd[1]));
149 phi[i+kvolHalo*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
150
151 phi_s[igorkovPP*nc+1]+=-dk4ms*(conj(-u12s)*(ru[0]-rgu[0]) +conj(u11s)*(ru[1]-rgu[1]))
152 -dk4psd*(conj(u12sd)*(rd[0]+rgd[0]) +u11sd*(rd[1]+rgd[1]));
153 phi[i+kvolHalo*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
154 }
155#endif
156 }
157 }
conj
__device__ __forceinline__ T conj(const T &z)
Complex Conjugation.
Definition cusu2hmc.cu:33
ndirac
#define ndirac
Dirac indices.
Definition sizes.h:186

References Complex_f, conj(), kvol, kvolHalo, nc, ndirac, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuDslash_f()

void cuDslash_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
Complex_f jqq,
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 668 of file cumatrices.cu.

670 {
671 const char funcname[] = "Dslash_f";
672 int cuCpyStat=0;
673 for(unsigned short j=0;j<nc*ngorkov;j++)
674 if((cuCpyStat=cudaMemcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex_f),cudaMemcpyDefault))){
675 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
676 CPYERROR,funcname,cuCpyStat);
677 exit(cuCpyStat);
678 }
679 Kernels::cuDslash<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],jqq,akappa);
680 return;
681}
Complex_f
#define Complex_f
Single precision complex number.
Definition sizes.h:62

References Complex_f, CPYERROR, Kernels::cuDslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuDslashd() [1/2]

void cuDslashd ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
Complex_f jqq,
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 623 of file cumatrices.cu.

625 {
626 const char funcname[] = "Dslashd";
627 int cuCpyStat=0;
628 for(unsigned short j=0;j<nc*ngorkov;j++)
629 if((cuCpyStat=cudaMemcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex),cudaMemcpyDefault))){
630 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
631 CPYERROR,funcname,cuCpyStat);
632 exit(cuCpyStat);
633 }
634 Kernels::cuDslashd<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],jqq,akappa);
635 return;
636}
Kernels::cuDslashd
__global__ void cuDslashd(complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, const unsigned int *iu, const unsigned int *id, complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa)
Evaluates .
Definition cumatrices.cu:175

References Complex, Complex_f, CPYERROR, Kernels::cuDslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuDslashd() [2/2]

template<typename T>
__global__ void Kernels::cuDslashd ( complex< T > * phi,
const complex< T > * r,
const complex< T > * u11t,
const complex< T > * u12t,
const unsigned int * iu,
const unsigned int * id,
complex< T > gamval[20],
const unsigned short gamin[16],
const T * dk4m,
const T * dk4p,
const Complex_f jqq,
const float akappa )

Evaluates \(\Phi=M^\dagger r\).

Parameters
[in,out]phiThe product
[in]rThe array being acted on by M
[in]u11t,u12tGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk4m,dk4p\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result added to phi

Definition at line 175 of file cumatrices.cu.

176 {
177 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
178 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
179 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
180 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
181 const unsigned int gthreadId= blockId * bsize+bthreadId;
182
183 for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
184 complex<T> ru[nc]; complex<T> rd[nc];
185 complex<T> rgu[nc]; complex<T> rgd[nc];
186 complex<T> phi_s[ngorkov*nc];
187 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
188 unsigned short igork = ((idirac>>1)+4)<<1;
189 unsigned int ind_d =4*ndirac+(idirac>>1);
190 complex<T> a_1=-conj(jqq)*gamval[ind_d];
191 complex<T> a_2=jqq*gamval[ind_d];
192 ind_d=i+kvol*(idirac); unsigned int ind_g=i+kvol*(igork);
193 phi_s[idirac]=phi[ind_d]+a_1*r[ind_g];
194 phi_s[igork]=phi[ind_g]+a_2*r[ind_d];
195 ind_d+=kvol; ind_g+=kvol;
196 phi_s[idirac+1]=phi[ind_d]+a_1*r[ind_g];
197 phi_s[igork+1]=phi[ind_g]+a_2*r[ind_d];
198 }
199 complex<T> u11s; complex<T> u12s;
200 complex<T> u11sd; complex<T> u12sd;
201 unsigned int ind;
202 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
203#ifndef NO_SPACE
204 for(unsigned short mu = 0; mu <3; mu++){
205 ind = i+kvol*mu;
206 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
207 ind = i+kvolHalo*mu;
208 u11s=u11t[ind]; u12s=u12t[ind];
209 ind = did+kvolHalo*mu;
210 u11sd=u11t[ind]; u12sd=u12t[ind];
211 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
212 unsigned short idirac=igorkov&3;
213 const complex<T> gam=gamval[mu*ndirac+idirac];
214 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing.
215 unsigned short igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
216 for(unsigned short c=0;c<nc;c++){
217 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
218 rgd[c]=r[did+kvolHalo*(igork1*nc+c)]; rgu[c]=r[uid+kvolHalo*(igork1*nc+c)];
219 }
220 //Wilson + Dirac term in that order. Definitely easier
221 phi_s[igorkov*nc]-= akappa*(u11s*ru[0] +u12s*ru[1]
222 +conj(u11sd)*rd[0] -u12sd *rd[1]);
223
224 //Dirac term
225 phi_s[igorkov*nc]-=gam* (u11s*rgu[0] +u12s*rgu[1]
226 -conj(u11sd)*rgd[0] +u12sd *rgd[1]);
227
228 phi_s[igorkov*nc+1]-= akappa*(-conj(u12s)*ru[0] +conj(u11s)*ru[1]
229 +conj(u12sd)*rd[0] +u11sd *rd[1]);
230 //Dirac term
231 phi_s[igorkov*nc+1]-=gam* (-conj(u12s)*rgu[0] +conj(u11s)*rgu[1]
232 -conj(u12sd)*rgd[0] -u11sd *rgd[1]);
233
234 }
235 }
236#endif
237 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
238 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
239 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
240 //Under dagger, dk4p and dk4m get swapped and the dirac component flips sign.
241#ifndef NO_TIME
242 ind=i+kvolHalo*3;
243 u11s=u11t[ind]; u12s=u12t[ind];
244 const T dk4ms=dk4m[i]; const T dk4ps=dk4p[i];
245 ind = i+kvol*3;
246 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
247 ind=did+kvolHalo*3;
248 u11sd=u11t[ind]; u12sd=u12t[ind];
249 const T dk4msd=dk4m[did]; const T dk4psd=dk4p[did];
250 for(unsigned short igorkov=0; igorkov<ndirac; igorkov++){
251 unsigned short igork1 = gamin[3*ndirac+igorkov];
252 for(unsigned short c=0;c<nc;c++){
253 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
254 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
255 }
256 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
257 phi_s[igorkov*nc]+=
258 -dk4ms*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
259 -dk4psd*(conj(u11sd)*(rd[0]-rgd[0]) -u12sd *(rd[1]-rgd[1]));
260 phi[i+kvol*(igorkov*nc)]=phi_s[igorkov*nc];
261
262 phi_s[igorkov*nc+1]+=
263 -dk4ms*(-conj(u12s)*(ru[0]+rgu[0]) +conj(u11s)*(ru[1]+rgu[1]))
264 -dk4psd*(conj(u12sd)*(rd[0]-rgd[0]) +u11sd *(rd[1]-rgd[1]));
265 phi[i+kvol*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
266 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
267 //the FORTRAN code did it.
268 igork1 += 4;
269 for(unsigned short c=0;c<nc;c++){
270 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
271 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
272 }
273 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
274 phi_s[igorkovPP*nc]+=-dk4ps*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
275 -dk4msd*(conj(u11sd)*(rd[0]-rgd[0]) -u12sd*(rd[1]-rgd[1]));
276 phi[i+kvol*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
277
278 phi_s[igorkovPP*nc+1]+=dk4ps*(conj(u12s)*(ru[0]+rgu[0]) -conj(u11s)*(ru[1]+rgu[1]))
279 -dk4msd*(conj(u12sd)*(rd[0]-rgd[0]) +u11sd*(rd[1]-rgd[1]));
280 phi[i+kvol*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
281 }
282#endif
283 }
284 }

References Complex_f, conj(), kvol, kvolHalo, nc, ndirac, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuDslashd_f()

void cuDslashd_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
Complex_f jqq,
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 682 of file cumatrices.cu.

684 {
685 const char funcname[] = "Dslashd_f";
686 int cuCpyStat=0;
687 for(unsigned short j=0;j<nc*ngorkov;j++)
688 if((cuCpyStat=cudaMemcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex_f),cudaMemcpyDefault))){
689 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
690 CPYERROR,funcname,cuCpyStat);
691 exit(cuCpyStat);
692 }
693 Kernels::cuDslashd<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],jqq,akappa);
694 return;
695}

References Complex_f, CPYERROR, Kernels::cuDslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslash() [1/2]

void cuHdslash ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 637 of file cumatrices.cu.

639 {
640 const char funcname[] = "Hdslash";
641 int cuCpyStat=0;
642 for(unsigned short j=0;j<nc*ndirac;j++)
643 if((cuCpyStat=cudaMemcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex),cudaMemcpyDefault))){
644 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
645 CPYERROR,funcname,cuCpyStat);
646 exit(cuCpyStat);
647 }
648 Kernels::cuHdslash<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],akappa);
649 return;
650}
Kernels::cuHdslash
__global__ void cuHdslash(complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, unsigned int *iu, unsigned int *id, __constant__ complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const __grid_constant__ float akappa)
Evaluates using up/down partitioning.
Definition cumatrices.cu:302

References Complex, CPYERROR, Kernels::cuHdslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslash() [2/2]

template<typename T>
__global__ void Kernels::cuHdslash ( complex< T > * phi,
const complex< T > * r,
const complex< T > * u11t,
const complex< T > * u12t,
unsigned int * iu,
unsigned int * id,
__constant__ complex< T > gamval[20],
const unsigned short gamin[16],
const T * dk4m,
const T * dk4p,
const __grid_constant__ float akappa )

Evaluates \(\Phi=Mr\) using up/down partitioning.

Parameters
[in,out]phiThe product
[in]rThe array being acted on by M
[in]u11t,u12tGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk4m,dk4p\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
Postcondition
Result added to phi

Definition at line 302 of file cumatrices.cu.

303 {
304 /*
305 * Half Dslash T precision
306 */
307 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
308 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
309 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
310 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
311 const unsigned int gthreadId= blockId * bsize+bthreadId;
312
313 //Right. Time to prefetch
314 complex<T> ru[2]; complex<T> rd[2];
315 complex<T> rgu[2]; complex<T> rgd[2];
316 complex<T> phi_s[ndirac*nc];
317 for(unsigned int i=gthreadId;i<kvol;i+=bsize*gsize){
318#pragma unroll
319 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
320#pragma unroll
321 for(unsigned short c=0; c<nc; c++)
322 //NOTE: idirac is increasing by nc each time. So should be read as idirac*nc
323 phi_s[idirac+c]=phi[i+kvolHalo*(c+idirac)];
324
325 //#pragma unroll
326 for(unsigned short mu = 0; mu <ndim; mu++){
327 unsigned int ind=i+kvolHalo*mu;
328 const complex<T> u11s=u11t[ind]; const complex<T> u12s=u12t[ind];
329 ind = i+kvol*mu;
330 const int did=id[ind]; const int uid = iu[ind];
331 ind=did+kvolHalo*mu;
332 const complex<T> u11sd=u11t[ind]; const complex<T> u12sd=u12t[ind];
333#pragma unroll
334 for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
335 const unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
336#pragma unroll
337 for(unsigned short c=0;c<nc;c++){
338 ind =kvolHalo*(idirac+c);
339 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
340 ind =kvolHalo*(igork1+c);
341 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
342 }
343 //Can manually vectorise with a pragma?
344 //Wilson + Dirac term in that order. Definitely easier
345 //to read when split into different loops, but should be faster this way
346 //Spacelike terms
347 if(mu<3){
348 const complex<T> gam=gamval[mu*ndirac+(idirac>>1)];
349 phi_s[idirac]+=-akappa*(u11s*ru[0]+u12s*ru[1]+\
350 conj(u11sd)*rd[0]-u12sd*rd[1]);
351 //Dirac term
352 phi_s[idirac]+=gam*(u11s*rgu[0]+u12s*rgu[1]-\
353 conj(u11sd)*rgd[0]+ u12sd*rgd[1]);
354
355 phi_s[idirac+1]+=-akappa*(-conj(u12s)*ru[0]+ conj(u11s)*ru[1]+\
356 conj(u12sd)*rd[0]+ u11sd*rd[1]);
357 //Dirac term
358 phi_s[idirac+1]+=gam*(-conj(u12s)*rgu[0]+ conj(u11s)*rgu[1]-\
359 conj(u12sd)*rgd[0]- u11sd*rgd[1]);
360 }
361 //Timelike terms
362 else{
363 const T dk4ms=dk4m[did]; const T dk4ps=dk4p[i];
364 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
365
366 phi_s[idirac+0]-= dk4ps*(u11s*(ru[0]-rgu[0])
367 +u12s*(ru[1]-rgu[1]));
368 phi_s[idirac+0]-= dk4ms*(conj(u11sd)*(rd[0]+rgd[0])
369 -u12sd *(rd[1]+rgd[1]));
370 phi[i+kvolHalo*(0+idirac)]=phi_s[idirac+0];
371
372 phi_s[idirac+1]-= dk4ps*(-conj(u12s)*(ru[0]-rgu[0])
373 +conj(u11s)*(ru[1]-rgu[1]));
374 phi_s[idirac+1]-= dk4ms*(conj(u12sd)*(rd[0]+rgd[0])
375 +u11sd *(rd[1]+rgd[1]));
376 phi[i+kvolHalo*(1+idirac)]=phi_s[idirac+1];
377 }
378 }
379 }
380 }
381 }
ndim
#define ndim
Dimensions.
Definition sizes.h:188

References conj(), kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslash_f()

void cuHdslash_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 696 of file cumatrices.cu.

697 {
698 const char funcname[] = "Hdslash_f";
699 int cuCpyStat=0;
700 for(unsigned short j=0;j<nc*ndirac;j++)
701 if((cuCpyStat=cudaMemcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex_f),cudaMemcpyDefault))){
702 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
703 CPYERROR,funcname,cuCpyStat);
704 exit(cuCpyStat);
705 }
706 const int bsize=dimGrid.x*dimGrid.y*dimGrid.z;
707 const int shareSize= ndim*bsize*nc*sizeof(Complex_f);
708 Kernels::cuHdslash<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],akappa);
709 return;
710}
dimGrid
dim3 dimGrid
Default grid size. First component is normally nt. Second and third depend whatever is needed to get ...
Definition cusu2hmc.cu:27

References Complex_f, CPYERROR, Kernels::cuHdslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslashd() [1/2]

void cuHdslashd ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 651 of file cumatrices.cu.

653 {
654 const char funcname[] = "Hdslashd";
655 //Spacelike term
656 int cuCpyStat=0;
657 for(unsigned short j=0;j<nc*ndirac;j++)
658 if((cuCpyStat=cudaMemcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex),cudaMemcpyDefault))){
659 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
660 CPYERROR,funcname,cuCpyStat);
661 exit(cuCpyStat);
662 }
663 Kernels::cuHdslashd<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],akappa);
664 return;
665}
Kernels::cuHdslashd
__global__ void cuHdslashd(complex< T > *phi, const complex< T > *r, const complex< T > *u11t, const complex< T > *u12t, unsigned int *iu, unsigned int *id, __constant__ complex< T > gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const __grid_constant__ float akappa)
Evaluates using up/down partitioning.
Definition cumatrices.cu:398

References Complex, CPYERROR, Kernels::cuHdslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslashd() [2/2]

template<typename T>
__global__ void Kernels::cuHdslashd ( complex< T > * phi,
const complex< T > * r,
const complex< T > * u11t,
const complex< T > * u12t,
unsigned int * iu,
unsigned int * id,
__constant__ complex< T > gamval[20],
const unsigned short gamin[16],
const T * dk4m,
const T * dk4p,
const __grid_constant__ float akappa )

Evaluates \(\Phi=M^\dagger r\) using up/down partitioning.

Parameters
[in,out]phiThe product
[in]rThe array being acted on by M
[in]u11t,u12tGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk4m,dk4p\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
Postcondition
Result added to phi

Definition at line 398 of file cumatrices.cu.

399 {
400 /*
401 * Half Dslash Dagger T precision
402 */
403 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
404 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
405 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
406 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
407 const unsigned int gthreadId= blockId * bsize+bthreadId;
408
409 //Right. Time to prefetch
410 for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
411 complex<T> phi_s[ndirac*nc];
412#pragma unroll
413 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
414#pragma unroll
415 for(unsigned short c=0; c<nc; c++)
416 //NOTE: idirac is increasing by nc each time. So should be read as idirac*nc
417 phi_s[idirac+c]=phi[i+kvol*(c+idirac)];
418
419 //#pragma unroll
420 for(unsigned short mu = 0; mu <ndim; mu++){
421 unsigned int ind=i+kvolHalo*mu;
422 const complex<T> u11s=u11t[ind]; const complex<T> u12s=u12t[ind];
423 ind = i+kvol*mu;
424 const int did=id[ind]; const int uid = iu[ind];
425 ind=did+kvolHalo*mu;
426 const complex<T> u11sd=u11t[ind]; const complex<T> u12sd=u12t[ind];
427#pragma unroll
428 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc){
429 const unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
430 complex<T> ru[2]; complex<T> rd[2];
431 complex<T> rgu[2]; complex<T> rgd[2];
432#pragma unroll
433 for(unsigned short c=0;c<nc;c++){
434 ind =kvolHalo*(idirac+c);
435 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
436 ind =kvolHalo*(igork1+c);
437 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
438 }
439 //Can manually vectorise with a pragma?
440 //Wilson + Dirac term in that order. Definitely easier
441 //to read when split into different loops, but should be faster this way
442 //Spacelike terms
443 if(mu<3){
444 const complex<T> gam=gamval[mu*ndirac+(idirac>>1)];
445 phi_s[idirac]-=akappa*(u11s*ru[0] +u12s*ru[1]
446 +conj(u11sd)*rd[0] -u12sd *rd[1]);
447 //Dirac term
448 phi_s[idirac]-=gam* (u11s*rgu[0] +u12s*rgu[1]
449 -conj(u11sd)*rgd[0] +u12sd *rgd[1]);
450
451 phi_s[idirac+1]-=akappa*(-conj(u12s)*ru[0] +conj(u11s)*ru[1]
452 +conj(u12sd)*rd[0] +u11sd *rd[1]);
453 //Dirac term
454 phi_s[idirac+1]-=gam*(-conj(u12s)*rgu[0] +conj(u11s)*rgu[1]
455 -conj(u12sd)*rgd[0] -u11sd *rgd[1]);
456 }
457 //Timelike terms
458 else{
459 const T dk4ms=dk4m[i]; const T dk4ps=dk4p[did];
460 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
461
462 phi_s[idirac]+= -dk4ms*(u11s*(ru[0]+rgu[0])
463 +u12s*(ru[1]+rgu[1]));
464 phi_s[idirac]+= -dk4ps*(conj(u11sd)*(rd[0]-rgd[0])
465 -u12sd *(rd[1]-rgd[1]));
466 phi[i+kvol*(0+idirac)]=phi_s[idirac+0];
467
468 phi_s[idirac+1]-= dk4ms*(-conj(u12s)*(ru[0]+rgu[0])
469 +conj(u11s)*(ru[1]+rgu[1]));
470 phi_s[idirac+1]-= +dk4ps*(conj(u12sd)*(rd[0]-rgd[0])
471 +u11sd *(rd[1]-rgd[1]));
472 phi[i+kvol*(1+idirac)]=phi_s[idirac+1];
473 }
474 }
475 }
476 }
477 }

References conj(), kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cuHdslashd_f()

void cuHdslashd_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
float akappa,
dim3 dimGrid,
dim3 dimBlock )

GPU calling wrapper for \(\Phi=M^\dagger r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
[in]dimGrid,dimBlockCUDA grid/block
Postcondition
Result written to phi

Definition at line 711 of file cumatrices.cu.

712 {
713 const char funcname[] = "Hdslashd_f";
714 int cuCpyStat=0;
715 for(unsigned short j=0;j<nc*ndirac;j++)
716 if((cuCpyStat=cudaMemcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex_f),cudaMemcpyDefault))){
717 fprintf(stderr,"Error %d in %s: Cuda failed to copy managed r into device Phi with code %d.\nExiting,,,\n\n",\
718 CPYERROR,funcname,cuCpyStat);
719 exit(cuCpyStat);
720 }
721 Kernels::cuHdslashd<<<dimGrid,dimBlock>>>(phi,r,ut[0],ut[1],iu,id,gamval,gamin,dk[0],dk[1],akappa);
722 return;
723}

References Complex_f, CPYERROR, Kernels::cuHdslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Dslash()

int Dslash ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
Complex_f jqq,
float akappa )

Evaluates \(\Phi=M r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 16 of file matrices.c.

17 {
18 const char funcname[] = "Dslash";
19 //Get the halos in order
20#if(nproc>1)
21 ZHalo_swap_all(r, 16);
22#endif
23
24 //Mass term
25 //Diquark Term (antihermitian)
26#ifdef USE_GPU
27 cuDslash(phi,r,ut,iu,id,gamval,gamin,dk,jqq,akappa,dimGrid,dimBlock);
28#else
29 for(unsigned short j=0;j<nc*ngorkov;j++)
30 memcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex));
31#pragma omp parallel for simd
32 for(unsigned int i=0;i<kvol;i++){
33 Complex ru[nc]; Complex rd[nc];
34 Complex rgu[nc]; Complex rgd[nc];
35 Complex phi_s[ngorkov*nc];
36 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
37 unsigned short igork = ((idirac>>1)+4)<<1;
38 unsigned int ind_d =4*ndirac+(idirac>>1);
39 Complex a_1=conj(jqq)*gamval[ind_d];
40 //We subtract a_2, hence the minus
41 Complex a_2=-jqq*gamval[ind_d];
42 ind_d=i+kvolHalo*(idirac); unsigned int ind_g=i+kvolHalo*(igork);
43 phi_s[idirac]=phi[ind_d]+a_1*r[ind_g];
44 phi_s[igork]=phi[ind_g]+a_2*r[ind_d];
45 ind_d+=kvolHalo; ind_g+=kvolHalo;
46 phi_s[idirac+1]=phi[ind_d]+a_1*r[ind_g];
47 phi_s[igork+1]=phi[ind_g]+a_2*r[ind_d];
48 }
49 Complex u11s; Complex u12s;
50 Complex u11sd; Complex u12sd;
51 unsigned int ind;
52 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
53#ifndef NO_SPACE
54 for(unsigned short mu = 0; mu <3; mu++){
55 ind = i+kvol*mu;
56 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
57 ind = i+kvolHalo*mu;
58 u11s=ut[0][ind]; u12s=ut[1][ind];
59 ind = did+kvolHalo*mu;
60 u11sd=ut[0][ind]; u12sd=ut[1][ind];
61 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
62 unsigned short idirac=igorkov&3;
63 unsigned short gind=mu*ndirac+idirac;
64 const Complex gam=gamval[gind];
65 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing in the dirac term.
66 unsigned short igork1 = (igorkov<4) ? gamin[gind] : gamin[gind]+4;
67 for(unsigned short c=0;c<nc;c++){
68 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
69 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
70 }
71 //Wilson + Dirac term in that order. Definitely easier
72 phi_s[igorkov*nc]+=-akappa*(u11s*ru[0]+ u12s*ru[1]+\
73 conj(u11sd)*rd[0]- u12sd*rd[1]);
74 //Dirac term
75 phi_s[igorkov*nc]+=gam*(u11s*rgu[0]+ u12s*rgu[1]-\
76 conj(u11sd)*rgd[0]+ u12sd*rgd[1]);
77
78 phi_s[igorkov*nc+1]+=-akappa*(-conj(u12s)*ru[0]+ conj(u11s)*ru[1]+\
79 conj(u12sd)*rd[0]+ u11sd*rd[1]);
80 //Dirac term
81 phi_s[igorkov*nc+1]+=gam*(-conj(u12s)*rgu[0]+ conj(u11s)*rgu[1]-\
82 conj(u12sd)*rgd[0]- u11sd*rgd[1]);
83 }
84 }
85 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
86 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
87 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
88#endif
89#ifndef NO_TIME
90 ind=i+kvolHalo*3;
91 u11s=ut[0][ind]; u12s=ut[1][ind];
92 const double dk4ms=dk[0][i]; const double dk4ps=dk[1][i];
93 ind=i+kvol*3;
94 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
95 ind=did+kvolHalo*3;
96 u11sd=ut[0][ind]; u12sd=ut[1][ind];
97 const double dk4msd=dk[0][did]; const double dk4psd=dk[1][did];
98 for(unsigned short igorkov=0;igorkov<ndirac;igorkov++){
99 unsigned short igork1 = gamin[3*ndirac+igorkov];
100 for(unsigned short c=0;c<nc;c++){
101 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
102 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
103 }
104 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
105 phi_s[igorkov*nc]+=
106 -dk4ps*(u11s*(ru[0]-rgu[0]) +u12s*(ru[1]-rgu[1]))
107 -dk4msd*(conj(u11sd)*(rd[0]+rgd[0]) -u12sd *(rd[1]+rgd[1]));
108 phi[i+kvolHalo*(igorkov*nc)]=phi_s[igorkov*nc];
109
110 phi_s[igorkov*nc+1]+=
111 -dk4ps*(-conj(u12s)*(ru[0]-rgu[0]) +conj(u11s)*(ru[1]-rgu[1]))
112 -dk4msd*(conj(u12sd)*(rd[0]+rgd[0]) +u11sd *(rd[1]+rgd[1]));
113 phi[i+kvolHalo*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
114 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
115 //the FORTRAN code did it.
116 igork1 += 4;
117 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
118 for(unsigned short c=0;c<nc;c++){
119 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
120 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
121 }
122 phi_s[igorkovPP*nc]+=-dk4ms*(u11s*(ru[0]-rgu[0])+ u12s*(ru[1]-rgu[1]))-
123 dk4psd*(conj(u11sd)*(rd[0]+rgd[0])- u12sd*(rd[1]+rgd[1]));
124 phi[i+kvolHalo*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
125
126 phi_s[igorkovPP*nc+1]+=-dk4ms*(conj(-u12s)*(ru[0]-rgu[0]) +conj(u11s)*(ru[1]-rgu[1]))
127 -dk4psd*(conj(u12sd)*(rd[0]+rgd[0]) +u11sd*(rd[1]+rgd[1]));
128 phi[i+kvolHalo*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
129 }
130#endif
131 }
132#endif
133 return 0;
134}
cuDslash
void cuDslash(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:609
ZHalo_swap_all
int ZHalo_swap_all(Complex *z, int ncpt)
Calls the functions to send data to both the up and down halos.
dimBlock
dim3 dimBlock
Default block size. Usually 128.
Definition cusu2hmc.cu:25

References Complex, Complex_f, conj(), cuDslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndirac, ngorkov, and ZHalo_swap_all().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Dslash_f()

int Dslash_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
Complex_f jqq,
float akappa )

Evaluates \(\Phi=M r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 425 of file matrices.c.

426 {
427 const char funcname[] = "Dslash_f";
428 //Get the halos in order
429#if(nproc>1)
430 CHalo_swap_all(r, 16);
431#endif
432
433 //Mass term
434 //Diquark Term (antihermitian)
435#ifdef USE_GPU
436 cuDslash_f(phi,r,ut,iu,id,gamval,gamin,dk,jqq,akappa,dimGrid,dimBlock);
437#else
438 for(unsigned short j=0;j<nc*ngorkov;j++)
439 memcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex_f));
440#pragma omp parallel for simd
441 for(unsigned int i=0;i<kvol;i++){
442 Complex_f ru[nc]; Complex_f rd[nc];
443 Complex_f rgu[nc]; Complex_f rgd[nc];
444 Complex_f phi_s[ngorkov*nc];
445 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
446 unsigned short igork = ((idirac>>1)+4)<<1;
447 unsigned int ind_d =4*ndirac+(idirac>>1);
448 Complex_f a_1=conjf(jqq)*gamval[ind_d];
449 //We subtract a_2, hence the minus
450 Complex_f a_2=-jqq*gamval[ind_d];
451 ind_d=i+kvolHalo*(idirac); unsigned int ind_g=i+kvolHalo*(igork);
452 phi_s[idirac]=phi[ind_d]+a_1*r[ind_g];
453 phi_s[igork]=phi[ind_g]+a_2*r[ind_d];
454 ind_d+=kvolHalo; ind_g+=kvolHalo;
455 phi_s[idirac+1]=phi[ind_d]+a_1*r[ind_g];
456 phi_s[igork+1]=phi[ind_g]+a_2*r[ind_d];
457 }
458 Complex_f u11s; Complex_f u12s;
459 Complex_f u11sd; Complex_f u12sd;
460 unsigned int ind;
461 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
462#ifndef NO_SPACE
463 for(unsigned short mu = 0; mu <3; mu++){
464 ind = i+kvol*mu;
465 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
466 ind = i+kvolHalo*mu;
467 u11s=ut[0][ind]; u12s=ut[1][ind];
468 ind = did+kvolHalo*mu;
469 u11sd=ut[0][ind]; u12sd=ut[1][ind];
470 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
471 unsigned short idirac=igorkov&3;
472 unsigned short gind=mu*ndirac+idirac;
473 const Complex_f gam=gamval[gind];
474 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing in the dirac term.
475 unsigned short igork1 = (igorkov<4) ? gamin[gind] : gamin[gind]+4;
476 for(unsigned short c=0;c<nc;c++){
477 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
478 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
479 }
480 //Wilson + Dirac term in that order. Definitely easier
481 phi_s[igorkov*nc]+=-akappa*(u11s*ru[0]+ u12s*ru[1]+\
482 conjf(u11sd)*rd[0]- u12sd*rd[1]);
483 //Dirac term
484 phi_s[igorkov*nc]+=gam*(u11s*rgu[0]+ u12s*rgu[1]-\
485 conjf(u11sd)*rgd[0]+ u12sd*rgd[1]);
486
487 phi_s[igorkov*nc+1]+=-akappa*(-conjf(u12s)*ru[0]+ conjf(u11s)*ru[1]+\
488 conjf(u12sd)*rd[0]+ u11sd*rd[1]);
489 //Dirac term
490 phi_s[igorkov*nc+1]+=gam*(-conjf(u12s)*rgu[0]+ conjf(u11s)*rgu[1]-\
491 conjf(u12sd)*rgd[0]- u11sd*rgd[1]);
492 }
493 }
494 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
495 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
496 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
497#endif
498#ifndef NO_TIME
499 ind=i+kvolHalo*3;
500 u11s=ut[0][ind]; u12s=ut[1][ind];
501 const float dk4ms=dk[0][i]; const float dk4ps=dk[1][i];
502 ind = i+kvol*3;
503 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
504 ind=did+kvolHalo*3;
505 u11sd=ut[0][ind]; u12sd=ut[1][ind];
506 const float dk4msd=dk[0][did]; const float dk4psd=dk[1][did];
507 for(unsigned short igorkov=0;igorkov<ndirac;igorkov++){
508 unsigned short igork1 = gamin[3*ndirac+igorkov];
509 for(unsigned short c=0;c<nc;c++){
510 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
511 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
512 }
513 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
514 phi_s[igorkov*nc]+=
515 -dk4ps*(u11s*(ru[0]-rgu[0]) +u12s*(ru[1]-rgu[1]))
516 -dk4msd*(conjf(u11sd)*(rd[0]+rgd[0]) -u12sd *(rd[1]+rgd[1]));
517 phi[i+kvolHalo*(igorkov*nc)]=phi_s[igorkov*nc];
518
519 phi_s[igorkov*nc+1]+=
520 -dk4ps*(-conjf(u12s)*(ru[0]-rgu[0]) +conjf(u11s)*(ru[1]-rgu[1]))
521 -dk4msd*(conjf(u12sd)*(rd[0]+rgd[0]) +u11sd *(rd[1]+rgd[1]));
522 phi[i+kvolHalo*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
523 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
524 //the FORTRAN code did it.
525 igork1 += 4;
526 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
527 for(unsigned short c=0;c<nc;c++){
528 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
529 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
530 }
531 phi_s[igorkovPP*nc]+=-dk4ms*(u11s*(ru[0]-rgu[0])+ u12s*(ru[1]-rgu[1]))-
532 dk4psd*(conjf(u11sd)*(rd[0]+rgd[0])- u12sd*(rd[1]+rgd[1]));
533 phi[i+kvolHalo*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
534
535 phi_s[igorkovPP*nc+1]+=-dk4ms*(conjf(-u12s)*(ru[0]-rgu[0]) +conjf(u11s)*(ru[1]-rgu[1]))
536 -dk4psd*(conjf(u12sd)*(rd[0]+rgd[0]) +u11sd*(rd[1]+rgd[1]));
537 phi[i+kvolHalo*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
538 }
539#endif
540 }
541#endif
542 return 0;
543}
cuDslash_f
void cuDslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:668
CHalo_swap_all
int CHalo_swap_all(Complex_f *c, int ncpt)
Calls the functions to send data to both the up and down halos.

References CHalo_swap_all(), Complex_f, cuDslash_f(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndirac, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Dslashd()

int Dslashd ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
Complex_f jqq,
float akappa )

Evaluates \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 135 of file matrices.c.

136 {
137 const char funcname[] = "Dslashd";
138 //Get the halos in order
139#if(nproc>1)
140 ZHalo_swap_all(r, 16);
141#endif
142
143 //Mass term
144#ifdef USE_GPU
145 cuDslashd(phi,r,ut,iu,id,gamval,gamin,dk,jqq,akappa,dimGrid,dimBlock);
146#else
147 for(unsigned short j=0;j<nc*ngorkov;j++)
148 memcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex));
149#pragma omp parallel for simd
150 for(unsigned int i=0;i<kvol;i++){
151 Complex ru[nc]; Complex rd[nc];
152 Complex rgu[nc]; Complex rgd[nc];
153 Complex phi_s[ngorkov*nc];
154 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
155 unsigned short igork = ((idirac>>1)+4)<<1;
156 unsigned int ind_d =4*ndirac+(idirac>>1);
157 Complex a_1=-conj(jqq)*gamval[ind_d];
158 Complex a_2=jqq*gamval[ind_d];
159 //ind_d=i+kvolHalo*(idirac); unsigned int ind_g=i+kvolHalo*(igork);
160 phi_s[idirac]=phi[i+kvol*idirac]+a_1*r[i+kvolHalo*igork];
161 phi_s[igork]=phi[i+kvol*igork]+a_2*r[i+kvolHalo*idirac];
162 //ind_d+=kvolHalo; ind_g+=kvolHalo;
163 phi_s[idirac+1]=phi[i+kvol*(idirac+1)]+a_1*r[i+kvolHalo*(igork+1)];
164 phi_s[igork+1]=phi[i+kvol*(igork+1)]+a_2*r[i+kvolHalo*(idirac+1)];
165 }
166 Complex u11s; Complex u12s;
167 Complex u11sd; Complex u12sd;
168 unsigned int ind;
169 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
170#ifndef NO_SPACE
171 for(unsigned short mu = 0; mu <3; mu++){
172 ind = i+kvol*mu;
173 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
174 ind = i+kvolHalo*mu;
175 u11s=ut[0][ind]; u12s=ut[1][ind];
176 ind = did+kvolHalo*mu;
177 u11sd=ut[0][ind]; u12sd=ut[1][ind];
178 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
179 unsigned short idirac=igorkov&3;
180 const Complex gam=gamval[mu*ndirac+idirac];
181 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing.
182 unsigned short igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
183 for(unsigned short c=0;c<nc;c++){
184 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
185 rgd[c]=r[did+kvolHalo*(igork1*nc+c)]; rgu[c]=r[uid+kvolHalo*(igork1*nc+c)];
186 }
187 //Wilson + Dirac term in that order. Definitely easier
188 phi_s[igorkov*nc]-= akappa*(u11s*ru[0] +u12s*ru[1]
189 +conj(u11sd)*rd[0] -u12sd *rd[1]);
190
191 //Dirac term
192 phi_s[igorkov*nc]-=gam* (u11s*rgu[0] +u12s*rgu[1]
193 -conj(u11sd)*rgd[0] +u12sd *rgd[1]);
194
195 phi_s[igorkov*nc+1]-= akappa*(-conj(u12s)*ru[0] +conj(u11s)*ru[1]
196 +conj(u12sd)*rd[0] +u11sd *rd[1]);
197 //Dirac term
198 phi_s[igorkov*nc+1]-=gam* (-conj(u12s)*rgu[0] +conj(u11s)*rgu[1]
199 -conj(u12sd)*rgd[0] -u11sd *rgd[1]);
200
201 }
202 }
203#endif
204 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
205 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
206 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
207 //Under dagger, dk4p and dk4m get swapped and the dirac component flips sign.
208#ifndef NO_TIME
209 ind=i+kvolHalo*3;
210 u11s=ut[0][ind]; u12s=ut[1][ind];
211 const double dk4ms=dk[0][i]; const double dk4ps=dk[1][i];
212 ind = i+kvol*3;
213 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
214 ind=did+kvolHalo*3;
215 u11sd=ut[0][ind]; u12sd=ut[1][ind];
216 const double dk4msd=dk[0][did]; const double dk4psd=dk[1][did];
217 for(unsigned short igorkov=0; igorkov<ndirac; igorkov++){
218 unsigned short igork1 = gamin[3*ndirac+igorkov];
219 for(unsigned short c=0;c<nc;c++){
220 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
221 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
222 }
223 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
224 phi_s[igorkov*nc]+=
225 -dk4ms*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
226 -dk4psd*(conj(u11sd)*(rd[0]-rgd[0]) -u12sd *(rd[1]-rgd[1]));
227 phi[i+kvol*(igorkov*nc)]=phi_s[igorkov*nc];
228
229 phi_s[igorkov*nc+1]+=
230 -dk4ms*(-conj(u12s)*(ru[0]+rgu[0]) +conj(u11s)*(ru[1]+rgu[1]))
231 -dk4psd*(conj(u12sd)*(rd[0]-rgd[0]) +u11sd *(rd[1]-rgd[1]));
232 phi[i+kvol*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
233 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
234 //the FORTRAN code did it.
235 igork1 += 4;
236 for(unsigned short c=0;c<nc;c++){
237 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
238 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
239 }
240 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
241 phi_s[igorkovPP*nc]+=-dk4ps*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
242 -dk4msd*(conj(u11sd)*(rd[0]-rgd[0]) -u12sd*(rd[1]-rgd[1]));
243 phi[i+kvol*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
244
245 phi_s[igorkovPP*nc+1]+=dk4ps*(conj(u12s)*(ru[0]+rgu[0]) -conj(u11s)*(ru[1]+rgu[1]))
246 -dk4msd*(conj(u12sd)*(rd[0]-rgd[0]) +u11sd*(rd[1]-rgd[1]));
247 phi[i+kvol*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
248 }
249#endif
250 }
251#endif
252 return 0;
253}
cuDslashd
void cuDslashd(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:623

References Complex, Complex_f, conj(), cuDslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndirac, ngorkov, and ZHalo_swap_all().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Dslashd_f()

int Dslashd_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
Complex_f jqq,
float akappa )

Evaluates \(\Phi=M^\dagger r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]jqqDiquark source
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 544 of file matrices.c.

545 {
546 const char funcname[] = "Dslashd_f";
547 //Get the halos in order
548#if(nproc>1)
549 CHalo_swap_all(r, 16);
550#endif
551
552 //Mass term
553#ifdef USE_GPU
554 cuDslashd_f(phi,r,ut,iu,id,gamval,gamin,dk,jqq,akappa,dimGrid,dimBlock);
555#else
556 for(unsigned short j=0;j<nc*ngorkov;j++)
557 memcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex_f));
558#pragma omp parallel for simd
559 for(unsigned int i=0;i<kvol;i++){
560 Complex_f ru[nc]; Complex_f rd[nc];
561 Complex_f rgu[nc]; Complex_f rgd[nc];
562 Complex_f phi_s[ngorkov*nc];
563 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
564 unsigned short igork = ((idirac>>1)+4)<<1;
565 unsigned int ind_d =4*ndirac+(idirac>>1);
566 Complex_f a_1=-conjf(jqq)*gamval[ind_d];
567 Complex_f a_2=jqq*gamval[ind_d];
568 //ind_d=i+kvolHalo*(idirac); unsigned int ind_g=i+kvolHalo*(igork);
569 phi_s[idirac]=phi[i+kvol*idirac]+a_1*r[i+kvolHalo*igork];
570 phi_s[igork]=phi[i+kvol*igork]+a_2*r[i+kvolHalo*idirac];
571 //ind_d+=kvolHalo; ind_g+=kvolHalo;
572 phi_s[idirac+1]=phi[i+kvol*(idirac+1)]+a_1*r[i+kvolHalo*(igork+1)];
573 phi_s[igork+1]=phi[i+kvol*(igork+1)]+a_2*r[i+kvolHalo*(idirac+1)];
574 }
575 Complex_f u11s; Complex_f u12s;
576 Complex_f u11sd; Complex_f u12sd;
577 unsigned int ind;
578 //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
579#ifndef NO_SPACE
580 for(unsigned short mu = 0; mu <3; mu++){
581 ind = i+kvol*mu;
582 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
583 ind = i+kvolHalo*mu;
584 u11s=ut[0][ind]; u12s=ut[1][ind];
585 ind = did+kvolHalo*mu;
586 u11sd=ut[0][ind]; u12sd=ut[1][ind];
587 for(unsigned short igorkov=0; igorkov<ngorkov; igorkov++){
588 unsigned short idirac=igorkov&3;
589 const Complex_f gam=gamval[mu*ndirac+idirac];
590 //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing.
591 unsigned short igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
592 for(unsigned short c=0;c<nc;c++){
593 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
594 rgd[c]=r[did+kvolHalo*(igork1*nc+c)]; rgu[c]=r[uid+kvolHalo*(igork1*nc+c)];
595 }
596 //Wilson + Dirac term in that order. Definitely easier
597 phi_s[igorkov*nc]-= akappa*(u11s*ru[0] +u12s*ru[1]
598 +conjf(u11sd)*rd[0] -u12sd *rd[1]);
599
600 //Dirac term
601 phi_s[igorkov*nc]-=gam* (u11s*rgu[0] +u12s*rgu[1]
602 -conjf(u11sd)*rgd[0] +u12sd *rgd[1]);
603
604 phi_s[igorkov*nc+1]-= akappa*(-conjf(u12s)*ru[0] +conjf(u11s)*ru[1]
605 +conjf(u12sd)*rd[0] +u11sd *rd[1]);
606 //Dirac term
607 phi_s[igorkov*nc+1]-=gam* (-conjf(u12s)*rgu[0] +conjf(u11s)*rgu[1]
608 -conjf(u12sd)*rgd[0] -u11sd *rgd[1]);
609
610 }
611 }
612#endif
613 //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
614 //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7 separately
615 //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
616 //Under dagger, dk4p and dk4m get swapped and the dirac component flips sign.
617#ifndef NO_TIME
618 ind=i+kvolHalo*3;
619 u11s=ut[0][ind]; u12s=ut[1][ind];
620 const float dk4ms=dk[0][i]; const float dk4ps=dk[1][i];
621 ind = i+kvol*3;
622 const unsigned int did=id[ind]; const unsigned int uid = iu[ind];
623 ind=did+kvolHalo*3;
624 u11sd=ut[0][ind]; u12sd=ut[1][ind];
625 const float dk4msd=dk[0][did]; const float dk4psd=dk[1][did];
626 for(unsigned short igorkov=0; igorkov<ndirac; igorkov++){
627 unsigned short igork1 = gamin[3*ndirac+igorkov];
628 for(unsigned short c=0;c<nc;c++){
629 ru[c]=r[uid+kvolHalo*(igorkov*nc+c)]; rd[c]=r[did+kvolHalo*(igorkov*nc+c)];
630 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
631 }
632 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
633 phi_s[igorkov*nc]+=
634 -dk4ms*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
635 -dk4psd*(conjf(u11sd)*(rd[0]-rgd[0]) -u12sd *(rd[1]-rgd[1]));
636 phi[i+kvol*(igorkov*nc)]=phi_s[igorkov*nc];
637
638 phi_s[igorkov*nc+1]+=
639 -dk4ms*(-conjf(u12s)*(ru[0]+rgu[0]) +conjf(u11s)*(ru[1]+rgu[1]))
640 -dk4psd*(conjf(u12sd)*(rd[0]-rgd[0]) +u11sd *(rd[1]-rgd[1]));
641 phi[i+kvol*(igorkov*nc+1)]=phi_s[igorkov*nc+1];
642 const unsigned short igorkovPP=igorkov+4; //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
643 //the FORTRAN code did it.
644 igork1 += 4;
645 for(unsigned short c=0;c<nc;c++){
646 ru[c]=r[uid+kvolHalo*(igorkovPP*nc+c)]; rd[c]=r[did+kvolHalo*(igorkovPP*nc+c)];
647 rgu[c]=r[uid+kvolHalo*(igork1*nc+c)]; rgd[c]=r[did+kvolHalo*(igork1*nc+c)];
648 }
649 //And the Gor'kov terms. Note that dk4p and dk4m swap positions compared to the above
650 phi_s[igorkovPP*nc]+=-dk4ps*(u11s*(ru[0]+rgu[0]) +u12s*(ru[1]+rgu[1]))
651 -dk4msd*(conjf(u11sd)*(rd[0]-rgd[0]) -u12sd*(rd[1]-rgd[1]));
652 phi[i+kvol*(igorkovPP*nc)]=phi_s[igorkovPP*nc];
653
654 phi_s[igorkovPP*nc+1]+=dk4ps*(conjf(u12s)*(ru[0]+rgu[0]) -conjf(u11s)*(ru[1]+rgu[1]))
655 -dk4msd*(conjf(u12sd)*(rd[0]-rgd[0]) +u11sd*(rd[1]-rgd[1]));
656 phi[i+kvol*(igorkovPP*nc+1)]=phi_s[igorkovPP*nc+1];
657 }
658#endif
659 }
660#endif
661 return 0;
662}
cuDslashd_f
void cuDslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:682

References CHalo_swap_all(), Complex_f, cuDslashd_f(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndirac, and ngorkov.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Hdslash()

int Hdslash ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
float akappa )

Evaluates \(\Phi=M r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge trial field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\)
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 254 of file matrices.c.

255 {
256 const char funcname[] = "Hdslash";
257 //Get the halos in order
258#if(nproc>1)
259 ZHalo_swap_all(r, 8);
260#endif
261
262 //Mass term
263 //Spacelike term
264#ifdef USE_GPU
265 cuHdslash(phi,r,ut,iu,id,gamval,gamin,dk,akappa,dimGrid,dimBlock);
266#else
267 for(unsigned short j=0;j<nc*ndirac;j++)
268 memcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex));
269#pragma omp parallel for simd
270 for(unsigned int i=0;i<kvol;i++){
271 Complex ru[nc]; Complex rd[nc];
272 Complex rgu[nc]; Complex rgd[nc];
273 Complex phi_s[ndirac*nc];
274 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
275#pragma unroll
276 for(unsigned short c=0; c<nc; c++)
277 //NOTE: idirac is increasing by nc each time. So should be read as idirac*nc in a Dirac-counted loop
278 phi_s[idirac+c]=phi[i+kvolHalo*(c+idirac)];
279
280 //#pragma unroll
281 for(unsigned short mu = 0; mu <ndim; mu++){
282 unsigned int ind=i+kvolHalo*mu;
283 const Complex u11s=ut[0][ind]; const Complex u12s=ut[1][ind];
284 ind = i+kvol*mu;
285 const int did=id[ind]; const int uid = iu[ind];
286 ind=did+kvolHalo*mu;
287 const Complex u11sd=ut[0][ind]; const Complex u12sd=ut[1][ind];
288 for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
289 const unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
290#pragma unroll
291 for(unsigned short c=0;c<nc;c++){
292 ind =kvolHalo*(idirac+c);
293 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
294 ind =kvolHalo*(igork1+c);
295 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
296 }
297 //Can manually vectorise with a pragma?
298 //Wilson + Dirac term in that order. Definitely easier
299 //to read when split into different loops, but should be faster this way
300 //Spacelike terms
301 if(mu<3){
302 const Complex gam=gamval[mu*ndirac+(idirac>>1)];
303 phi_s[idirac]+=-akappa*(u11s*ru[0]+u12s*ru[1]+\
304 conj(u11sd)*rd[0]-u12sd*rd[1]);
305 //Dirac term
306 phi_s[idirac]+=gam*(u11s*rgu[0]+u12s*rgu[1]-\
307 conj(u11sd)*rgd[0]+ u12sd*rgd[1]);
308
309 phi_s[idirac+1]+=-akappa*(-conj(u12s)*ru[0]+ conj(u11s)*ru[1]+\
310 conj(u12sd)*rd[0]+ u11sd*rd[1]);
311 //Dirac term
312 phi_s[idirac+1]+=gam*(-conj(u12s)*rgu[0]+ conj(u11s)*rgu[1]-\
313 conj(u12sd)*rgd[0]- u11sd*rgd[1]);
314 }
315 //Timelike terms
316 else{
317 const double dk4ms=dk[0][did]; const double dk4ps=dk[1][i];
318 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
319
320 phi_s[idirac+0]-= dk4ps*(u11s*(ru[0]-rgu[0])
321 +u12s*(ru[1]-rgu[1]));
322 phi_s[idirac+0]-= dk4ms*(conj(u11sd)*(rd[0]+rgd[0])
323 -u12sd *(rd[1]+rgd[1]));
324 phi[i+kvolHalo*(0+idirac)]=phi_s[idirac+0];
325
326 phi_s[idirac+1]-= dk4ps*(-conj(u12s)*(ru[0]-rgu[0])
327 +conj(u11s)*(ru[1]-rgu[1]));
328 phi_s[idirac+1]-= dk4ms*(conj(u12sd)*(rd[0]+rgd[0])
329 +u11sd *(rd[1]+rgd[1]));
330 phi[i+kvolHalo*(1+idirac)]=phi_s[idirac+1];
331 }
332 }
333 }
334 }
335#endif
336 return 0;
337}
cuHdslash
void cuHdslash(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:637

References Complex, conj(), cuHdslash(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndim, ndirac, and ZHalo_swap_all().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Hdslash_f()

int Hdslash_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
float akappa )

Evaluates \(\Phi=M r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 663 of file matrices.c.

664 {
665 const char funcname[] = "Hdslash_f";
666 //Get the halos in order
667#if(nproc>1)
668 CHalo_swap_all(r, 8);
669#endif
670#ifdef USE_GPU
671 cuHdslash_f(phi,r,ut,iu,id,gamval,gamin,dk,akappa,dimGrid,dimBlock);
672#else
673 //Mass term
674 for(unsigned short j=0;j<nc*ndirac;j++)
675 memcpy(phi+j*kvolHalo, r+j*kvolHalo, kvol*sizeof(Complex_f));
676#pragma omp parallel for simd
677 for(unsigned int i=0;i<kvol;i++){
678 Complex_f ru[nc]; Complex_f rd[nc];
679 Complex_f rgu[nc]; Complex_f rgd[nc];
680 Complex_f phi_s[ndirac*nc];
681 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
682#pragma unroll
683 for(unsigned short c=0; c<nc; c++)
684 //NOTE: idirac is increasing by nc each time.
685 //So should be read as idirac*nc in a Dirac-counted loop
686 phi_s[idirac+c]=phi[i+kvolHalo*(c+idirac)];
687
688 //#pragma unroll
689 for(unsigned short mu = 0; mu <ndim; mu++){
690 unsigned int ind=i+kvolHalo*mu;
691 const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
692 ind = i+kvol*mu;
693 const int did=id[ind]; const int uid = iu[ind];
694 ind=did+kvolHalo*mu;
695 const Complex_f u11sd=ut[0][ind]; const Complex_f u12sd=ut[1][ind];
696 for(unsigned short idirac=0; idirac<ndirac*nc; idirac+=nc){
697 const unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
698#pragma unroll
699 for(unsigned short c=0;c<nc;c++){
700 ind =kvolHalo*(idirac+c);
701 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
702 ind =kvolHalo*(igork1+c);
703 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
704 }
705 //Can manually vectorise with a pragma?
706 //Wilson + Dirac term in that order. Definitely easier
707 //to read when split into different loops, but should be faster this way
708 //Spacelike terms
709 if(mu<3){
710 const Complex_f gam=gamval[mu*ndirac+(idirac>>1)];
711 phi_s[idirac]+=-akappa*(u11s*ru[0]+u12s*ru[1]+\
712 conjf(u11sd)*rd[0]-u12sd*rd[1]);
713 //Dirac term
714 phi_s[idirac]+=gam*(u11s*rgu[0]+u12s*rgu[1]-\
715 conjf(u11sd)*rgd[0]+ u12sd*rgd[1]);
716
717 phi_s[idirac+1]+=-akappa*(-conjf(u12s)*ru[0]+ conjf(u11s)*ru[1]+\
718 conjf(u12sd)*rd[0]+ u11sd*rd[1]);
719 //Dirac term
720 phi_s[idirac+1]+=gam*(-conjf(u12s)*rgu[0]+ conjf(u11s)*rgu[1]-\
721 conjf(u12sd)*rgd[0]- u11sd*rgd[1]);
722 }
723 //Timelike terms
724 else{
725 const float dk4ms=dk[0][did]; const float dk4ps=dk[1][i];
726 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
727
728 phi_s[idirac+0]-= dk4ps*(u11s*(ru[0]-rgu[0])
729 +u12s*(ru[1]-rgu[1]));
730 phi_s[idirac+0]-= dk4ms*(conjf(u11sd)*(rd[0]+rgd[0])
731 -u12sd *(rd[1]+rgd[1]));
732 phi[i+kvolHalo*(0+idirac)]=phi_s[idirac+0];
733
734 phi_s[idirac+1]-= dk4ps*(-conjf(u12s)*(ru[0]-rgu[0])
735 +conjf(u11s)*(ru[1]-rgu[1]));
736 phi_s[idirac+1]-= dk4ms*(conjf(u12sd)*(rd[0]+rgd[0])
737 +u11sd *(rd[1]+rgd[1]));
738 phi[i+kvolHalo*(1+idirac)]=phi_s[idirac+1];
739 }
740 }
741 }
742 }
743#endif
744 return 0;
745}
cuHdslash_f
void cuHdslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in single precision.
Definition cumatrices.cu:696

References CHalo_swap_all(), Complex_f, cuHdslash_f(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Hdslashd()

int Hdslashd ( Complex * phi,
Complex * r,
Complex * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex gamval[20],
const unsigned short gamin[16],
double * dk[nc],
float akappa )

Evaluates \(\Phi=M^\dagger r\) in double precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 338 of file matrices.c.

339 {
340 const char funcname[] = "Hdslashd";
341 //Get the halos in order.
342#if(nproc>1)
343 ZHalo_swap_all(r, 8);
344#endif
345
346 //Mass term
347#ifdef USE_GPU
348 cuHdslashd(phi,r,ut,iu,id,gamval,gamin,dk,akappa,dimGrid,dimBlock);
349#else
350 for(unsigned short j=0;j<nc*ndirac;j++)
351 memcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex));
352 //Spacelike term
353#pragma omp parallel for simd
354 for(unsigned int i=0;i<kvol;i++){
355 //Right. Time to prefetch
356 Complex ru[nc]; Complex rd[nc];
357 Complex rgu[nc]; Complex rgd[nc];
358 Complex phi_s[ndirac*nc];
359 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
360#pragma unroll
361 for(unsigned short c=0; c<nc; c++)
362 //NOTE: idirac is increasing by nc each time. So should be read as idirac*nc in a Dirac-counted loop
363 phi_s[idirac+c]=phi[i+kvol*(c+idirac)];
364
365 //#pragma unroll
366 for(unsigned short mu = 0; mu <ndim; mu++){
367 unsigned int ind=i+kvolHalo*mu;
368 const Complex u11s=ut[0][ind]; const Complex u12s=ut[1][ind];
369 ind = i+kvol*mu;
370 const int did=id[ind]; const int uid = iu[ind];
371 ind=did+kvolHalo*mu;
372 const Complex u11sd=ut[0][ind]; const Complex u12sd=ut[1][ind];
373 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc){
374 unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
375#pragma unroll
376 for(unsigned short c=0;c<nc;c++){
377 ind =kvolHalo*(idirac+c);
378 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
379 ind =kvolHalo*(igork1+c);
380 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
381 }
382 //Can manually vectorise with a pragma?
383 //Wilson + Dirac term in that order. Definitely easier
384 //to read when split into different loops, but should be faster this way
385 //Spacelike terms
386 if(mu<3){
387 const Complex gam=gamval[mu*ndirac+(idirac>>1)];
388 phi_s[idirac]-=akappa*(u11s*ru[0] +u12s*ru[1]
389 +conj(u11sd)*rd[0] -u12sd *rd[1]);
390 //Dirac term
391 phi_s[idirac]-=gam* (u11s*rgu[0] +u12s*rgu[1]
392 -conj(u11sd)*rgd[0] +u12sd *rgd[1]);
393
394 phi_s[idirac+1]-=akappa*(-conj(u12s)*ru[0] +conj(u11s)*ru[1]
395 +conj(u12sd)*rd[0] +u11sd *rd[1]);
396 //Dirac term
397 phi_s[idirac+1]-=gam*(-conj(u12s)*rgu[0] +conj(u11s)*rgu[1]
398 -conj(u12sd)*rgd[0] -u11sd *rgd[1]);
399 }
400 //Timelike terms
401 else{
402 const double dk4ms=dk[0][i]; const double dk4ps=dk[1][did];
403 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
404
405 phi_s[idirac]+= -dk4ms*(u11s*(ru[0]+rgu[0])
406 +u12s*(ru[1]+rgu[1]));
407 phi_s[idirac]+= -dk4ps*(conj(u11sd)*(rd[0]-rgd[0])
408 -u12sd *(rd[1]-rgd[1]));
409 phi[i+kvol*(0+idirac)]=phi_s[idirac+0];
410
411 phi_s[idirac+1]-= dk4ms*(-conj(u12s)*(ru[0]+rgu[0])
412 +conj(u11s)*(ru[1]+rgu[1]));
413 phi_s[idirac+1]-= +dk4ps*(conj(u12sd)*(rd[0]-rgd[0])
414 +u11sd *(rd[1]-rgd[1]));
415 phi[i+kvol*(1+idirac)]=phi_s[idirac+1];
416 }
417 }
418 }
419 }
420#endif
421 return 0;
422}
cuHdslashd
void cuHdslashd(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in double precision.
Definition cumatrices.cu:651

References Complex, conj(), cuHdslashd(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndim, ndirac, and ZHalo_swap_all().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Hdslashd_f()

int Hdslashd_f ( Complex_f * phi,
Complex_f * r,
Complex_f * ut[nc],
unsigned int * iu,
unsigned int * id,
Complex_f gamval[20],
const unsigned short gamin[16],
float * dk[nc],
float akappa )

Evaluates \(\Phi=M^\dagger r\) in single precision.

Parameters
[out]phiThe product
[in]rThe array being acted on by M
[in]utGauge field
[in]iu,idUpper/lower halo indices
[in]gamvalGamma matrices rescaled by kappa
[in]gaminIndices for dirac terms
[in]dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1+\gamma_0\right)e^{+\mu}\)
[in]akappaHopping parameter
Postcondition
Result written to phi
Returns
Zero on success, integer error code otherwise

Definition at line 746 of file matrices.c.

747 {
748 const char funcname[] = "Hdslashd_f";
749 //Get the halos in order. Because C is row major, we need to extract the correct
750 //terms for each halo first. Changing the indices was considered but that caused
751 //issues with the BLAS routines.
752#if(nproc>1)
753 CHalo_swap_all(r, 8);
754#endif
755
756 //Mass term
757#ifdef USE_GPU
758 cuHdslashd_f(phi,r,ut,iu,id,gamval,gamin,dk,akappa,dimGrid,dimBlock);
759#else
760 for(unsigned short j=0;j<nc*ndirac;j++)
761 memcpy(phi+j*kvol, r+j*kvolHalo, kvol*sizeof(Complex_f));
762
763 //Spacelike term
764#pragma omp parallel for simd
765 for(unsigned int i=0;i<kvol;i++){
766 //Right. Time to prefetch
767 Complex_f ru[nc]; Complex_f rd[nc];
768 Complex_f rgu[nc]; Complex_f rgd[nc];
769 Complex_f phi_s[ndirac*nc];
770 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc)
771#pragma unroll
772 for(unsigned short c=0; c<nc; c++)
773 //NOTE: idirac is increasing by nc each time. So should be read as idirac*nc in a Dirac counted loop
774 phi_s[idirac+c]=phi[i+kvol*(c+idirac)];
775
776 //#pragma unroll
777 for(unsigned short mu = 0; mu <ndim; mu++){
778 unsigned int ind=i+kvolHalo*mu;
779 const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
780 ind = i+kvol*mu;
781 const int did=id[ind]; const int uid = iu[ind];
782 ind=did+kvolHalo*mu;
783 const Complex_f u11sd=ut[0][ind]; const Complex_f u12sd=ut[1][ind];
784 for(unsigned short idirac=0; idirac<nc*ndirac; idirac+=nc){
785 unsigned short igork1 = gamin[mu*ndirac+(idirac>>1)] << (nc-1);
786#pragma unroll
787 for(unsigned short c=0;c<nc;c++){
788 ind =kvolHalo*(idirac+c);
789 ru[c]=r[uid+ind]; rd[c]=r[did+ind];
790 ind =kvolHalo*(igork1+c);
791 rgu[c]=r[uid+ind]; rgd[c]=r[did+ind];
792 }
793 //Can manually vectorise with a pragma?
794 //Wilson + Dirac term in that order. Definitely easier
795 //to read when split into different loops, but should be faster this way
796 //Spacelike terms
797 if(mu<3){
798 const Complex_f gam=gamval[mu*ndirac+(idirac>>1)];
799 phi_s[idirac]-=akappa*(u11s*ru[0] +u12s*ru[1]
800 +conjf(u11sd)*rd[0] -u12sd *rd[1]);
801 //Dirac term
802 phi_s[idirac]-=gam* (u11s*rgu[0] +u12s*rgu[1]
803 -conjf(u11sd)*rgd[0] +u12sd *rgd[1]);
804
805 phi_s[idirac+1]-=akappa*(-conjf(u12s)*ru[0] +conjf(u11s)*ru[1]
806 +conjf(u12sd)*rd[0] +u11sd *rd[1]);
807 //Dirac term
808 phi_s[idirac+1]-=gam*(-conjf(u12s)*rgu[0] +conjf(u11s)*rgu[1]
809 -conjf(u12sd)*rgd[0] -u11sd *rgd[1]);
810 }
811 //Timelike terms
812 else{
813 const float dk4ms=dk[0][i]; const float dk4ps=dk[1][did];
814 //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
815
816 phi_s[idirac]+= -dk4ms*(u11s*(ru[0]+rgu[0])
817 +u12s*(ru[1]+rgu[1]));
818 phi_s[idirac]+= -dk4ps*(conjf(u11sd)*(rd[0]-rgd[0])
819 -u12sd *(rd[1]-rgd[1]));
820 phi[i+kvol*(0+idirac)]=phi_s[idirac+0];
821
822 phi_s[idirac+1]-= dk4ms*(-conjf(u12s)*(ru[0]+rgu[0])
823 +conjf(u11s)*(ru[1]+rgu[1]));
824 phi_s[idirac+1]-= +dk4ps*(conjf(u12sd)*(rd[0]-rgd[0])
825 +u11sd *(rd[1]-rgd[1]));
826 phi[i+kvol*(1+idirac)]=phi_s[idirac+1];
827 }
828 }
829 }
830 }
831#endif
832 return 0;
833}
cuHdslashd_f
void cuHdslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock)
GPU calling wrapper for in single precision.
Definition cumatrices.cu:711

References CHalo_swap_all(), Complex_f, cuHdslashd_f(), dimBlock, dimGrid, kvol, kvolHalo, nc, ndim, and ndirac.

Here is the call graph for this function:
Here is the caller graph for this function: