su2hmc
Loading...
Searching...
No Matches
cumatrices.cu
Go to the documentation of this file.
1
6#include <assert.h>
7#include <su2hmc.h>
8#include <matrices.h>
9#include <thrust_complex.h>
10namespace Device{
20 template <typename T,unsigned int bsize>
21 __device__ void warpReduce_sum(volatile T* sdata, const unsigned int tid){
22 if(bsize >= 64) sdata[tid] += sdata[tid + 32];
23 if(bsize >= 32) sdata[tid] += sdata[tid + 16];
24 if(bsize >= 16) sdata[tid] += sdata[tid + 8];
25 if(bsize >= 8) sdata[tid] += sdata[tid + 4];
26 if(bsize >= 4) sdata[tid] += sdata[tid + 2];
27 if(bsize >= 2) sdata[tid] += sdata[tid + 1];
28 }
29}
30namespace Kernels{
47 template <typename T>
48 __global__ void cuDslash(complex<T> *phi, complex<T> *r, complex<T> *u11t, complex<T> *u12t,const unsigned int *iu, const unsigned int *id,\
49 complex<T> gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa){
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 }
158
174 template <typename T>
175 __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,\
176 complex<T> gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const Complex_f jqq, const float akappa){
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 }
285
301 template <typename T>
302 __global__ void cuHdslash(complex<T> *phi, const complex<T> *r, const complex<T> *u11t, const complex<T> *u12t,unsigned int *iu, unsigned int *id,\
303 __constant__ complex<T> gamval[20], const unsigned short gamin[16], const T *dk4m, const T *dk4p, const __grid_constant__ float akappa){
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 }
382
397 template <typename T>
398 __global__ void cuHdslashd(complex<T> *phi, const complex<T>* r, const complex<T>* u11t, const complex<T>* u12t,unsigned int* iu, unsigned int* id,\
399 __constant__ complex<T> gamval[20], const unsigned short gamin[16], const T* dk4m, const T* dk4p, const __grid_constant__ float akappa){
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 }
478
489 template <typename T>
490 __global__ void Transpose(T *out, const T *in, const int fast_in, const int fast_out){
491 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
492 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
493 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
494 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
495 const unsigned int gthreadId= blockId * bsize+bthreadId;
496
497 //The if/else here is only to ensure we maximise GPU bandwidth
498 //Typically this is used to write back to the AoS/Coalseced format
499 if(fast_out>fast_in){
500 for(unsigned int x=gthreadId;x<fast_out;x+=gsize*bsize)
501 for(unsigned int y=0; y<fast_in;y++)
502 out[y*fast_out+x]=in[x*fast_in+y];
503 }
504 //Typically this is used to write back to the SoA/saved config format
505 else{
506 for(unsigned int x=0; x<fast_out;x++)
507 for(unsigned int y=gthreadId;y<fast_in;y+=gsize*bsize)
508 out[y*fast_out+x]=in[x*fast_in+y];
509 }
510 }
511
520 __global__ void Mixed_Sumto(double *d, float *f, const unsigned int n){
521 const unsigned int gsize = gridDim.x*gridDim.y*gridDim.z;
522 const unsigned int bsize = blockDim.x*blockDim.y*blockDim.z;
523 const unsigned int blockId = blockIdx.x+ blockIdx.y * gridDim.x+ gridDim.x * gridDim.y * blockIdx.z;
524 const unsigned int bthreadId= (threadIdx.z * blockDim.y+ threadIdx.y)* blockDim.x+ threadIdx.x;
525 const unsigned int gthreadId= blockId * bsize+bthreadId;
526
527 for(unsigned int i=gthreadId; i<n;i+=bsize*gsize)
528 d[i]+=(double)f[i];
529 return;
530 }
531
542 template <typename T,unsigned int bsize>
543 __global__ void reduce_sum(T *g_in_data, T *g_out_data, const unsigned int n){
544 extern __shared__ T sdata[]; // stored in the shared memory
545
546 // Each thread loading one element from global onto shared memory
547 const unsigned short tid = threadIdx.x;
548 unsigned int i = blockIdx.x*(bsize*2) + tid;
549 const unsigned int gridSize = blockDim.x * 2 * gridDim.x;
550 sdata[tid] = 0;
551 while (i < n) {
552 sdata[tid] += g_in_data[i];
553 if (i + bsize < n) {
554 sdata[tid] += g_in_data[i + bsize];
555 }
556 i += gridSize;
557 }
558 __syncthreads();
559
560 // Perform reductions in steps, reducing thread synchronization
561 // CUDA warps aren't guaranteed to be 32, and AMD vary. warpSize from the compiler fixes that.
562 for(unsigned int s=bsize/2;s>=warpSize;s>>=1){
563 if (tid < s)
564 sdata[tid] += sdata[tid + s]; __syncthreads();
565 }
566#if defined(__HIP_PLATFORM_AMD__) || defined(__HIP_PLATFORM_HCC__)
567#define SU2_WARP_MASK 0xffffffffffffffffULL // 64-bit lane mask (HIP/AMD)
568#else
569#define SU2_WARP_MASK 0xffffffffU // 32-bit lane mask (CUDA)
570#endif
571 if (tid < warpSize)
572 //Device::warpReduce_sum<T,bsize>(sdata, tid);
573 {
574 T val = sdata[tid];
575 for (int offset = warpSize/2; offset > 0; offset >>= 1)
576 val += __shfl_down_sync(SU2_WARP_MASK, val, offset);
577 if (tid == 0) sdata[0] = val;
578 }
579
580 if (tid == 0){
581 g_out_data[blockIdx.x] = sdata[0];
582 }
583 }
584}
585//Calling Functions
586//================
587double cureduce_sum_d(double *input, const unsigned int n,const unsigned short stream){
588 const unsigned int bsize=256;
589 unsigned int gsize=(n + (2 * bsize) - 1) / (2 * bsize);
590 double *cachein, *cacheout;
591 cudaMallocAsync(&cacheout,gsize*sizeof(double),streams[stream]);
592 Kernels::reduce_sum<double,bsize><<<gsize,bsize,bsize*sizeof(double),streams[stream]>>>(input,cacheout,n);
593 while(gsize>1){
594 cudaMallocAsync(&cachein,gsize*sizeof(double),streams[stream]);
595 cudaMemcpyAsync(cachein,cacheout,gsize*sizeof(double),cudaMemcpyDefault,streams[stream]);
596 cudaFreeAsync(cacheout,streams[stream]);
597 gsize>>=1;
598 cudaMallocAsync(&cacheout,gsize*sizeof(double),streams[stream]);
599 Kernels::reduce_sum<double,bsize><<<gsize,bsize,bsize*sizeof(double),streams[stream]>>>(cachein,cacheout,gsize);
600 cudaFreeAsync(cachein,streams[stream]);
601 }
602 double output=0;
603 cudaStreamSynchronize(streams[stream]);
604 cudaMemcpyAsync(&output,cacheout,sizeof(double),cudaMemcpyDefault,streams[stream]);
605 cudaStreamSynchronize(streams[stream]);
606 cudaFreeAsync(cacheout,streams[stream]);
607 return output;
608}
609void cuDslash(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,
610 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa,
611 dim3 dimGrid, dim3 dimBlock){
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}
623void cuDslashd(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,
624 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa,
625 dim3 dimGrid, dim3 dimBlock){
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}
637void cuHdslash(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,
638 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa,
639 dim3 dimGrid, dim3 dimBlock){
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}
651void cuHdslashd(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,
652 Complex gamval[20], const unsigned short gamin[16],double *dk[nc], float akappa,
653 dim3 dimGrid, dim3 dimBlock){
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}
666
667//Float editions
668void cuDslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,
669 Complex_f gamval[20],const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa,
670 dim3 dimGrid, dim3 dimBlock){
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}
682void cuDslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,
683 Complex_f gamval[20],const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa,
684 dim3 dimGrid, dim3 dimBlock){
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}
696void cuHdslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id, Complex_f gamval[20],
697 const unsigned short gamin[16], float *dk[nc], float akappa, dim3 dimGrid, dim3 dimBlock){
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}
711void cuHdslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,
712 Complex_f gamval[20],const unsigned short gamin[16],float *dk[nc], float akappa,dim3 dimGrid, dim3 dimBlock){
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}
724
725void cuTranspose_z(Complex *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
726 Complex *holder;
727 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(Complex));
728 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(Complex),cudaMemcpyDefault);
729 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
730 cudaFree(holder);
731}
732void cuTranspose_c(Complex_f *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
733 Complex_f *holder;
734 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(Complex_f));
735 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(Complex_f),cudaMemcpyDefault);
736 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
737 cudaDeviceSynchronise();
738 cudaFree(holder);
739}
740void cuTranspose_d(double *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
741 double *holder;
742 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(double));
743 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(double),cudaMemcpyDefault);
744 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
745 cudaFree(holder);
746}
747void cuTranspose_f(float *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
748 float *holder;
749 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(float));
750 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(float),cudaMemcpyDefault);
751 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
752 cudaFree(holder);
753}
754void cuTranspose_I(int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
755 int *holder;
756 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(int));
757 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(int),cudaMemcpyDefault);
758 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
759 cudaFree(holder);
760}
761void cuTranspose_U(unsigned int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock){
762 unsigned int *holder;
763 cudaMalloc((void **)&holder,fast_in*fast_out*sizeof(unsigned int));
764 cudaMemcpy(holder,out,fast_in*fast_out*sizeof(unsigned int),cudaMemcpyDefault);
765 Kernels::Transpose<<<dimGrid,dimBlock>>>(out,holder,fast_in,fast_out);
766 cudaFree(holder);
767}
768
769void cuMixed_Sumto(double *d, float *f,const unsigned int n,const dim3 dimGrid,const dim3 dimBlock){
770 Kernels::Mixed_Sumto<<<dimGrid,dimBlock>>>(d,f,n);
771 return;
772}
CPYERROR
#define CPYERROR
Copy failed.
Definition errorcodes.h:66
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
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
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
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
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
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
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
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
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
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
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
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
cureduce_sum_d
double cureduce_sum_d(double *input, const unsigned int n, const unsigned short stream)
Sum all terms in an array of doubles.
Definition cumatrices.cu:587
cuTranspose_f
void cuTranspose_f(float *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:747
cuTranspose_d
void cuTranspose_d(double *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:740
cuTranspose_U
void cuTranspose_U(unsigned int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:761
Device::warpReduce_sum
__device__ void warpReduce_sum(volatile T *sdata, const unsigned int tid)
Performs a warp reduction for sum.
Definition cumatrices.cu:21
Kernels::reduce_sum
__global__ void reduce_sum(T *g_in_data, T *g_out_data, const unsigned int n)
Performs a block reduction for sum.
Definition cumatrices.cu:543
conj
__device__ __forceinline__ T conj(const T &z)
Complex Conjugation.
Definition cusu2hmc.cu:33
cuTranspose_z
void cuTranspose_z(Complex *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:725
Kernels::Mixed_Sumto
__global__ void Mixed_Sumto(double *d, float *f, const unsigned int n)
Sums a float array into a double array.
Definition cumatrices.cu:520
cuTranspose_c
void cuTranspose_c(Complex_f *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:732
cuTranspose_I
void cuTranspose_I(int *out, const int fast_in, const int fast_out, const dim3 dimGrid, const dim3 dimBlock)
In place transpose used to convert from AoS to SoA memory layout.
Definition cumatrices.cu:754
Kernels::Transpose
__global__ void Transpose(T *out, const T *in, const int fast_in, const int fast_out)
Swaps the order of the gauge field so that it is now SoA instead of AoS and it is nice and coalesced ...
Definition cumatrices.cu:490
cuMixed_Sumto
void cuMixed_Sumto(double *d, float *f, const unsigned int n, const dim3 dimGrid, const dim3 dimBlock)
Add a single to a double value, and save the output in the double array For complex valued arrays,...
Definition cumatrices.cu:769
matrices.h
Matrix multiplication and related declarations.
Device
CUDA Device code.
Definition cubosonic.cu:13
Kernels
CUDA Kernels.
Definition cubosonic.cu:47
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
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
thrust_complex.h
Complex Header for CUDA. Sets macros for C compatability.