su2hmc
Loading...
Searching...
No Matches
cuforce.cu
Go to the documentation of this file.
1
7#include <matrices.h>
8#include <su2hmc.h>
9//CUDA Kernels
10namespace Kernels{
23 __global__ void Plus_staple(const int mu, const int nu,unsigned int *iu, Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11t, Complex_f *u12t){
24 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
25 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
26 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
27 const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
28 for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
29 const unsigned int uidm = iu[mu*kvol+i];
30 unsigned int indn=uidm+kvolHalo*nu;
31 const unsigned int uidn = iu[nu*kvol+i];
32 unsigned int indm=uidn+kvolHalo*mu;
33 Complex_f a11=u11t[indn]*conj(u11t[indm])+\
34 u12t[indn]*conj(u12t[indm]);
35 Complex_f a12=-u11t[indn]*u12t[indm]+\
36 u12t[indn]*u11t[indm];
37 indn=i+kvolHalo*nu;
38 Sigma11[i]+=a11*conj(u11t[indn])+a12*conj(u12t[indn]);
39 Sigma12[i]+=-a11*u12t[indn]+a12*u11t[indn];
40 }
41 }
42
55 __global__ void Minus_staple(const int mu,const int nu,unsigned int *iu,unsigned int *id, Complex_f *Sigma11, Complex_f *Sigma12,\
56 Complex_f *u11sh, Complex_f *u12sh, Complex_f *u11t, Complex_f *u12t){
57 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
58 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
59 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
60 const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
61 for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
62 const unsigned int uidm = iu[mu*kvol+i];
63 const unsigned int didn = id[nu*kvol+i];
64 //uidm is correct here
65 unsigned int ind=didn+kvolHalo*mu;
66 Complex_f u11s=u11t[ind]; Complex_f u12s=u12t[ind];
67 Complex_f a11=conj(u11sh[uidm])*conj(u11s)-\
68 u12sh[uidm]*conj(u12s);
69 Complex_f a12=-conj(u11sh[uidm])*u12s-\
70 u12sh[uidm]*u11s;
71 ind=didn+kvolHalo*nu;
72 u11s=u11t[ind]; u12s=u12t[ind];
73 Sigma11[i]+=a11*u11s-a12*conj(u12s);
74 Sigma12[i]+=a11*u12s+a12*conj(u11s);
75 }
76 }
77
89 __global__ void cuGaugeForce(int mu, Complex_f *Sigma11, Complex_f *Sigma12,double* dSdpi,Complex_f *u11t, Complex_f *u12t, float beta){
90 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
91 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
92 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
93 const unsigned int threadId= blockId * bsize+(threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
94 for(unsigned int i=threadId;i<kvol;i+=gsize*bsize){
95 const unsigned int ind = i+kvolHalo*mu;
96 Complex_f a11 = u11t[ind]*Sigma12[i]+u12t[ind]*conj(Sigma11[i]);
97 Complex_f a12 = u11t[ind]*Sigma11[i]+conj(u12t[ind])*Sigma12[i];
98 //Not worth splitting into different streams, before we get ideas...
99 dSdpi[i+kvol*mu]=beta*a11.imag();
100 dSdpi[i+kvol*(1*ndim+mu)]=beta*a11.real();
101 dSdpi[i+kvol*(2*ndim+mu)]=beta*a12.imag();
102 }
103 }
104
115 template <typename T>
116 __global__ void Gather(T *x, T *y, const unsigned int n, unsigned int *table, const unsigned short mu)
117 {
118 //FORTRAN had a second parameter m giving the size of y (kvol+halo) normally
119 //Pointers mean that's not an issue for us so I'm leaving it out
120 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
121 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
122 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
123 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
124 const unsigned int gthreadId= blockId * bsize+bthreadId;
125 const unsigned int kvbmu=kvolHalo*mu;
126 for(unsigned int i = gthreadId; i<kvol;i+=gsize*bsize)
127 x[i]=y[table[i+kvbmu]+kvbmu];
128 }
129
146 __global__ void cuForce_s(double *dSdpi, Complex_f *u11t, Complex_f *u12t, Complex_f *X1, Complex_f *X2, Complex_f gamval[20],\
147 unsigned int *iu, const unsigned short gamin[16],float akappa, const unsigned short mu){
148 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
149 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
150 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
151 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
152 const unsigned int gthreadId= blockId * bsize+bthreadId;
153
154 for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
155 const unsigned int ind=i+kvolHalo*mu;
156 const Complex_f u11s=u11t[ind]; const Complex_f u12s=u12t[ind];
157 const unsigned int uid = iu[i+kvol*mu];
158 //Similarly to Hdslash we always see idirac*nc so we do that here too.
159 for(unsigned short idirac=0;idirac<nc*ndirac;idirac+=nc){
160 Complex_f X1s[nc]; Complex_f X1su[nc];
161 Complex_f X2s[nc]; Complex_f X2su[nc];
162
163 X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
164 X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
165 X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
166 X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
167
168 // Need to be double to avoid accumulation errors
169 double dSdpis[3];
170 //Careful!! cant use ind here as dSdpi has no halo!
171 dSdpis[0]=dSdpi[i+kvol*mu];
172 //Multiplying by i and taking the real component is the same as taking the negative imaginary component
173 //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
174 //generators by hand
175 dSdpis[0]+=-akappa*(
176 conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
177 +conj(X1s[1])*(u11s*X2su[0]+u12s*X2su[1])
178 +conj(X1su[0])*(u12s*X2s[0]-conj(u11s)*X2s[1])
179 +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1])).imag();
180
181 dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
182 dSdpis[1]+=akappa*(
183 (conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
184 +conj(X1s[1])*(-u11s*X2su[0]-u12s*X2su[1])
185 +conj(X1su[0])*(-u12s*X2s[0]-conj(u11s)*X2s[1])
186 +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1]))).real();
187
188 dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
189 dSdpis[2]+=-akappa*(
190 conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
191 +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
192 +conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
193 +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1])).imag();
194
195 const unsigned short gindex=mu*ndirac+(idirac>>1);
196 const Complex_f gamval_c=gamval[gindex];
197 //Rescaling gind by nc
198 const unsigned short gind = gamin[gindex]<<1;
199 X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
200 X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
201
202 //If you are asked to rederive the force from Montvay and Munster you'll notice that it should be kappa*gamma
203 //but below is only gamma. We rescaled gamma by kappa already when we defined it so that's where it has gone
204 dSdpis[0]+=-(gamval_c*
205 (conj(X1s[0])* (-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
206 +conj(X1s[1])* (u11s *X2su[0]+u12s *X2su[1])
207 +conj(X1su[0])* (-u12s *X2s[0] +conj(u11s)*X2s[1])
208 +conj(X1su[1])*(u11s *X2s[0] +conj(u12s)*X2s[1]))).imag();
209 dSdpi[i+kvol*mu]=dSdpis[0];
210
211 dSdpis[1]+=(gamval_c*
212 (conj(X1s[0])* (-conj(u12s)*X2su[0] +conj(u11s)*X2su[1])
213 +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1])
214 +conj(X1su[0])* (u12s *X2s[0]+conj(u11s)*X2s[1])
215 +conj(X1su[1])* (-u11s *X2s[0]+conj(u12s)*X2s[1]))).real();
216 dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
217
218 dSdpis[2]+=-(gamval_c*
219 (conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
220 +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
221 +conj(X1su[0])*(conj(u11s)*X2s[0]+u12s *X2s[1])
222 +conj(X1su[1])*(conj(u12s)*X2s[0]-u11s *X2s[1]))).imag();
223 dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
224 }
225 }
226 }
227
243 __global__ void cuForce_t(double *dSdpi, Complex_f *u11t, Complex_f *u12t,Complex_f *X1, Complex_f *X2, Complex_f gamval[20],\
244 float *dk4m, float *dk4p, unsigned int *iu, const unsigned short gamin[16],float akappa){
245 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
246 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
247 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
248 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
249 const unsigned int gthreadId= blockId * bsize+bthreadId;
250
251 const unsigned short mu=3;
252 for(unsigned int i=gthreadId;i<kvol;i+=gsize*bsize){
253 const unsigned int ind=i+kvolHalo*mu;
254 const Complex_f u11s=u11t[ind]; const Complex_f u12s=u12t[ind];
255 const float dk4ms=dk4m[i]; const float dk4ps=dk4p[i];
256 //Up indices
257 const unsigned int uid = iu[i+kvol*mu];
258 //Similarly to Hdslash we always see idirac*nc so we do that here too.
259 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
260 Complex_f X1s[nc]; Complex_f X1su[nc];
261 Complex_f X2s[nc]; Complex_f X2su[nc];
262
263 X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
264 X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
265 X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
266 X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
267
268 // Need to be double to avoid accumulation errors
269 double dSdpis[3];
270 dSdpis[0]=dSdpi[i+kvol*mu];
271 //Multiplying by i and taking the real component is the same as taking the negative imaginary component
272 //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
273 //generators by hand
274 dSdpis[0]+=-(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
275 +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
276 +dk4ps*(conj(X1su[0])*(+u12s*X2s[0]-conj(u11s)*X2s[1])
277 +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1]))).imag();
278
279 dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
280 dSdpis[1]+=(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
281 +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1]))
282 +dk4ps*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
283 +conj(X1su[1])*( u11s *X2s[0]-conj(u12s)*X2s[1]))).real();
284
285 dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
286 dSdpis[2]+=-(dk4ms* (conj(X1s[0])* (u11s *X2su[0]+u12s *X2su[1])
287 +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
288 +dk4ps*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
289 +conj(X1su[1])* (-conj(u12s)*X2s[0]+u11s *X2s[1]))).imag();
290
291 const unsigned short gindex=mu*ndirac+(idirac>>1);
292 //Rescaling gind by nc
293 const unsigned short gind = gamin[gindex]<<1;
294 X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
295 X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
296
297 dSdpis[0]+=-(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
298 +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
299 -dk4ps*(conj(X1su[0])* (u12s *X2s[0]-conj(u11s)*X2s[1])
300 +conj(X1su[1])*(-u11s *X2s[0]-conj(u12s)*X2s[1]))).imag();
301 dSdpi[i+kvol*mu]=dSdpis[0];
302
303 dSdpis[1]+=(dk4ms*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
304 +conj(X1s[1])*(-u11s*X2su[0]-u12s *X2su[1]))
305 -dk4ps*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
306 +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1]))).real();
307 dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
308
309 dSdpis[2]+=-(dk4ms*(conj(X1s[0])*(u11s*X2su[0] +u12s *X2su[1])
310 +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
311 -dk4ps*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
312 +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1]))).imag();
313 dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
314 }
315 }
316 }
317}
318//Calling functions
319void cuGauge_force(Complex_f *ut[2],double *dSdpi,float beta,unsigned int *iu,unsigned int *id,dim3 dimGrid, dim3 dimBlock){
320 const char funcname[] = "Gauge_force";
321 int device=-1;
322 cudaGetDevice(&device);
323 Complex_f *Sigma[ndim][2], *ush[ndim][2];
324 for(unsigned short i=0;i<ndim;i++){
325#ifdef _DEBUG
326 cudaMallocManaged((void **)&Sigma[i][0],kvol*sizeof(Complex_f),cudaMemAttachGlobal);
327 cudaMallocManaged((void **)&Sigma[i][1],kvol*sizeof(Complex_f),cudaMemAttachGlobal);
328 cudaMallocManaged((void **)&ush[i][0],kvolHalo*sizeof(Complex_f),cudaMemAttachGlobal);
329 cudaMallocManaged((void **)&ush[i][1],kvolHalo*sizeof(Complex_f),cudaMemAttachGlobal);
330#else
331 cudaMallocAsync((void **)&Sigma[i][0],kvol*sizeof(Complex_f),streams[i]);
332 cudaMallocAsync((void **)&Sigma[i][1],kvol*sizeof(Complex_f),streams[i]);
333 cudaMallocAsync((void **)&ush[i][0],kvolHalo*sizeof(Complex_f),streams[i]);
334 cudaMallocAsync((void **)&ush[i][1],kvolHalo*sizeof(Complex_f),streams[i]);
335#endif
336 }
337 for(unsigned short mu=0; mu<ndim; mu++){
338 cudaMemsetAsync(Sigma[mu][0],0, kvol*sizeof(Complex_f),streams[mu]);
339 cudaMemsetAsync(Sigma[mu][1],0, kvol*sizeof(Complex_f),streams[mu]);
340 for(unsigned short nu=0; nu<ndim; nu++)
341 if(nu!=mu){
342 //The @f$-\nu@f$ Staple
343 Kernels::Plus_staple<<<dimGrid,dimBlock,0,streams[mu]>>>(mu, nu, iu, Sigma[mu][0], Sigma[mu][1],ut[0],ut[1]);
344 Kernels::Gather<<<dimGrid,dimBlock,0,streams[mu]>>>(ush[mu][0], ut[0], kvol, id, nu);
345 Kernels::Gather<<<dimGrid,dimBlock,0,streams[mu]>>>(ush[mu][1], ut[1], kvol, id, nu);
346
347#if(nproc>1)
348 //Prefetch to the CPU for until we get NCCL working
349 //cudaMemPrefetchAsync(ush[0], kvolHalo*sizeof(Complex_f),cudaCpuDeviceId,streams[0]);
350 //cudaMemPrefetchAsync(ush[1], kvolHalo*sizeof(Complex_f),cudaCpuDeviceId,streams[1]);
351 CHalo_swap_dir(ush[mu][0], 1, mu, DOWN); CHalo_swap_dir(ush[mu][1], 1, mu, DOWN);
352 //cudaMemPrefetchAsync(ush[0]+kvol, halo*sizeof(Complex_f),device,streams[0]);
353 //cudaMemPrefetchAsync(ush[1]+kvol, halo*sizeof(Complex_f),device,streams[1]);
354#endif
355 //Next up, the @f$-\nu@f$ staple
356 Kernels::Minus_staple<<<dimGrid,dimBlock,0,streams[mu]>>>(mu, nu, iu, id,Sigma[mu][0],Sigma[mu][1],\
357 ush[mu][0],ush[mu][1],ut[0],ut[1]);
358 }
359 //Now get the gauge force acting in the @f$\mu@f$ direction
360 Kernels::cuGaugeForce<<<dimGrid,dimBlock,0,streams[mu]>>>(mu,Sigma[mu][0],Sigma[mu][1],dSdpi,ut[0],ut[1],beta);
361 }
362 for(unsigned short i=0;i<ndim;i++){
363#ifdef _DEBUG
364 cudaFree(Sigma[i][0]); cudaFree(Sigma[i][1]);
365 cudaFree(ush[i][0]); cudaFree(ush[i][1]);
366#else
367 cudaFreeAsync(Sigma[i][0],streams[i]); cudaFreeAsync(Sigma[i][1],streams[i]);
368 cudaFreeAsync(ush[i][0],streams[i]); cudaFreeAsync(ush[i][1],streams[i]);
369#endif
370 }
371 cudaDeviceSynchronise();
372}
373void cuForce(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, \
374 Complex_f gamval[20],float *dk[2],unsigned int *iu,const unsigned short gamin[16],\
375 float akappa, dim3 dimGrid, dim3 dimBlock){
376 const char *funcname = "Force";
377 //X1=(M†M)^{1} Phi
378 // Transpose_z(X1,ndirac*nc,kvol); Transpose_z(X2,ndirac*nc,kvol);
379 cudaDeviceSynchronise();
380#pragma unroll
381 for(unsigned short mu=0;mu<3;mu++){
382 Kernels::cuForce_s<<<dimGrid,dimBlock,0,streams[mu]>>>(dSdpi,ut[0],ut[1],X1,X2,gamval,iu,gamin,akappa,mu);
383 }
384 //Set stream for time direction
385 unsigned short mu=3;
386 Kernels::cuForce_t<<<dimGrid,dimBlock,0,streams[mu]>>>(dSdpi,ut[0],ut[1],X1,X2,gamval,dk[0],dk[1],iu,gamin,akappa);
387 cudaDeviceSynchronise();
388}
conj
__device__ __forceinline__ T conj(const T &z)
Complex Conjugation.
Definition cusu2hmc.cu:33
Kernels::Gather
__global__ void Gather(T *x, T *y, const unsigned int n, unsigned int *table, const unsigned short mu)
Extracts all the single precision gauge links in the direction only.
Definition cuforce.cu:116
cuForce
void cuForce(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk[2], unsigned int *iu, const unsigned short gamin[16], float akappa, dim3 dimGrid, dim3 dimBlock)
Calculates the force at each intermediate time.
Definition cuforce.cu:373
cuGauge_force
void cuGauge_force(Complex_f *ut[2], double *dSdpi, float beta, unsigned int *iu, unsigned int *id, dim3 dimGrid, dim3 dimBlock)
Calculate the gauge contribution to the force.
Definition cuforce.cu:319
Kernels::Plus_staple
__global__ void Plus_staple(const int mu, const int nu, unsigned int *iu, Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11t, Complex_f *u12t)
Calculates the staple in the positive direction.
Definition cuforce.cu:23
Kernels::cuForce_s
__global__ void cuForce_s(double *dSdpi, Complex_f *u11t, Complex_f *u12t, Complex_f *X1, Complex_f *X2, Complex_f gamval[20], unsigned int *iu, const unsigned short gamin[16], float akappa, const unsigned short mu)
Calculates the force at each intermediate time.
Definition cuforce.cu:146
Kernels::cuGaugeForce
__global__ void cuGaugeForce(int mu, Complex_f *Sigma11, Complex_f *Sigma12, double *dSdpi, Complex_f *u11t, Complex_f *u12t, float beta)
Calculates the gauge force due to the Wilson Action at each intermediate time.
Definition cuforce.cu:89
Kernels::Minus_staple
__global__ void Minus_staple(const int mu, const int nu, unsigned int *iu, unsigned int *id, Complex_f *Sigma11, Complex_f *Sigma12, Complex_f *u11sh, Complex_f *u12sh, Complex_f *u11t, Complex_f *u12t)
Calculates the staple in the negative direction.
Definition cuforce.cu:55
Kernels::cuForce_t
__global__ void cuForce_t(double *dSdpi, Complex_f *u11t, Complex_f *u12t, Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk4m, float *dk4p, unsigned int *iu, const unsigned short gamin[16], float akappa)
Calculates the force at each intermediate time.
Definition cuforce.cu:243
CHalo_swap_dir
int CHalo_swap_dir(Complex_f *c, int ncpt, int idir, int layer)
Swaps the halos along the axis given by idir in the direction given by layer.
matrices.h
Matrix multiplication and related declarations.
Kernels
CUDA Kernels.
Definition cubosonic.cu:47
DOWN
#define DOWN
Flag for send down.
Definition par_mpi.h:37
nc
#define nc
Colours.
Definition sizes.h:182
kvol
#define kvol
Sublattice volume.
Definition sizes.h:163
ndirac
#define ndirac
Dirac indices.
Definition sizes.h:186
cudaDeviceSynchronise
#define cudaDeviceSynchronise()
Get rid of that bastardised yankee English.
Definition sizes.h:53
Complex_f
#define Complex_f
Single precision complex number.
Definition sizes.h:62
dimGrid
dim3 dimGrid
Default grid size. First component is normally nt. Second and third depend whatever is needed to get ...
Definition cusu2hmc.cu:27
ndim
#define ndim
Dimensions.
Definition sizes.h:188
kvolHalo
#define kvolHalo
Subvolume + halo size.
Definition sizes.h:234
dimBlock
dim3 dimBlock
Default block size. Usually 128.
Definition cusu2hmc.cu:25
su2hmc.h
Function declarations for most of the routines.
streams
cudaStream_t streams[ndirac *ndim *nadj]
An array of concurrent GPU streams to keep it busy.
Definition cusu2hmc.cu:29