su2hmc
Loading...
Searching...
No Matches
matrices.c
Go to the documentation of this file.
1
12#include <assert.h>
13#include <matrices.h>
14//TODO: Check and see are there any terms we are evaluating twice in the same loop
15//and use a variable to hold them instead to reduce the number of evaluations.
16int Dslash(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu,unsigned int *id,\
17 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa){
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}
135int Dslashd(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,\
136 Complex gamval[20], const unsigned short gamin[16], double *dk[nc],Complex_f jqq, float akappa){
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}
254int Hdslash(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,\
255 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa){
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}
338int Hdslashd(Complex *phi, Complex *r, Complex *ut[nc],unsigned int *iu,unsigned int *id,\
339 Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa){
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}
423//Float Versions
424//int Dslash_f(Complex_f *phi, Complex_f *r){
425int Dslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu, unsigned int *id,\
426 Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa){
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}
544int Dslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,\
545 Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa){
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}
663int Hdslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,\
664 Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa){
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}
746int Hdslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc],unsigned int *iu,unsigned int *id,\
747 Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa){
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}
834
835
836inline void Transpose_c(Complex_f *out, const int fast_in, const int fast_out){
837 const volatile char funcname[]="Transpose_c";
838
839#ifdef USE_GPU
840 cuTranspose_c(out,fast_in,fast_out,dimGrid,dimBlock);
841#else
842 Complex_f *in = (Complex_f *)aligned_alloc(AVX,fast_in*fast_out*sizeof(Complex_f));
843 memcpy(in,out,fast_in*fast_out*sizeof(Complex_f));
844 //Typically this is used to write back to the AoS/Coalseced format
845 if(fast_out>fast_in){
846 for(int x=0;x<fast_out;x++)
847 for(int y=0; y<fast_in;y++)
848 out[y*fast_out+x]=in[x*fast_in+y];
849 }
850 //Typically this is used to write back to the SoA/saved config format
851 else{
852 for(int x=0; x<fast_out;x++)
853 for(int y=0;y<fast_in;y++)
854 out[y*fast_out+x]=in[x*fast_in+y];
855 }
856 free(in);
857#endif
858}
859inline void Transpose_z(Complex *out, const int fast_in, const int fast_out){
860 const volatile char funcname[]="Transpose_c";
861
862#ifdef USE_GPU
863 cuTranspose_z(out,fast_in,fast_out,dimGrid,dimBlock);
864#else
865 Complex *in = (Complex *)aligned_alloc(AVX,fast_in*fast_out*sizeof(Complex));
866 memcpy(in,out,fast_in*fast_out*sizeof(Complex));
867 //Typically this is used to write back to the AoS/Coalseced format
868 if(fast_out>fast_in){
869 for(int x=0;x<fast_out;x++)
870 for(int y=0; y<fast_in;y++)
871 out[y*fast_out+x]=in[x*fast_in+y];
872 }
873 //Typically this is used to write back to the SoA/saved config format
874 else{
875 for(int x=0; x<fast_out;x++)
876 for(int y=0;y<fast_in;y++)
877 out[y*fast_out+x]=in[x*fast_in+y];
878 }
879 free(in);
880#endif
881}
882inline void Transpose_f(float *out, const int fast_in, const int fast_out){
883 const char funcname[]="Transpose_f";
884
885#ifdef USE_GPU
886 cuTranspose_f(out,fast_in,fast_out,dimGrid,dimBlock);
887#else
888 float *in = (float *)aligned_alloc(AVX,fast_in*fast_out*sizeof(float));
889 memcpy(in,out,fast_in*fast_out*sizeof(float));
890 //Typically this is used to write back to the AoS/Coalseced format
891 if(fast_out>fast_in){
892 for(int x=0;x<fast_out;x++)
893 for(int y=0; y<fast_in;y++)
894 out[y*fast_out+x]=in[x*fast_in+y];
895 }
896 //Typically this is used to write back to the SoA/saved config format
897 else{
898 for(int x=0; x<fast_out;x++)
899 for(int y=0;y<fast_in;y++)
900 out[y*fast_out+x]=in[x*fast_in+y];
901 }
902 free(in);
903#endif
904}
905inline void Transpose_d(double *out, const int fast_in, const int fast_out){
906 const char funcname[]="Transpose_f";
907
908#ifdef USE_GPU
909 cuTranspose_d(out,fast_in,fast_out,dimGrid,dimBlock);
910#else
911 double *in = (double *)aligned_alloc(AVX,fast_in*fast_out*sizeof(double));
912 memcpy(in,out,fast_in*fast_out*sizeof(double));
913 //Typically this is used to write back to the AoS/Coalseced format
914 if(fast_out>fast_in){
915 for(int x=0;x<fast_out;x++)
916 for(int y=0; y<fast_in;y++)
917 out[y*fast_out+x]=in[x*fast_in+y];
918 }
919 //Typically this is used to write back to the SoA/saved config format
920 else{
921 for(int x=0; x<fast_out;x++)
922 for(int y=0;y<fast_in;y++)
923 out[y*fast_out+x]=in[x*fast_in+y];
924 }
925 free(in);
926#endif
927}
928inline void Transpose_I(int *out, const int fast_in, const int fast_out){
929 const char funcname[]="Transpose_I";
930
931#ifdef USE_GPU
932 cuTranspose_I(out,fast_in,fast_out,dimGrid,dimBlock);
933#else
934 int *in = (int *)aligned_alloc(AVX,fast_in*fast_out*sizeof(int));
935 memcpy(in,out,fast_in*fast_out*sizeof(int));
936 //Typically this is used to write back to the AoS/Coalseced format
937 if(fast_out>fast_in){
938 for(int x=0;x<fast_out;x++)
939 for(int y=0; y<fast_in;y++)
940 out[y*fast_out+x]=in[x*fast_in+y];
941 }
942 //Typically this is used to write back to the SoA/saved config format
943 else{
944 for(int x=0; x<fast_out;x++)
945 for(int y=0;y<fast_in;y++)
946 out[y*fast_out+x]=in[x*fast_in+y];
947 }
948 free(in);
949#endif
950}
951inline void Transpose_U(unsigned int *out, const int fast_in, const int fast_out){
952 const char funcname[]="Transpose_I";
953
954#ifdef USE_GPU
955 cuTranspose_U(out,fast_in,fast_out,dimGrid,dimBlock);
956#else
957 unsigned int *in = (unsigned int *)aligned_alloc(AVX,fast_in*fast_out*sizeof(unsigned int));
958 memcpy(in,out,fast_in*fast_out*sizeof(unsigned int));
959 //Typically this is used to write back to the AoS/Coalseced format
960 if(fast_out>fast_in){
961 for(unsigned int x=0;x<fast_out;x++)
962 for(unsigned int y=0; y<fast_in;y++)
963 out[y*fast_out+x]=in[x*fast_in+y];
964 }
965 //Typically this is used to write back to the SoA/saved config format
966 else{
967 for(unsigned int x=0; x<fast_out;x++)
968 for(unsigned int y=0;y<fast_in;y++)
969 out[y*fast_out+x]=in[x*fast_in+y];
970 }
971 free(in);
972#endif
973}
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
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 in single precision.
Definition matrices.c:425
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
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 in single precision.
Definition matrices.c:544
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
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 in single precision.
Definition matrices.c:746
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 in double precision.
Definition matrices.c:338
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
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
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 in double precision.
Definition matrices.c:254
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 in double precision.
Definition matrices.c:135
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
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 in single precision.
Definition matrices.c:663
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 in double precision.
Definition matrices.c:16
Transpose_c
void Transpose_c(Complex_f *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:836
Transpose_U
void Transpose_U(unsigned int *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:951
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
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
Transpose_f
void Transpose_f(float *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:882
Transpose_z
void Transpose_z(Complex *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:859
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
Transpose_d
void Transpose_d(double *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:905
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
Transpose_I
void Transpose_I(int *out, const int fast_in, const int fast_out)
In place transpose used to convert from AoS to SoA memory layout.
Definition matrices.c:928
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.
ZHalo_swap_all
int ZHalo_swap_all(Complex *z, int ncpt)
Calls the functions to send data to both the up and down halos.
matrices.h
Matrix multiplication and related declarations.
AVX
#define AVX
Alignment of arrays. 64 for AVX-512, 32 for AVX/AVX2. 16 for SSE. Since AVX is standard on modern x86...
Definition sizes.h:279
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
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