Bridge++  Ver.2.1.3
mult_Domainwall_5din_LUinv_openacc-inc.h
Go to the documentation of this file.
1 
10 #ifndef MULT_DOMAINWALL_5DIN_LUINV_ACC_INCLUDED
11 #define MULT_DOMAINWALL_5DIN_LUINV_ACC_INCLUDED
12 
13 //====================================================================
14 void mult_domainwall_5din_LUinv_dirac(
15  real_t *RESTRICT vp, real_t *RESTRICT wp,
16  int Ns, int *Nsize,
17  real_t *e, real_t *f,
18  real_t *dpinv, real_t *dm, real_t alpha)
19 {
20  int Nx = Nsize[0];
21  int Ny = Nsize[1];
22  int Nz = Nsize[2];
23  int Nt = Nsize[3];
24  int Nst = Nx * Ny * Nz * Nt;
25  int Nst_pad = CEIL_NWP(Nst);
26 
27  int Nin5 = NVCD * Ns;
28  int size = Nin5 * Nst_pad;
29 
30 #pragma acc data present(vp[0:size], wp[0:size]) \
31  copyin(Nst, Nst_pad, Nin5, Ns, e[0:Ns-1], f[0:Ns-1], \
32  dpinv[0:Ns], dm[0:Ns], alpha)
33  {
34 
35 #pragma acc parallel num_workers(NUM_WORKERS) vector_length(VECTOR_LENGTH)
36  {
37 
38 #pragma acc loop gang worker vector
39  for (int idx = 0; idx < NVC * Nst_pad; ++idx) {
40  int idx2_wp = idx / NWP;
41  int idx_in = idx % NWP;
42  int ivc = idx2_wp % NVC;
43  int idx_out = idx2_wp / NVC;
44  int site = idx_in + NWP * idx_out;
45  if(site < Nst){
46 
47  real_t vt[ND], yt[ND], xt[ND];
48 
49  int is = 0;
50  for(int id = 0; id < ND; ++id){
51  int ivcd = ivc + NVC * id;
52  vt[id] = wp[IDX2(Nin5, (ivcd + NVCD * is), site)];
53  }
54  for(int id = 0; id < ND; ++id){
55  int ivcd = ivc + NVC * id;
56  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
57  }
58 
59  for(int id = 0; id < ND; ++id){
60  yt[id] = e[0] * vt[id];
61  }
62 
63  for (int is = 1; is < Ns-1; ++is) {
64 
65  for(int id = 0; id < ND; ++id){
66  xt[id] = vt[id];
67  }
68 
69  for(int id = 0; id < ND; ++id){
70  int ivcd = ivc + NVC * id;
71  vt[id] = wp[IDX2(Nin5, (ivcd + NVCD * is), site)];
72  }
73 
74  real_t a = real_t(0.5) * dm[is] * dpinv[is-1];
75 
76  vt[0] += a * (xt[0] + xt[2]);
77  vt[1] += a * (xt[1] + xt[3]);
78  vt[2] += a * (xt[2] + xt[0]);
79  vt[3] += a * (xt[3] + xt[1]);
80 
81  for(int id = 0; id < ND; ++id){
82  int ivcd = ivc + NVC * id;
83  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
84  }
85 
86  for(int id = 0; id < ND; ++id){
87  yt[id] += e[is] * vt[id];
88  }
89 
90  }
91 
92  is = Ns-1;
93 
94  for(int id = 0; id < ND; ++id){
95  xt[id] = vt[id];
96  }
97 
98  for(int id = 0; id < ND; ++id){
99  int ivcd = ivc + NVC * id;
100  vt[id] = wp[IDX2(Nin5, (ivcd + NVCD * is), site)];
101  }
102 
103  real_t a = real_t(0.5) * dm[is] * dpinv[is-1];
104 
105  vt[0] += a * (xt[0] + xt[2]);
106  vt[1] += a * (xt[1] + xt[3]);
107  vt[2] += a * (xt[2] + xt[0]);
108  vt[3] += a * (xt[3] + xt[1]);
109 
110  vt[0] += -0.5 * (yt[0] - yt[2]);
111  vt[1] += -0.5 * (yt[1] - yt[3]);
112  vt[2] += -0.5 * (yt[2] - yt[0]);
113  vt[3] += -0.5 * (yt[3] - yt[1]);
114 
115  for(int id = 0; id < ND; ++id){
116  int ivcd = ivc + NVC * id;
117  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
118  }
119  // L_inv completed
120 
121  is = Ns-1;
122 
123  a = dpinv[Ns-1];
124  real_t f1 = 0.5 * (1.0 + alpha);
125  real_t f2 = 0.5 * (-1.0 + alpha);
126 
127  real_t vt1, vt2, vt3, vt4;
128  vt1 = vp[IDX2(Nin5, (ID1 + ivc + NVCD * is), site)];
129  vt2 = vp[IDX2(Nin5, (ID2 + ivc + NVCD * is), site)];
130  vt3 = vp[IDX2(Nin5, (ID3 + ivc + NVCD * is), site)];
131  vt4 = vp[IDX2(Nin5, (ID4 + ivc + NVCD * is), site)];
132  vt[0] = a * (f1 * vt1 + f2 * vt3);
133  vt[1] = a * (f1 * vt2 + f2 * vt4);
134  vt[2] = a * (f1 * vt3 + f2 * vt1);
135  vt[3] = a * (f1 * vt4 + f2 * vt2);
136 
137  for(int id = 0; id < ND; ++id){
138  int ivcd = ivc + NVC * id;
139  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
140  }
141 
142  yt[0] = 0.5 * (vt[0] + vt[2]);
143  yt[1] = 0.5 * (vt[1] + vt[3]);
144  yt[2] = 0.5 * (vt[2] + vt[0]);
145  yt[3] = 0.5 * (vt[3] + vt[1]);
146 
147  for (int is = Ns-2; is >= 0; --is) {
148 
149  for(int id = 0; id < ND; ++id){
150  xt[id] = vt[id];
151  }
152 
153  for(int id = 0; id < ND; ++id){
154  int ivcd = ivc + NVC * id;
155  vt[id] = vp[IDX2(Nin5, (ivcd + NVCD * is), site)];
156  }
157 
158  real_t a = real_t(0.5) * dm[is];
159 
160  vt[0] += a * (xt[0] - xt[2]);
161  vt[1] += a * (xt[1] - xt[3]);
162  vt[2] += a * (xt[2] - xt[0]);
163  vt[3] += a * (xt[3] - xt[1]);
164 
165  for(int id = 0; id < ND; ++id){
166  vt[id] += - f[is] * yt[id];
167  }
168 
169  real_t aa = dpinv[is];
170 
171  for(int id = 0; id < ND; ++id){
172  vt[id] *= aa;
173  }
174 
175  if(is == 0){
176  real_t f1 = 0.5 * (1.0 + alpha);
177  real_t f2 = 0.5 * (1.0 - alpha);
178  vt1 = f1 * vt[0] + f2 * vt[2];
179  vt2 = f1 * vt[1] + f2 * vt[3];
180  vt3 = f1 * vt[2] + f2 * vt[0];
181  vt4 = f1 * vt[3] + f2 * vt[1];
182  vp[IDX2(Nin5, (ID1 + ivc + NVCD * is), site)] = vt1;
183  vp[IDX2(Nin5, (ID2 + ivc + NVCD * is), site)] = vt2;
184  vp[IDX2(Nin5, (ID3 + ivc + NVCD * is), site)] = vt3;
185  vp[IDX2(Nin5, (ID4 + ivc + NVCD * is), site)] = vt4;
186  }else{
187  for(int id = 0; id < ND; ++id){
188  int ivcd = ivc + NVC * id;
189  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
190  }
191  }
192 
193  }
194 
195  }
196  } // idx loop
197 
198  }
199  }
200 
201 }
202 
203 //====================================================================
204 void mult_domainwall_5din_LUdaginv_dirac(
205  real_t *RESTRICT vp, real_t *RESTRICT wp,
206  int Ns, int *Nsize,
207  real_t *e, real_t *f,
208  real_t *dpinv, real_t *dm, real_t alpha)
209 {
210  int Nx = Nsize[0];
211  int Ny = Nsize[1];
212  int Nz = Nsize[2];
213  int Nt = Nsize[3];
214  int Nst = Nx * Ny * Nz * Nt;
215  int Nst_pad = CEIL_NWP(Nst);
216 
217  int Nin5 = NVCD * Ns;
218  int size = Nin5 * Nst_pad;
219 
220 #pragma acc data present(vp[0:size], wp[0:size]) \
221  copyin(Nst, Nst_pad, Nin5, Ns, e[0:Ns-1], f[0:Ns-1], \
222  dpinv[0:Ns], dm[0:Ns], alpha)
223  {
224 
225 #pragma acc parallel num_workers(NUM_WORKERS) vector_length(VECTOR_LENGTH)
226  {
227 
228 #pragma acc loop gang worker vector
229  for (int idx = 0; idx < NVC * Nst_pad; ++idx) {
230  int idx2_wp = idx / NWP;
231  int idx_in = idx % NWP;
232  int ivc = idx2_wp % NVC;
233  int idx_out = idx2_wp / NVC;
234  int site = idx_in + NWP * idx_out;
235  if(site < Nst){
236 
237  real_t vt[ND], yt[ND], xt[ND];
238 
239  int is = 0;
240 
241  real_t a = dpinv[0];
242  real_t f1 = 0.5 * (1.0 + alpha);
243  real_t f2 = 0.5 * (1.0 - alpha);
244  {
245  real_t vt1, vt2, vt3, vt4;
246  vt1 = wp[IDX2(Nin5, (ID1 + ivc + NVCD * is), site)];
247  vt2 = wp[IDX2(Nin5, (ID2 + ivc + NVCD * is), site)];
248  vt3 = wp[IDX2(Nin5, (ID3 + ivc + NVCD * is), site)];
249  vt4 = wp[IDX2(Nin5, (ID4 + ivc + NVCD * is), site)];
250  vt[0] = a * (f1 * vt1 + f2 * vt3);
251  vt[1] = a * (f1 * vt2 + f2 * vt4);
252  vt[2] = a * (f1 * vt3 + f2 * vt1);
253  vt[3] = a * (f1 * vt4 + f2 * vt2);
254  }
255 
256  for(int id = 0; id < ND; ++id){
257  int ivcd = ivc + NVC * id;
258  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
259  }
260 
261  for(int id = 0; id < ND; ++id){
262  yt[id] = f[0] * vt[id];
263  }
264 
265  for (int is = 1; is < Ns-1; ++is) {
266 
267  for(int id = 0; id < ND; ++id){
268  xt[id] = vt[id];
269  }
270 
271  for(int id = 0; id < ND; ++id){
272  int ivcd = ivc + NVC * id;
273  vt[id] = wp[IDX2(Nin5, (ivcd + NVCD * is), site)];
274  }
275 
276  real_t a = real_t(0.5) * dm[is - 1];
277 
278  vt[0] += a * (xt[0] - xt[2]);
279  vt[1] += a * (xt[1] - xt[3]);
280  vt[2] += a * (xt[2] - xt[0]);
281  vt[3] += a * (xt[3] - xt[1]);
282 
283  real_t aa = dpinv[is];
284 
285  for(int id = 0; id < ND; ++id){
286  vt[id] *= aa;
287  }
288 
289  for(int id = 0; id < ND; ++id){
290  int ivcd = ivc + NVC * id;
291  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
292  }
293 
294  for(int id = 0; id < ND; ++id){
295  yt[id] += f[is] * vt[id];
296  }
297 
298  }
299 
300  is = Ns-1;
301 
302  for(int id = 0; id < ND; ++id){
303  xt[id] = vt[id];
304  }
305 
306  for(int id = 0; id < ND; ++id){
307  int ivcd = ivc + NVC * id;
308  vt[id] = wp[IDX2(Nin5, (ivcd + NVCD * is), site)];
309  }
310 
311  a = real_t(0.5) * dm[is - 1];
312 
313  vt[0] += a * (xt[0] - xt[2]);
314  vt[1] += a * (xt[1] - xt[3]);
315  vt[2] += a * (xt[2] - xt[0]);
316  vt[3] += a * (xt[3] - xt[1]);
317 
318  vt[0] += -0.5 * (yt[0] + yt[2]);
319  vt[1] += -0.5 * (yt[1] + yt[3]);
320  vt[2] += -0.5 * (yt[2] + yt[0]);
321  vt[3] += -0.5 * (yt[3] + yt[1]);
322 
323  real_t aa = dpinv[is];
324 
325  for(int id = 0; id < ND; ++id){
326  vt[id] *= aa;
327  }
328 
329  real_t ff1 = 0.5 * ( 1.0 + alpha);
330  real_t ff2 = 0.5 * (-1.0 + alpha);
331  {
332  real_t vt1, vt2, vt3, vt4;
333  vt1 = ff1 * vt[0] + ff2 * vt[2];
334  vt2 = ff1 * vt[1] + ff2 * vt[3];
335  vt3 = ff1 * vt[2] + ff2 * vt[0];
336  vt4 = ff1 * vt[3] + ff2 * vt[1];
337  vp[IDX2(Nin5, (ID1 + ivc + NVCD*is), site)] = vt1;
338  vp[IDX2(Nin5, (ID2 + ivc + NVCD*is), site)] = vt2;
339  vp[IDX2(Nin5, (ID3 + ivc + NVCD*is), site)] = vt3;
340  vp[IDX2(Nin5, (ID4 + ivc + NVCD*is), site)] = vt4;
341  }
342  // Udag_inv completed
343 
344  is = Ns-1;
345 
346  for(int id = 0; id < ND; ++id){
347  int ivcd = ivc + NVC * id;
348  vt[id] = vp[IDX2(Nin5, (ivcd + NVCD * is), site)];
349  }
350 
351  for(int id = 0; id < ND; ++id){
352  int ivcd = ivc + NVC * id;
353  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
354  }
355 
356  yt[0] = 0.5 * (vt[0] - vt[2]);
357  yt[1] = 0.5 * (vt[1] - vt[3]);
358  yt[2] = 0.5 * (vt[2] - vt[0]);
359  yt[3] = 0.5 * (vt[3] - vt[1]);
360 
361  for (int is = Ns-2; is >= 0; --is) {
362 
363  for(int id = 0; id < ND; ++id){
364  xt[id] = vt[id];
365  }
366 
367  for(int id = 0; id < ND; ++id){
368  int ivcd = ivc + NVC * id;
369  vt[id] = vp[IDX2(Nin5, (ivcd + NVCD * is), site)];
370  }
371 
372  real_t a = real_t(0.5) * dm[is + 1] * dpinv[is];
373 
374  vt[0] += a * (xt[0] + xt[2]);
375  vt[1] += a * (xt[1] + xt[3]);
376  vt[2] += a * (xt[2] + xt[0]);
377  vt[3] += a * (xt[3] + xt[1]);
378 
379  for(int id = 0; id < ND; ++id){
380  vt[id] += -e[is] * yt[id];
381  }
382 
383  for(int id = 0; id < ND; ++id){
384  int ivcd = ivc + NVC * id;
385  vp[IDX2(Nin5, (ivcd + NVCD * is), site)] = vt[id];
386  }
387  }
388 
389  }
390  } // idx loop end
391 
392  }
393  }
394 
395 }
396 
397 #endif
398 //============================================================END=====
ID1
#define ID1
Definition: fopr_Wilson_impl_SU2-inc.h:18
NVCD
#define NVCD
Definition: define_params_SU3.h:20
mult_domainwall_5din_LUdaginv_dirac
void mult_domainwall_5din_LUdaginv_dirac(real_t *RESTRICT vp, real_t *RESTRICT wp, int Ns, int *Nsize, real_t *e, real_t *f, real_t *dpinv, real_t *dm, real_t alpha)
Definition: mult_Domainwall_5din_LUinv_openacc-inc.h:204
vt4
real_t vt4
Definition: mult_Staggered_uvdn1_openacc-inc.h:8
CEIL_NWP
#define CEIL_NWP(nst)
Definition: define_params.h:47
ID2
#define ID2
Definition: fopr_Wilson_impl_SU2-inc.h:19
ID4
#define ID4
Definition: fopr_Wilson_impl_SU2-inc.h:21
vt1
real_t vt1
Definition: mult_Staggered_uvdn1_openacc-inc.h:8
AIndex_eo_accel::idx
int idx(const int in, const int Nin, const int ist, const int leo, const int Nvol2, const int ex)
Definition: aindex_eo.h:28
ND
#define ND
Definition: field_F_imp_SU2-inc.h:18
real_t
double real_t
Definition: bridgeACC_AField_double.cpp:14
ID3
#define ID3
Definition: fopr_Wilson_impl_SU2-inc.h:20
IDX2
#define IDX2(nin, in, ist)
Definition: define_index.h:28
NVC
#define NVC
Definition: fopr_Wilson_impl_SU2-inc.h:15
mult_domainwall_5din_LUinv_dirac
void mult_domainwall_5din_LUinv_dirac(real_t *RESTRICT vp, real_t *RESTRICT wp, int Ns, int *Nsize, real_t *e, real_t *f, real_t *dpinv, real_t *dm, real_t alpha)
Definition: mult_Domainwall_5din_LUinv_openacc-inc.h:14
NWP
#define NWP
Definition: define_params.h:32
vt2
real_t vt2
Definition: mult_Staggered_uvdn1_openacc-inc.h:8
vt3
real_t vt3
Definition: mult_Staggered_uvdn1_openacc-inc.h:8