CHROMA
eoprec_ovlap_contfrac5d_linop_array_opt_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Optimized version of 5D continued frac. linop
3  */
4 
5 
6 #include "chromabase.h"
7 #include "actions/ferm/linop/eoprec_ovlap_contfrac5d_linop_array_opt_w.h"
8 
9 
10 namespace Chroma
11 {
12 
13 
14 #if 0
15  // FOR SOME UNKNOWN REASON, THIS CODE IS MUCH SLOWER THAN THE QDP VERSION!!
16 
17 
18  //! Apply the even-even (odd-odd) coupling piece of the domain-wall fermion operator
19  /*!
20  * \ingroup linop
21  *
22  * The operator acts on the entire lattice
23  *
24  * \param psi Pseudofermion field (Read)
25  * \param isign Flag ( PLUS | MINUS ) (Read)
26  * \param cb checkerboard ( 0 | 1 ) (Read)
27  *
28  *
29  * Flopcount: N5*6NcNs + (N5-2)*4NcNs = NcNs( 6N5 +4(N5-2) ) = (10N5-8) Nc Ns / cb_site
30  */
31  void
32  OptEvenOddPrecOvlapContFrac5DLinOpArray::applyDiag(
33  multi1d<LatticeFermion>& chi,
34  const multi1d<LatticeFermion>& psi,
35  enum PlusMinus isign,
36  const int cb) const
37  {
38  if( chi.size() != N5 ) chi.resize(N5);
39 
40  // We don't care about the isign because our operator is Hermitian
41  // Apply matrix
42  // [ A_0 B_0 0 ... ] [ psi_0 ]
43  // [ B_0 A_1 B_1 ... ... ] [ psi_1 ]
44  // [ 0 ... ... ... ... ] [ psi_2 ]
45  // [ ... ... 0 B_N5-3 A_N5-2 B_N5-2 ] [ psi_N5-2 ]
46  // [ ... ... ... 0 B_N5-2 A_N5-1 ] [ psi_N5-1 ]
47 
48  // With A_i = gamma_5 a_i = a_i gamma_5
49  // and B_i = b_i I = alpha_i I
50 
51  // First 0ne
52  // Operation: chi[0][rb[cb]] = a[0] G5 psi[0] + alpha[0]*psi[1]
53  //
54  // Useful flops: 6Nc Ns / site
55  // tmp[rb[cb]] = Gamma(G5)*psi[0];
56  // chi[0][rb[cb]] = a[0]*tmp;
57 
58  const int start = rb[cb].start();
59  const int end = rb[cb].end();
60 
61  if( N5 > 1 ) {
62  //chi[0][rb[cb]] = alpha[0]*psi[1] + a[0]*(GammaConst<Ns,Ns*Ns-1>()*psi[0]);
63 
64  REAL scalep = alpha[0].elem().elem().elem().elem();
65  REAL scalep2 = a[0].elem().elem().elem().elem();
66  REAL *InScale = (REAL *) &(psi[1].elem(start).elem(0).elem(0).real());
67  REAL *Add = (REAL *) &(psi[0].elem(start).elem(0).elem(0).real());
68  REAL *Out = (REAL *) &(chi[0].elem(start).elem(0).elem(0).real());
69 
70  // (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*g5*std::vector)Add)
71  //void axpbyz_g5(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
72  register int index_x = 0;
73  register int index_y = 0;
74  register int index_z = 0;
75 
76  for(int site=start; site <= end; ++site)
77  {
78  // Spin Component 0 (AXPY3)
79  REAL x0r = InScale[index_x++];
80  REAL y0r = Add[index_y++];
81  REAL z0r = scalep*x0r ;
82  z0r += scalep2*y0r;
83  Out[index_z++] = z0r;
84 
85  REAL x0i = InScale[index_x++];
86  REAL y0i = Add[index_y++];
87  REAL z0i = scalep*x0i;
88  z0i += scalep2*y0i;
89  Out[index_z++] = z0i;
90 
91  REAL x1r = InScale[index_x++];
92  REAL y1r = Add[index_y++];
93  REAL z1r = scalep*x1r ;
94  z1r += scalep2*y1r;
95  Out[index_z++] = z1r;
96 
97  REAL x1i = InScale[index_x++];
98  REAL y1i = Add[index_y++];
99  REAL z1i = scalep*x1i;
100  z1i += scalep2*y1i;
101  Out[index_z++] = z1i;
102 
103  REAL x2r = InScale[index_x++];
104  REAL y2r = Add[index_y++];
105  REAL z2r = scalep*x2r ;
106  z2r += scalep2*y2r;
107  Out[index_z++] = z2r;
108 
109  REAL x2i = InScale[index_x++];
110  REAL y2i = Add[index_y++];
111  REAL z2i = scalep*x2i ;
112  z2i += scalep2*y2i;
113  Out[index_z++] = z2i;
114 
115  // Spin Component 1
116  x0r = InScale[index_x++];
117  y0r = Add[index_y++];
118  z0r = scalep*x0r;
119  z0r += scalep2*y0r;
120  Out[index_z++] = z0r;
121 
122  x0i = InScale[index_x++];
123  y0i = Add[index_y++];
124  z0i = scalep*x0i;
125  z0i += scalep2*y0i;
126  Out[index_z++] = z0i;
127 
128  x1r = InScale[index_x++];
129  y1r = Add[index_y++];
130  z1r = scalep*x1r ;
131  z1r += scalep2*y1r;
132  Out[index_z++] = z1r;
133 
134  x1i = InScale[index_x++];
135  y1i = Add[index_y++];
136  z1i = scalep*x1i;
137  z1i += scalep2*y1i;
138  Out[index_z++] = z1i;
139 
140  x2r = InScale[index_x++];
141  y2r = Add[index_y++];
142  z2r = scalep*x2r;
143  z2r += scalep2*y2r;
144  Out[index_z++] = z2r;
145 
146  x2i = InScale[index_x++];
147  y2i = Add[index_y++];
148  z2i = scalep*x2i;
149  z2i += scalep2*y2i;
150  Out[index_z++] = z2i;
151 
152  // Spin Component 2 (AXPY3)
153  x0r = InScale[index_x++];
154  y0r = Add[index_y++];
155  z0r = scalep*x0r ;
156  z0r -= scalep2*y0r;
157  Out[index_z++] = z0r;
158 
159  x0i = InScale[index_x++];
160  y0i = Add[index_y++];
161  z0i = scalep*x0i;
162  z0i -= scalep2*y0i;
163  Out[index_z++] = z0i;
164 
165  x1r = InScale[index_x++];
166  y1r = Add[index_y++];
167  z1r = scalep*x1r ;
168  z1r -= scalep2*y1r;
169  Out[index_z++] = z1r;
170 
171  x1i = InScale[index_x++];
172  y1i = Add[index_y++];
173  z1i = scalep*x1i;
174  z1i -= scalep2*y1i;
175  Out[index_z++] = z1i;
176 
177  x2r = InScale[index_x++];
178  y2r = Add[index_y++];
179  z2r = scalep*x2r ;
180  z2r -= scalep2*y2r;
181  Out[index_z++] = z2r;
182 
183  x2i = InScale[index_x++];
184  y2i = Add[index_y++];
185  z2i = scalep*x2i ;
186  z2i -= scalep2*y2i;
187  Out[index_z++] = z2i;
188 
189  // Spin Component 3
190  x0r = InScale[index_x++];
191  y0r = Add[index_y++];
192  z0r = scalep*x0r;
193  z0r -= scalep2*y0r;
194  Out[index_z++] = z0r;
195 
196  x0i = InScale[index_x++];
197  y0i = Add[index_y++];
198  z0i = scalep*x0i;
199  z0i -= scalep2*y0i;
200  Out[index_z++] = z0i;
201 
202  x1r = InScale[index_x++];
203  y1r = Add[index_y++];
204  z1r = scalep*x1r ;
205  z1r -= scalep2*y1r;
206  Out[index_z++] = z1r;
207 
208  x1i = InScale[index_x++];
209  y1i = Add[index_y++];
210  z1i = scalep*x1i;
211  z1i -= scalep2*y1i;
212  Out[index_z++] = z1i;
213 
214  x2r = InScale[index_x++];
215  y2r = Add[index_y++];
216  z2r = scalep*x2r;
217  z2r -= scalep2*y2r;
218  Out[index_z++] = z2r;
219 
220  x2i = InScale[index_x++];
221  y2i = Add[index_y++];
222  z2i = scalep*x2i;
223  z2i -= scalep2*y2i;
224  Out[index_z++] = z2i;
225  }
226  }
227  else
228  {
229  //chi[0][rb[cb]] = a[0]*(GammaConst<Ns,Ns*Ns-1>()*psi[0]);
230  REAL scalep = a[0].elem().elem().elem().elem();
231  REAL *In = (REAL *) &(psi[0].elem(start).elem(0).elem(0).real());
232  REAL *Out = (REAL *) &(chi[0].elem(start).elem(0).elem(0).real());
233 
234  // (Vector) out = (Scalar) (*scalep) * (Vector) In
235  //void scal_g5(REAL *Out, REAL *scalep, REAL *In, int n_4vec)
236  register int index_x = 0;
237  register int index_z = 0;
238 
239  for(int site=start; site <= end; ++site)
240  {
241  // Spin Component 0
242  REAL x0r = In[index_x++];
243  REAL z0r = scalep*x0r;
244  Out[index_z++] = z0r;
245 
246  REAL x0i = In[index_x++];
247  REAL z0i = scalep*x0i;
248  Out[index_z++] = z0i;
249 
250  REAL x1r = In[index_x++];
251  REAL z1r = scalep*x1r;
252  Out[index_z++] = z1r;
253 
254  REAL x1i = In[index_x++];
255  REAL z1i = scalep*x1i;
256  Out[index_z++] = z1i;
257 
258  REAL x2r = In[index_x++];
259  REAL z2r = scalep*x2r;
260  Out[index_z++] = z2r;
261 
262  REAL x2i = In[index_x++];
263  REAL z2i = scalep*x2i;
264  Out[index_z++] = z2i;
265 
266  // Spin Component 1
267  x0r = In[index_x++];
268  z0r = scalep*x0r;
269  Out[index_z++] = z0r;
270 
271  x0i = In[index_x++];
272  z0i = scalep*x0i;
273  Out[index_z++] = z0i;
274 
275  x1r = In[index_x++];
276  z1r = scalep*x1r;
277  Out[index_z++] = z1r;
278 
279  x1i = In[index_x++];
280  z1i = scalep*x1i;
281  Out[index_z++] = z1i;
282 
283  x2r = In[index_x++];
284  z2r = scalep*x2r;
285  Out[index_z++] = z2r;
286 
287  x2i = In[index_x++];
288  z2i = scalep*x2i;
289  Out[index_z++] = z2i;
290 
291  // Spin Component 2
292  x0r = In[index_x++];
293  z0r = scalep*x0r;
294  Out[index_z++] =- z0r;
295 
296  x0i = In[index_x++];
297  z0i = scalep*x0i;
298  Out[index_z++] =- z0i;
299 
300  x1r = In[index_x++];
301  z1r = scalep*x1r;
302  Out[index_z++] =-z1r;
303 
304  x1i = In[index_x++];
305  z1i = scalep*x1i;
306  Out[index_z++] =-z1i;
307 
308  x2r = In[index_x++];
309  z2r = scalep*x2r;
310  Out[index_z++] =-z2r;
311 
312  x2i = In[index_x++];
313  z2i = scalep*x2i;
314  Out[index_z++] =-z2i;
315 
316  // Spin Component 3
317  x0r = In[index_x++];
318  z0r = scalep*x0r;
319  Out[index_z++] =- z0r;
320 
321  x0i = In[index_x++];
322  z0i = scalep*x0i;
323  Out[index_z++] =- z0i;
324 
325  x1r = In[index_x++];
326  z1r = scalep*x1r;
327  Out[index_z++] =-z1r;
328 
329  x1i = In[index_x++];
330  z1i = scalep*x1i;
331  Out[index_z++] =-z1i;
332 
333  x2r = In[index_x++];
334  z2r = scalep*x2r;
335  Out[index_z++] =-z2r;
336 
337  x2i = In[index_x++];
338  z2i = scalep*x2i;
339  Out[index_z++] =-z2i;
340  }
341  }
342 
343  // All the rest
344  for(int i=1; i < N5; i++)
345  {
346  // Operation:
347  // N5 - 1 times:
348  // chi[i] = alpha[i-1]*psi[i-1] + a[i] Gamma_5 *psi[i]
349  // N5 - 2 times:
350  // chi[i] += alpha[i]*psi[i+1];
351  // Useful flops = (N5-1) * 6NcNs + (N5-2)*4Nc*Ns
352 
353  /*
354  // B_{i-1} psi_[i-1]
355  chi[i][rb[cb]] = alpha[i-1]*psi[i-1];
356 
357  // A_{i} psi[i] = a_{i} g_5 psi[i]
358  tmp[rb[cb]] = Gamma(G5)*psi[i];
359  chi[i][rb[cb]] += a[i]*tmp;
360  */
361  {
362  //chi[i][rb[cb]] = alpha[i-1]*psi[i-1] + a[i]*(GammaConst<Ns,Ns*Ns-1>()*psi[i]);
363 
364  REAL scalep = alpha[i-1].elem().elem().elem().elem();
365  REAL scalep2 = a[i].elem().elem().elem().elem();
366  REAL *InScale = (REAL *) &(psi[i-1].elem(start).elem(0).elem(0).real());
367  REAL *Add = (REAL *) &(psi[i].elem(start).elem(0).elem(0).real());
368  REAL *Out = (REAL *) &(chi[i].elem(start).elem(0).elem(0).real());
369 
370  // (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*g5*std::vector)Add)
371  //void axpbyz_g5(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
372  register int index_x = 0;
373  register int index_y = 0;
374  register int index_z = 0;
375 
376  for(int site=start; site <= end; ++site)
377  {
378  // Spin Component 0 (AXPY3)
379  REAL x0r = InScale[index_x++];
380  REAL y0r = Add[index_y++];
381  REAL z0r = scalep*x0r ;
382  z0r += scalep2*y0r;
383  Out[index_z++] = z0r;
384 
385  REAL x0i = InScale[index_x++];
386  REAL y0i = Add[index_y++];
387  REAL z0i = scalep*x0i;
388  z0i += scalep2*y0i;
389  Out[index_z++] = z0i;
390 
391  REAL x1r = InScale[index_x++];
392  REAL y1r = Add[index_y++];
393  REAL z1r = scalep*x1r ;
394  z1r += scalep2*y1r;
395  Out[index_z++] = z1r;
396 
397  REAL x1i = InScale[index_x++];
398  REAL y1i = Add[index_y++];
399  REAL z1i = scalep*x1i;
400  z1i += scalep2*y1i;
401  Out[index_z++] = z1i;
402 
403  REAL x2r = InScale[index_x++];
404  REAL y2r = Add[index_y++];
405  REAL z2r = scalep*x2r ;
406  z2r += scalep2*y2r;
407  Out[index_z++] = z2r;
408 
409  REAL x2i = InScale[index_x++];
410  REAL y2i = Add[index_y++];
411  REAL z2i = scalep*x2i ;
412  z2i += scalep2*y2i;
413  Out[index_z++] = z2i;
414 
415  // Spin Component 1
416  x0r = InScale[index_x++];
417  y0r = Add[index_y++];
418  z0r = scalep*x0r;
419  z0r += scalep2*y0r;
420  Out[index_z++] = z0r;
421 
422  x0i = InScale[index_x++];
423  y0i = Add[index_y++];
424  z0i = scalep*x0i;
425  z0i += scalep2*y0i;
426  Out[index_z++] = z0i;
427 
428  x1r = InScale[index_x++];
429  y1r = Add[index_y++];
430  z1r = scalep*x1r ;
431  z1r += scalep2*y1r;
432  Out[index_z++] = z1r;
433 
434  x1i = InScale[index_x++];
435  y1i = Add[index_y++];
436  z1i = scalep*x1i;
437  z1i += scalep2*y1i;
438  Out[index_z++] = z1i;
439 
440  x2r = InScale[index_x++];
441  y2r = Add[index_y++];
442  z2r = scalep*x2r;
443  z2r += scalep2*y2r;
444  Out[index_z++] = z2r;
445 
446  x2i = InScale[index_x++];
447  y2i = Add[index_y++];
448  z2i = scalep*x2i;
449  z2i += scalep2*y2i;
450  Out[index_z++] = z2i;
451 
452  // Spin Component 2 (AXPY3)
453  x0r = InScale[index_x++];
454  y0r = Add[index_y++];
455  z0r = scalep*x0r ;
456  z0r -= scalep2*y0r;
457  Out[index_z++] = z0r;
458 
459  x0i = InScale[index_x++];
460  y0i = Add[index_y++];
461  z0i = scalep*x0i;
462  z0i -= scalep2*y0i;
463  Out[index_z++] = z0i;
464 
465  x1r = InScale[index_x++];
466  y1r = Add[index_y++];
467  z1r = scalep*x1r ;
468  z1r -= scalep2*y1r;
469  Out[index_z++] = z1r;
470 
471  x1i = InScale[index_x++];
472  y1i = Add[index_y++];
473  z1i = scalep*x1i;
474  z1i -= scalep2*y1i;
475  Out[index_z++] = z1i;
476 
477  x2r = InScale[index_x++];
478  y2r = Add[index_y++];
479  z2r = scalep*x2r ;
480  z2r -= scalep2*y2r;
481  Out[index_z++] = z2r;
482 
483  x2i = InScale[index_x++];
484  y2i = Add[index_y++];
485  z2i = scalep*x2i ;
486  z2i -= scalep2*y2i;
487  Out[index_z++] = z2i;
488 
489  // Spin Component 3
490  x0r = InScale[index_x++];
491  y0r = Add[index_y++];
492  z0r = scalep*x0r;
493  z0r -= scalep2*y0r;
494  Out[index_z++] = z0r;
495 
496  x0i = InScale[index_x++];
497  y0i = Add[index_y++];
498  z0i = scalep*x0i;
499  z0i -= scalep2*y0i;
500  Out[index_z++] = z0i;
501 
502  x1r = InScale[index_x++];
503  y1r = Add[index_y++];
504  z1r = scalep*x1r ;
505  z1r -= scalep2*y1r;
506  Out[index_z++] = z1r;
507 
508  x1i = InScale[index_x++];
509  y1i = Add[index_y++];
510  z1i = scalep*x1i;
511  z1i -= scalep2*y1i;
512  Out[index_z++] = z1i;
513 
514  x2r = InScale[index_x++];
515  y2r = Add[index_y++];
516  z2r = scalep*x2r;
517  z2r -= scalep2*y2r;
518  Out[index_z++] = z2r;
519 
520  x2i = InScale[index_x++];
521  y2i = Add[index_y++];
522  z2i = scalep*x2i;
523  z2i -= scalep2*y2i;
524  Out[index_z++] = z2i;
525  }
526  }
527 
528  // When i hits N5-1, we don't have the B_N5-1 term
529  if(i < N5-1)
530  {
531  //chi[i][rb[cb]] += alpha[i]*psi[i+1];
532 
533  REAL scalep = alpha[i].elem().elem().elem().elem();
534  REAL* InScale = (REAL *)&(psi[i+1].elem(start).elem(0).elem(0).real());
535  REAL* Add = (REAL *)&(chi[i].elem(start).elem(0).elem(0).real());
536  REAL* Out = (REAL *)&(chi[i].elem(start).elem(0).elem(0).real());
537 
538  //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
539  // (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (Vector) Add
540  //void vaxpy3(REAL *Out,REAL *scalep,REAL *InScale, REAL *Add,int n_3vec)
541  register int index_x = 0;
542  register int index_y = 0;
543  register int index_z = 0;
544 
545  for(int site=start; site <= end; ++site)
546  {
547  // Spin Component 0 (AXPY3)
548  REAL x0r = InScale[index_x++];
549  REAL y0r = Add[index_y++];
550  REAL z0r = scalep*x0r + y0r;
551  Out[index_z++] =(REAL) z0r;
552 
553  REAL x0i = InScale[index_x++];
554  REAL y0i = Add[index_y++];
555  REAL z0i = scalep*x0i + y0i;
556  Out[index_z++] =(REAL) z0i;
557 
558  REAL x1r = InScale[index_x++];
559  REAL y1r = Add[index_y++];
560  REAL z1r = scalep*x1r + y1r;
561  Out[index_z++] = (REAL)z1r;
562 
563  REAL x1i = InScale[index_x++];
564  REAL y1i = Add[index_y++];
565  REAL z1i = scalep*x1i + y1i;
566  Out[index_z++] = (REAL)z1i;
567 
568  REAL x2r = InScale[index_x++];
569  REAL y2r = Add[index_y++];
570  REAL z2r = scalep*x2r + y2r;
571  Out[index_z++] = (REAL)z2r;
572 
573  REAL x2i = InScale[index_x++];
574  REAL y2i = Add[index_y++];
575  REAL z2i = scalep*x2i + y2i;
576  Out[index_z++] = (REAL)z2i;
577 
578  // Spin Component 1 (AXPY3)
579  x0r = InScale[index_x++];
580  y0r = Add[index_y++];
581  z0r = scalep*x0r + y0r;
582  Out[index_z++] =(REAL) z0r;
583 
584  x0i = InScale[index_x++];
585  y0i = Add[index_y++];
586  z0i = scalep*x0i + y0i;
587  Out[index_z++] =(REAL) z0i;
588 
589  x1r = InScale[index_x++];
590  y1r = Add[index_y++];
591  z1r = scalep*x1r + y1r;
592  Out[index_z++] = (REAL)z1r;
593 
594  x1i = InScale[index_x++];
595  y1i = Add[index_y++];
596  z1i = scalep*x1i + y1i;
597  Out[index_z++] = (REAL)z1i;
598 
599  x2r = InScale[index_x++];
600  y2r = Add[index_y++];
601  z2r = scalep*x2r + y2r;
602  Out[index_z++] = (REAL)z2r;
603 
604  x2i = InScale[index_x++];
605  y2i = Add[index_y++];
606  z2i = scalep*x2i + y2i;
607  Out[index_z++] = (REAL)z2i;
608 
609  // Spin Component 2 (AXPY3)
610  x0r = InScale[index_x++];
611  y0r = Add[index_y++];
612  z0r = scalep*x0r + y0r;
613  Out[index_z++] =(REAL) z0r;
614 
615  x0i = InScale[index_x++];
616  y0i = Add[index_y++];
617  z0i = scalep*x0i + y0i;
618  Out[index_z++] =(REAL) z0i;
619 
620  x1r = InScale[index_x++];
621  y1r = Add[index_y++];
622  z1r = scalep*x1r + y1r;
623  Out[index_z++] = (REAL)z1r;
624 
625  x1i = InScale[index_x++];
626  y1i = Add[index_y++];
627  z1i = scalep*x1i + y1i;
628  Out[index_z++] = (REAL)z1i;
629 
630  x2r = InScale[index_x++];
631  y2r = Add[index_y++];
632  z2r = scalep*x2r + y2r;
633  Out[index_z++] = (REAL)z2r;
634 
635  x2i = InScale[index_x++];
636  y2i = Add[index_y++];
637  z2i = scalep*x2i + y2i;
638  Out[index_z++] = (REAL)z2i;
639 
640  // Spin Component 3 (AXPY3)
641  x0r = InScale[index_x++];
642  y0r = Add[index_y++];
643  z0r = scalep*x0r + y0r;
644  Out[index_z++] =(REAL) z0r;
645 
646  x0i = InScale[index_x++];
647  y0i = Add[index_y++];
648  z0i = scalep*x0i + y0i;
649  Out[index_z++] =(REAL) z0i;
650 
651  x1r = InScale[index_x++];
652  y1r = Add[index_y++];
653  z1r = scalep*x1r + y1r;
654  Out[index_z++] = (REAL)z1r;
655 
656  x1i = InScale[index_x++];
657  y1i = Add[index_y++];
658  z1i = scalep*x1i + y1i;
659  Out[index_z++] = (REAL)z1i;
660 
661  x2r = InScale[index_x++];
662  y2r = Add[index_y++];
663  z2r = scalep*x2r + y2r;
664  Out[index_z++] = (REAL)z2r;
665 
666  x2i = InScale[index_x++];
667  y2i = Add[index_y++];
668  z2i = scalep*x2i + y2i;
669  Out[index_z++] = (REAL)z2i;
670  }
671  }
672  }
673  }
674 #endif
675 
676 
677  //! Apply the inverse even-even (odd-odd)
678  /*!
679  * \ingroup linop
680  *
681  * Here we apply the LDU decomposition
682  *
683  * \param psi Pseudofermion field (Read)
684  * \param isign Flag ( PLUS | MINUS ) (Read)
685  * \param cb checkerboard ( 0 | 1 ) (Read)
686 
687  *
688  * Total flopcount: (N5-1)*4NcNs + 2NcNs + (N5-1)*6NcNs
689  * = (N5-1)*10NcNs + 2NcNs
690  * = (10N5-8) Nc Ns per_cb_site
691  */
692  void
693  OptEvenOddPrecOvlapContFrac5DLinOpArray::applyDiagInv(
694  multi1d<LatticeFermion>& chi,
695  const multi1d<LatticeFermion>& psi,
696  enum PlusMinus isign,
697  const int cb) const
698  {
699  if( chi.size() != N5 ) chi.resize(N5);
700 
701  multi1d<Fermion> y(N5);
702 
703  for(int site=rb[cb].start(); site <= rb[cb].end(); ++site)
704  {
705  // Solve L y = psi
706  {
707  //y[0] = psi[0];
708  REAL *In = (REAL *) &(psi[0].elem(site).elem(0).elem(0).real());
709  REAL* Out = (REAL *) &(y[0].elem().elem(0).elem(0).real());
710 
711  register int index_x = 0;
712  register int index_z = 0;
713 
714  Out[index_z++] = In[index_x++]; // 00r
715  Out[index_z++] = In[index_x++]; // 00i
716  Out[index_z++] = In[index_x++]; // 01r
717  Out[index_z++] = In[index_x++]; // 01i
718  Out[index_z++] = In[index_x++]; // 02r
719  Out[index_z++] = In[index_x++]; // 02i
720 
721  Out[index_z++] = In[index_x++]; // 10r
722  Out[index_z++] = In[index_x++]; // 10i
723  Out[index_z++] = In[index_x++]; // 11r
724  Out[index_z++] = In[index_x++]; // 11i
725  Out[index_z++] = In[index_x++]; // 12r
726  Out[index_z++] = In[index_x++]; // 12i
727 
728  Out[index_z++] = In[index_x++]; // 20r
729  Out[index_z++] = In[index_x++]; // 20i
730  Out[index_z++] = In[index_x++]; // 21r
731  Out[index_z++] = In[index_x++]; // 21i
732  Out[index_z++] = In[index_x++]; // 22r
733  Out[index_z++] = In[index_x++]; // 22i
734 
735  Out[index_z++] = In[index_x++]; // 30r
736  Out[index_z++] = In[index_x++]; // 30i
737  Out[index_z++] = In[index_x++]; // 31r
738  Out[index_z++] = In[index_x++]; // 31i
739  Out[index_z++] = In[index_x++]; // 32r
740  Out[index_z++] = In[index_x++]; // 32i
741  }
742 
743  // (N5-1)*4NcNs
744  for(int i = 1; i < N5; i++)
745  {
746  // tmp[rb[cb]] = Gamma(G5)*y[i-1];
747  // y[i][rb[cb]] = psi[i] - u[i-1]*tmp;
748  //y[i][rb[cb]] = psi[i] - u[i-1]*(GammaConst<Ns,Ns*Ns-1>()*y[i-1]);
749 
750  REAL scalep = u[i-1].elem().elem().elem().elem();
751  REAL *InScale = (REAL *) &(y[i-1].elem().elem(0).elem(0).real());
752  REAL *Add = (REAL *) &(psi[i].elem(site).elem(0).elem(0).real());
753  REAL* Out = (REAL *) &(y[i].elem().elem(0).elem(0).real());
754 
755  // (Vector) out = (Vector) Add + (Scalar) (*scalep) * (Vector) P{+} InScale
756  // void xmayz_g5(REAL *Out,REAL *scalep,REAL *Add, REAL *InScale,int n_4vec)
757  register int index_x = 0;
758  register int index_y = 0;
759  register int index_z = 0;
760 
761  // Spin Component 0 (AYPX)
762  REAL x0r = Add[index_x++];
763  REAL y0r = InScale[index_y++];
764  REAL z0r = x0r - scalep*y0r;
765  Out[index_z++] = z0r;
766 
767  REAL x0i = Add[index_x++];
768  REAL y0i = InScale[index_y++];
769  REAL z0i = x0i - scalep*y0i;
770  Out[index_z++] = z0i;
771 
772  REAL x1r = Add[index_x++];
773  REAL y1r = InScale[index_y++];
774  REAL z1r = x1r - scalep*y1r;
775  Out[index_z++] = z1r;
776 
777  REAL x1i = Add[index_x++];
778  REAL y1i = InScale[index_y++];
779  REAL z1i = x1i - scalep*y1i;
780  Out[index_z++] = z1i;
781 
782  REAL x2r = Add[index_x++];
783  REAL y2r = InScale[index_y++];
784  REAL z2r = x2r - scalep*y2r;
785  Out[index_z++] = z2r;
786 
787  REAL x2i = Add[index_x++];
788  REAL y2i = InScale[index_y++];
789  REAL z2i = x2i - scalep*y2i;
790  Out[index_z++] = z2i;
791 
792  // Spin Component 1 (AYPX)
793  x0r = Add[index_x++];
794  y0r = InScale[index_y++];
795  z0r = x0r - scalep*y0r;
796  Out[index_z++] = z0r;
797 
798  x0i = Add[index_x++];
799  y0i = InScale[index_y++];
800  z0i = x0i - scalep*y0i;
801  Out[index_z++] = z0i;
802 
803  x1r = Add[index_x++];
804  y1r = InScale[index_y++];
805  z1r = x1r - scalep*y1r;
806  Out[index_z++] = z1r;
807 
808  x1i = Add[index_x++];
809  y1i = InScale[index_y++];
810  z1i = x1i - scalep*y1i;
811  Out[index_z++] = z1i;
812 
813  x2r = Add[index_x++];
814  y2r = InScale[index_y++];
815  z2r = x2r - scalep*y2r;
816  Out[index_z++] = z2r;
817 
818  x2i = Add[index_x++];
819  y2i = InScale[index_y++];
820  z2i = x2i - scalep*y2i;
821  Out[index_z++] = z2i;
822 
823  // Spin Component 2 (AYPX)
824  x0r = Add[index_x++];
825  y0r = InScale[index_y++];
826  z0r = x0r + scalep*y0r;
827  Out[index_z++] = z0r;
828 
829  x0i = Add[index_x++];
830  y0i = InScale[index_y++];
831  z0i = x0i + scalep*y0i;
832  Out[index_z++] = z0i;
833 
834  x1r = Add[index_x++];
835  y1r = InScale[index_y++];
836  z1r = x1r + scalep*y1r;
837  Out[index_z++] = z1r;
838 
839  x1i = Add[index_x++];
840  y1i = InScale[index_y++];
841  z1i = x1i + scalep*y1i;
842  Out[index_z++] = z1i;
843 
844  x2r = Add[index_x++];
845  y2r = InScale[index_y++];
846  z2r = x2r + scalep*y2r;
847  Out[index_z++] = z2r;
848 
849  x2i = Add[index_x++];
850  y2i = InScale[index_y++];
851  z2i = x2i + scalep*y2i;
852  Out[index_z++] = z2i;
853 
854  // Spin Component 3 (AYPX)
855  x0r = Add[index_x++];
856  y0r = InScale[index_y++];
857  z0r = x0r + scalep*y0r;
858  Out[index_z++] = z0r;
859 
860  x0i = Add[index_x++];
861  y0i = InScale[index_y++];
862  z0i = x0i + scalep*y0i;
863  Out[index_z++] = z0i;
864 
865  x1r = Add[index_x++];
866  y1r = InScale[index_y++];
867  z1r = x1r + scalep*y1r;
868  Out[index_z++] = z1r;
869 
870  x1i = Add[index_x++];
871  y1i = InScale[index_y++];
872  z1i = x1i + scalep*y1i;
873  Out[index_z++] = z1i;
874 
875  x2r = Add[index_x++];
876  y2r = InScale[index_y++];
877  z2r = x2r + scalep*y2r;
878  Out[index_z++] = z2r;
879 
880  x2i = Add[index_x++];
881  y2i = InScale[index_y++];
882  z2i = x2i + scalep*y2i;
883  Out[index_z++] = z2i;
884  }
885 
886  // Invert diagonal piece y <- D^{-1} y
887  // N5 times: y = (1/d_i) gamma_5 y[i]
888  // Useful flops: N5 * 2NcNs flops / site
889  // Rolled this into the bottom loop
890 
891  // Backsubstitute U chi = y
892 
893  // 2NcNs
894 
895  {
896  //chi[N5-1][rb[cb]] = invd[N5-1]*(GammaConst<Ns,Ns*Ns-1>()*y[N5-1]);
897  REAL scalep = invd[N5-1].elem().elem().elem().elem();
898  REAL *In = (REAL *) &(y[N5-1].elem().elem(0).elem(0).real());
899  REAL *Out = (REAL *) &(chi[N5-1].elem(site).elem(0).elem(0).real());
900 
901  // (Vector) out = (Scalar) (*scalep) * (Vector) In
902  //void scal_g5(REAL *Out, REAL *scalep, REAL *In, int n_4vec)
903  register int index_x = 0;
904  register int index_z = 0;
905 
906  // Spin Component 0
907  REAL x0r = In[index_x++];
908  REAL z0r = scalep*x0r;
909  Out[index_z++] = z0r;
910 
911  REAL x0i = In[index_x++];
912  REAL z0i = scalep*x0i;
913  Out[index_z++] = z0i;
914 
915  REAL x1r = In[index_x++];
916  REAL z1r = scalep*x1r;
917  Out[index_z++] = z1r;
918 
919  REAL x1i = In[index_x++];
920  REAL z1i = scalep*x1i;
921  Out[index_z++] = z1i;
922 
923  REAL x2r = In[index_x++];
924  REAL z2r = scalep*x2r;
925  Out[index_z++] = z2r;
926 
927  REAL x2i = In[index_x++];
928  REAL z2i = scalep*x2i;
929  Out[index_z++] = z2i;
930 
931  // Spin Component 1
932  x0r = In[index_x++];
933  z0r = scalep*x0r;
934  Out[index_z++] = z0r;
935 
936  x0i = In[index_x++];
937  z0i = scalep*x0i;
938  Out[index_z++] = z0i;
939 
940  x1r = In[index_x++];
941  z1r = scalep*x1r;
942  Out[index_z++] = z1r;
943 
944  x1i = In[index_x++];
945  z1i = scalep*x1i;
946  Out[index_z++] = z1i;
947 
948  x2r = In[index_x++];
949  z2r = scalep*x2r;
950  Out[index_z++] = z2r;
951 
952  x2i = In[index_x++];
953  z2i = scalep*x2i;
954  Out[index_z++] = z2i;
955 
956  // Spin Component 2
957  x0r = In[index_x++];
958  z0r = scalep*x0r;
959  Out[index_z++] =- z0r;
960 
961  x0i = In[index_x++];
962  z0i = scalep*x0i;
963  Out[index_z++] =- z0i;
964 
965  x1r = In[index_x++];
966  z1r = scalep*x1r;
967  Out[index_z++] =-z1r;
968 
969  x1i = In[index_x++];
970  z1i = scalep*x1i;
971  Out[index_z++] =-z1i;
972 
973  x2r = In[index_x++];
974  z2r = scalep*x2r;
975  Out[index_z++] =-z2r;
976 
977  x2i = In[index_x++];
978  z2i = scalep*x2i;
979  Out[index_z++] =-z2i;
980 
981  // Spin Component 3
982  x0r = In[index_x++];
983  z0r = scalep*x0r;
984  Out[index_z++] =- z0r;
985 
986  x0i = In[index_x++];
987  z0i = scalep*x0i;
988  Out[index_z++] =- z0i;
989 
990  x1r = In[index_x++];
991  z1r = scalep*x1r;
992  Out[index_z++] =-z1r;
993 
994  x1i = In[index_x++];
995  z1i = scalep*x1i;
996  Out[index_z++] =-z1i;
997 
998  x2r = In[index_x++];
999  z2r = scalep*x2r;
1000  Out[index_z++] =-z2r;
1001 
1002  x2i = In[index_x++];
1003  z2i = scalep*x2i;
1004  Out[index_z++] =-z2i;
1005  }
1006 
1007  // N5-1 * 6NcNs
1008  for(int i = N5-2; i >= 0; i--)
1009  {
1010  // tmp[rb[cb]] = Gamma(G5)*chi[i+1]
1011  // chi[i][rb[cb]] = y[i] - u[i]*tmp;
1012  // y[i][rb[cb]] = invd[i]*(GammaConst<Ns,Ns*Ns-1>()*y[i]);
1013  // chi[i][rb[cb]] = y[i] - u[i]*(GammaConst<Ns,Ns*Ns-1>()*chi[i+1]);
1014 
1015  //chi[i][rb[cb]] = GammaConst<Ns,Ns*Ns-1>()*(invd[i]*y[i]-u[i]*chi[i+1]);
1016  REAL scalep = invd[i].elem().elem().elem().elem();
1017  REAL scalep2 = u[i].elem().elem().elem().elem();
1018  REAL *InScale = (REAL *) &(y[i].elem().elem(0).elem(0).real());
1019  REAL *Add = (REAL *) &(chi[i+1].elem(site).elem(0).elem(0).real());
1020  REAL* Out = (REAL *) &(chi[i].elem(site).elem(0).elem(0).real());
1021 
1022  // (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*g5*std::vector)Add)
1023  //g5_axmbyz(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
1024  register int index_x = 0;
1025  register int index_y = 0;
1026  register int index_z = 0;
1027 
1028  // Spin Component 0 (AXPY3)
1029  REAL x0r = InScale[index_x++];
1030  REAL y0r = Add[index_y++];
1031  REAL z0r = scalep*x0r ;
1032  z0r -= scalep2*y0r;
1033  Out[index_z++] = z0r;
1034 
1035  REAL x0i = InScale[index_x++];
1036  REAL y0i = Add[index_y++];
1037  REAL z0i = scalep*x0i;
1038  z0i -= scalep2*y0i;
1039  Out[index_z++] = z0i;
1040 
1041  REAL x1r = InScale[index_x++];
1042  REAL y1r = Add[index_y++];
1043  REAL z1r = scalep*x1r ;
1044  z1r -= scalep2*y1r;
1045  Out[index_z++] = z1r;
1046 
1047  REAL x1i = InScale[index_x++];
1048  REAL y1i = Add[index_y++];
1049  REAL z1i = scalep*x1i;
1050  z1i -= scalep2*y1i;
1051  Out[index_z++] = z1i;
1052 
1053  REAL x2r = InScale[index_x++];
1054  REAL y2r = Add[index_y++];
1055  REAL z2r = scalep*x2r ;
1056  z2r -= scalep2*y2r;
1057  Out[index_z++] = z2r;
1058 
1059  REAL x2i = InScale[index_x++];
1060  REAL y2i = Add[index_y++];
1061  REAL z2i = scalep*x2i ;
1062  z2i -= scalep2*y2i;
1063  Out[index_z++] = z2i;
1064 
1065  // Spin Component 1
1066  x0r = InScale[index_x++];
1067  y0r = Add[index_y++];
1068  z0r = scalep*x0r;
1069  z0r -= scalep2*y0r;
1070  Out[index_z++] = z0r;
1071 
1072  x0i = InScale[index_x++];
1073  y0i = Add[index_y++];
1074  z0i = scalep*x0i;
1075  z0i -= scalep2*y0i;
1076  Out[index_z++] = z0i;
1077 
1078  x1r = InScale[index_x++];
1079  y1r = Add[index_y++];
1080  z1r = scalep*x1r ;
1081  z1r -= scalep2*y1r;
1082  Out[index_z++] = z1r;
1083 
1084  x1i = InScale[index_x++];
1085  y1i = Add[index_y++];
1086  z1i = scalep*x1i;
1087  z1i -= scalep2*y1i;
1088  Out[index_z++] = z1i;
1089 
1090  x2r = InScale[index_x++];
1091  y2r = Add[index_y++];
1092  z2r = scalep*x2r;
1093  z2r -= scalep2*y2r;
1094  Out[index_z++] = z2r;
1095 
1096  x2i = InScale[index_x++];
1097  y2i = Add[index_y++];
1098  z2i = scalep*x2i;
1099  z2i -= scalep2*y2i;
1100  Out[index_z++] = z2i;
1101 
1102  // Spin Component 2 (AXPY3)
1103  x0r = InScale[index_x++];
1104  y0r = Add[index_y++];
1105  z0r = scalep2*y0r;
1106  z0r -= scalep*x0r ;
1107 
1108  Out[index_z++] = z0r;
1109 
1110  x0i = InScale[index_x++];
1111  y0i = Add[index_y++];
1112  z0i = scalep2*y0i;
1113  z0i -= scalep*x0i;
1114  Out[index_z++] = z0i;
1115 
1116  x1r = InScale[index_x++];
1117  y1r = Add[index_y++];
1118  z1r = scalep2*y1r;
1119  z1r -= scalep*x1r ;
1120  Out[index_z++] = z1r;
1121 
1122  x1i = InScale[index_x++];
1123  y1i = Add[index_y++];
1124  z1i = scalep2*y1i;
1125  z1i -= scalep*x1i;
1126  Out[index_z++] = z1i;
1127 
1128  x2r = InScale[index_x++];
1129  y2r = Add[index_y++];
1130  z2r = scalep2*y2r;
1131  z2r -= scalep*x2r ;
1132  Out[index_z++] = z2r;
1133 
1134  x2i = InScale[index_x++];
1135  y2i = Add[index_y++];
1136  z2i = scalep2*y2i;
1137  z2i -= scalep*x2i ;
1138  Out[index_z++] = z2i;
1139 
1140  // Spin Component 3 (AXPY3)
1141  x0r = InScale[index_x++];
1142  y0r = Add[index_y++];
1143  z0r = scalep2*y0r;
1144  z0r -= scalep*x0r ;
1145 
1146  Out[index_z++] = z0r;
1147 
1148  x0i = InScale[index_x++];
1149  y0i = Add[index_y++];
1150  z0i = scalep2*y0i;
1151  z0i -= scalep*x0i;
1152  Out[index_z++] = z0i;
1153 
1154  x1r = InScale[index_x++];
1155  y1r = Add[index_y++];
1156  z1r = scalep2*y1r;
1157  z1r -= scalep*x1r ;
1158  Out[index_z++] = z1r;
1159 
1160  x1i = InScale[index_x++];
1161  y1i = Add[index_y++];
1162  z1i = scalep2*y1i;
1163  z1i -= scalep*x1i;
1164  Out[index_z++] = z1i;
1165 
1166  x2r = InScale[index_x++];
1167  y2r = Add[index_y++];
1168  z2r = scalep2*y2r;
1169  z2r -= scalep*x2r ;
1170  Out[index_z++] = z2r;
1171 
1172  x2i = InScale[index_x++];
1173  y2i = Add[index_y++];
1174  z2i = scalep2*y2i;
1175  z2i -= scalep*x2i ;
1176  Out[index_z++] = z2i;
1177  }
1178  }
1179  }
1180 
1181 } // End Namespace Chroma
1182 
chromabase.h
Primary include file for CHROMA library code.
eoprec_ovlap_contfrac5d_linop_array_opt_w.h
Optimized Even-odd prec. 5D continued fraction linop.
Chroma::OptEvenOddPrecOvlapContFrac5DLinOpArray::applyDiagInv
void applyDiagInv(multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign, const int cb) const
Apply the inverse even-even (odd-odd) coupling piece of the domain-wall fermion operator.
Definition: eoprec_ovlap_contfrac5d_linop_array_opt_w.cc:693
Chroma::EvenOddPrecOvlapContFrac5DLinOpBaseArray::applyDiag
virtual void applyDiag(multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign, const int cb) const
Apply the even-even (odd-odd) coupling piece of the domain-wall fermion operator.
Definition: eoprec_ovlap_contfrac5d_linop_base_array_w.cc:135
y
int y
Definition: meslate.cc:35
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::chi
multi1d< LatticeFermion > chi(Ncb)
Chroma::psi
LatticeFermion psi
Definition: mespbg5p_w.cc:35
Chroma::cb
int cb
Definition: invbicg.cc:120
Generated by doxygen 1.9.1