CHROMA
eoprec_nef_general_linop_array_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief 4D-style even-odd preconditioned NEF domain-wall linear operator
3  */
4 
5 #include "chromabase.h"
6 #include "actions/ferm/linop/eoprec_nef_general_linop_array_w.h"
7 
8 using namespace QDP::Hints;
9 
10 namespace Chroma
11 {
12 
13  //! Creation routine
14  /*! \ingroup fermact
15  *
16  * \param u_ gauge field (Read)
17  * \param WilsonMass_ DWF height (Read)
18  * \param b5_ NEF parameter array (Read)
19  * \param c5_ NEF parameter array (Read)
20  * \param m_q_ quark mass (Read)
21  * \param N5_ extent of 5D (Read)
22  */
23  void
24  EvenOddPrecGenNEFDWLinOpArray::create(Handle< FermState<T,P,Q> > fs,
25  const Real& WilsonMass_,
26  const Real& m_q_,
27  const multi1d<Real>& b5_,
28  const multi1d<Real>& c5_,
29  int N5_)
30  {
31  START_CODE();
32 
33  WilsonMass = WilsonMass_;
34  m_q = m_q_;
35 
36  // Sanity checking:
37  if( b5_.size() != N5_ ) {
38  QDPIO::cerr << "b5 array size and N5 are inconsistent" << std::endl;
39  QDPIO::cerr << "b5_array.size() = " << b5_.size() << std::endl;
40  QDPIO::cerr << "N5_ = " << N5_ << std::endl << std::flush;
41  QDP_abort(1);
42  }
43  N5 = N5_;
44  D.create(fs, N5);
45 
46  b5.resize(N5);
47  c5.resize(N5);
48  for(int i=0; i < N5; i++) {
49  b5[i] = b5_[i];
50  c5[i] = c5_[i];
51  }
52 
53  f_plus.resize(N5);
54  f_minus.resize(N5);
55  for(int i=0; i < N5; i++) {
56  f_plus[i] = b5[i]*( Real(Nd) - WilsonMass ) + 1;
57  f_minus[i]= c5[i]*( Real(Nd) - WilsonMass ) - 1;
58  }
59 
60 
61  l.resize(N5-1);
62  r.resize(N5-1);
63 
64  l[0] = -m_q*f_minus[N5-1]/f_plus[0];
65  r[0] = -m_q*f_minus[0]/f_plus[0];
66 
67  for(int i=1; i < N5-1; i++) {
68  l[i] = -(f_minus[i-1]/f_plus[i])*l[i-1];
69  r[i] = -(f_minus[i]/f_plus[i])*r[i-1];
70  }
71 
72  a.resize(N5-1);
73  b.resize(N5-1);
74  for(int i=0; i < N5-1; i++) {
75  a[i] = f_minus[i+1]/f_plus[i];
76  b[i] = f_minus[i]/f_plus[i];
77  }
78 
79  d.resize(N5);
80  for(int i=0; i < N5; i++) {
81  d[i] = f_plus[i];
82  }
83  // Last bits of d can be computed 2 ways which should be equal
84  Real tmp1 = f_minus[N5-2]*l[N5-2];
85  Real tmp2 = f_minus[N5-1]*r[N5-2];
86 
87  // Sanity checking:
88  // QDPIO::cout << " The following 2 should be equal: ";
89  // QDPIO::cout << tmp1 << " and " << tmp2 << std::endl;
90 
91 
92  // Subtrace ONLY ONE of them onto d[N5-1]
93  d[N5-1] -= tmp1;
94 
95  END_CODE();
96  }
97 
98 
99  //! Apply the even-even (odd-odd) coupling piece of the domain-wall fermion operator
100  /*!
101  * \ingroup linop
102  *
103  * The operator acts on the entire lattice
104  *
105  * \param psi Pseudofermion field (Read)
106  * \param isign Flag ( PLUS | MINUS ) (Read)
107  * \param cb checkerboard ( 0 | 1 ) (Read)
108  */
109  void
110  EvenOddPrecGenNEFDWLinOpArray::applyDiag(multi1d<LatticeFermion>& chi,
111  const multi1d<LatticeFermion>& psi,
112  enum PlusMinus isign,
113  int cb) const
114  {
115  START_CODE();
116 
117  if( chi.size() != N5 ) chi.resize(N5);
118 
119  switch ( isign ) {
120 
121  case PLUS:
122  {
123  Real fact;
124 
125 
126  // fplus[0]*psi[0] + fminus[0]P_- psi[1] - m fminus[0] P_+ psi[N5-1]
127  // = fplus[0]*psi[0] + (fminus/2)[ (1-g5)psi[1] - m(1+g5) psi[N5-1] ]
128  // = fplus[0]*psi[0] + (fminus/2)[ psi[1] - m psi[N5-1]
129  // - g5( psi[1] + m psi[N5-1] ) ]
130  // = fplus[0]*psi[0] + (fminus/2)[ tmp_1 - g5 tmp2 ]
131  //
132  // with tmp1 = psi[1] - m psi[N5-1]
133  // tmp2 = psi[1] + m psi[N5-1]
134 
135  /* Old Code */
136  /*
137  LatticeFermion tmp1, tmp2;
138 
139  tmp1[rb[cb]] = psi[1] - m_q*psi[N5-1];
140  tmp2[rb[cb]] = psi[1] + m_q*psi[N5-1];
141  fact = Real(0.5)*f_minus[0];
142 
143  chi[0][rb[cb]] = f_plus[0]*psi[0] +
144  fact*( tmp1 - GammaConst<Ns, Ns*Ns-1>()*tmp2 );
145  */
146 
147  // Recoded using chiral projectors
148  fact = f_minus[0]*m_q;
149  chi[0][rb[cb]] = f_plus[0]*psi[0];
150  chi[0][rb[cb]] += f_minus[0]*chiralProjectMinus(psi[1]);
151  chi[0][rb[cb]] -= fact*chiralProjectPlus(psi[N5-1]);
152 
153 
154  // -m fminus[N5-1] P_- psi[0] + f_minus[N5-1] P_+ psi[N5-2]
155  // + f_plus[N5-1] psi[N5-1]
156  //
157  // fminus[N5-1] ( -m P_- psi[0] + P_+ psi[N5-2] )
158  // + f_plus[N5-1] psi[N5-1]
159  //
160  // fminus[N5-1]/2 ( (1 + g5) psi[N5-2] - m(1-g5) psi[0] )
161  // + f_plus[N5-1] psi[N5-1]
162  //
163  // fminus[N5-1]/2 ( psi[N5-2] - m psi[0] + g5( psi[N5-2] +m psi[0] ) )
164  // + f_plus[N5-1] psi[N5-1]
165  //
166  // fminus[N5-1]/2 ( tmp1 + g5 ( tmp2 )) + f_plus[N5-1] psi[N5-1]
167  //
168  // tmp1 = psi[N5-2] - m psi[0]
169  // tmp2 = psi[N5-2] + m psi[0]
170 
171  /* OLD CODE */
172  /*
173  fact = Real(0.5)*f_minus[N5-1];
174  tmp1[rb[cb]] = psi[N5-2] - m_q * psi[0];
175  tmp2[rb[cb]] = psi[N5-2] + m_q * psi[0];
176 
177  chi[N5-1][rb[cb]] = f_plus[N5-1]*psi[N5-1]
178  + fact*( tmp1 + GammaConst<Ns,Ns*Ns-1>() * tmp2);
179  */
180 
181  // Recoded using chiral projectors
182  // -m fminus[N5-1] P_- psi[0] + f_minus[N5-1] P_+ psi[N5-2]
183  // + f_plus[N5-1] psi[N5-1]
184  //
185  // fminus[N5-1] ( -m P_- psi[0] + P_+ psi[N5-2] )
186  // + f_plus[N5-1] psi[N5-1]
187  fact = f_minus[N5-1]*m_q;
188  chi[N5-1][rb[cb]] = f_plus[N5-1]*psi[N5-1];
189  chi[N5-1][rb[cb]] += f_minus[N5-1]*chiralProjectPlus(psi[N5-2]);
190  chi[N5-1][rb[cb]] -= fact*chiralProjectMinus(psi[0]);
191 
192 
193  for(int s=1; s<N5-1; s++) {
194 
195 
196  // fminus[s] P_+ psi[s-1] + fplus[s] psi[s] + fminus[s] P_- psi[s+1]
197  //= fplus[s] psi[s] + fminus[s]( P_+ psi[s-1] + P_-psi[s+1])
198  //= fplus[s] psi[s] + (fminus[s]/2)( psi[s-1] + psi[s+1] +
199  // g_5*{ psi[s-1] - psi[s + 1] } )
200 
201  // OLD CODE
202  /*
203  fact = Real(0.5)*f_minus[s];
204  tmp1[rb[cb]] = psi[s-1] + psi[s+1];
205  tmp2[rb[cb]] = psi[s-1] - psi[s+1];
206 
207  chi[s][rb[cb]] = f_plus[s]*psi[s] +
208  fact*( tmp1 + GammaConst<Ns,Ns*Ns-1>() * tmp2 );
209  */
210 
211  // Recoded using chiral projectors
212  chi[s][rb[cb]] = f_plus[s]*psi[s];
213  chi[s][rb[cb]] += f_minus[s]*chiralProjectMinus(psi[s+1]);
214  chi[s][rb[cb]] += f_minus[s]*chiralProjectPlus(psi[s-1]);
215 
216  }
217 
218 
219  }
220  break ;
221 
222  case MINUS:
223  {
224 
225  // Scarily different. Daggering in the 5th dimension
226  // changes the f_minuses from constant along a row to
227  // varying along a row...
228 
229  Real fact;
230 
231  // OLD CODE
232  /*
233  LatticeFermion tmp1, tmp2;
234 
235  // Fact is constant now.
236  fact = Real(0.5);
237 
238 
239  tmp1[rb[cb]] = f_minus[1]*psi[1] - m_q*f_minus[N5-1]*psi[N5-1];
240  tmp2[rb[cb]] = f_minus[1]*psi[1] + m_q*f_minus[N5-1]*psi[N5-1];
241 
242  chi[0][rb[cb]] = f_plus[0]*psi[0] +
243  fact*( tmp1 + GammaConst<Ns, Ns*Ns-1>()*tmp2 );
244  */
245 
246  // Recoded using chiral projectors
247  fact = m_q * f_minus[N5-1];
248  chi[0][rb[cb]] = f_plus[0]*psi[0];
249  chi[0][rb[cb]] += f_minus[1]*chiralProjectPlus(psi[1]);
250  chi[0][rb[cb]] -= fact*chiralProjectMinus(psi[N5-1]);
251 
252 
253  // OLD CODE
254  /*
255  tmp1[rb[cb]] = f_minus[N5-2]*psi[N5-2] - m_q * f_minus[0]* psi[0];
256  tmp2[rb[cb]] = f_minus[N5-2]*psi[N5-2] + m_q * f_minus[0]* psi[0];
257 
258  chi[N5-1][rb[cb]] = f_plus[N5-1]*psi[N5-1]
259  + fact*( tmp1 - GammaConst<Ns,Ns*Ns-1>()*tmp2 ) ;
260  */
261 
262  // Recoded using Chiral Projectors
263  fact = m_q * f_minus[0];
264  chi[N5-1][rb[cb]] = f_plus[N5-1]*psi[N5-1];
265  chi[N5-1][rb[cb]] += f_minus[N5-2]*chiralProjectMinus(psi[N5-2]);
266  chi[N5-1][rb[cb]] -= fact*chiralProjectPlus(psi[0]);
267 
268 
269  for(int s=1; s<N5-1; s++) {
270 
271  // OLD CODE
272  /*
273  tmp1[rb[cb]] = f_minus[s-1]*psi[s-1] + f_minus[s+1]*psi[s+1];
274  tmp2[rb[cb]] = f_minus[s-1]*psi[s-1] - f_minus[s+1]*psi[s+1];
275 
276  chi[s][rb[cb]] = f_plus[s]*psi[s] +
277  fact*( tmp1 - GammaConst<Ns,Ns*Ns-1>()*tmp2 );
278  */
279 
280  // Recoded using Chiral Projectors
281  chi[s][rb[cb]] = f_plus[s]*psi[s];
282  chi[s][rb[cb]] += f_minus[s-1]*chiralProjectMinus(psi[s-1]);
283  chi[s][rb[cb]] += f_minus[s+1]*chiralProjectPlus(psi[s+1]);
284 
285  }
286 
287  }
288  break ;
289  }
290 
291  END_CODE();
292  }
293 
294 
295  //! Apply the inverse even-even (odd-odd) coupling piece of the domain-wall fermion operator
296  /*!
297  * \ingroup linop
298  *
299  * The operator acts on the entire lattice
300  *
301  * \param psi Pseudofermion field (Read)
302  * \param isign Flag ( PLUS | MINUS ) (Read)
303  * \param cb checkerboard ( 0 | 1 ) (Read)
304  */
305  void
306  EvenOddPrecGenNEFDWLinOpArray::applyDiagInv(multi1d<LatticeFermion>& chi,
307  const multi1d<LatticeFermion>& psi,
308  enum PlusMinus isign,
309  int cb) const
310  {
311  START_CODE();
312 
313  if( chi.size() != N5 ) chi.resize(N5);
314 
315  // I use two temporaries
316  multi1d<LatticeFermion> z(N5); moveToFastMemoryHint(z);
317  multi1d<LatticeFermion> z_prime(N5); moveToFastMemoryHint(z_prime);
318  Real fact;
319 
320  switch ( isign ) {
321 
322  case PLUS:
323  {
324 
325  // First apply the inverse of Lm :
326  // Solve Lm z = chi
327  z[N5-1][rb[cb]] = psi[N5-1];
328  for(int s=0; s < N5-1; s++){
329  z[s][rb[cb]] = psi[s];
330 
331 
332  // OLD CODE
333  // The factor of 1/2 is for the projection expression
334  // fact = Real(0.5)*l[s];
335  // z[N5-1][rb[cb]] -= fact*(psi[s] - GammaConst<Ns,Ns*Ns-1>()*psi[s]) ;
336 
337  // Recoded with projector
338  z[N5-1][rb[cb]] -= l[s]*chiralProjectMinus(psi[s]);
339 
340  }
341 
342  //Now apply the inverse of L. Forward elimination
343  //
344  // L z' = z
345  z_prime[0][rb[cb]] = z[0];
346  for(int s = 0; s < N5-1; s++) {
347 
348  // OLD CODE
349  // The factor of 1/2 is for the projection part
350  // fact = Real(0.5)*a[s];
351  // z_prime[s+1][rb[cb]] = z[s+1] - fact*(z_prime[s] + GammaConst<Ns,Ns*Ns-1>()*z_prime[s]);
352 
353  // Recoded using chiralProjectors
354  z_prime[s+1][rb[cb]] = z[s+1] - a[s]*chiralProjectPlus(z_prime[s]);
355 
356  }
357 
358  // z = D^{-1} z'
359  for(int s=0; s < N5; s++) {
360  fact = Real(1)/d[s];
361  z[s][rb[cb]] = fact*z_prime[s];
362  }
363 
364  // z' = U^{-1} z => U z' = z
365  //The inverse of R. Back substitution...... Getting there!
366  z_prime[N5-1][rb[cb]] = z[N5-1];
367  for(int s=N5-2; s >=0; s-- ) {
368 
369  // OLD CODE
370  // fact = Real(0.5)*b[s];
371 
372  // z_prime[s][rb[cb]] = z[s] - fact*(z_prime[s+1] - GammaConst<Ns,Ns*Ns-1>()*z_prime[s+1]);
373 
374  // Recoded Using ChiralProjectors
375  z_prime[s][rb[cb]] = z[s] - b[s]*chiralProjectMinus(z_prime[s+1]);
376  }
377 
378  //Finally the inverse of Rm
379  chi[N5-1][rb[cb]] = z_prime[N5-1];
380  for(int s=0; s < N5-1; s++) {
381 
382  // OLD CODE
383  // fact = Real(0.5)*r[s];
384 
385  // chi[s][rb[cb]] = z_prime[s] - fact*(chi[N5-1] + GammaConst<Ns,Ns*Ns-1>()*chi[N5-1]);
386 
387  // Recoded using ChiralProjectors
388  chi[s][rb[cb]] = z_prime[s] - r[s]*chiralProjectPlus(chi[N5-1]);
389 
390  }
391  }
392  break ;
393 
394  case MINUS:
395  {
396 
397  // First apply the inverse of Rm :
398  // Solve (Rm)^{T} z = psi
399 
400  // This is the same as the Lm case above but with l P_ => r P+
401 
402  z[N5-1][rb[cb]] = psi[N5-1];
403  for(int s=0; s < N5-1; s++){
404  z[s][rb[cb]] = psi[s];
405 
406  // OLD CODE
407  // The factor of 1/2 is for the projection expression
408  // fact = Real(0.5)*r[s];
409  // z[N5-1][rb[cb]] -= fact*(psi[s] + GammaConst<Ns,Ns*Ns-1>()*psi[s]) ;
410  // Recoded Using Chiral Projectors
411  z[N5-1][rb[cb]] -= r[s]*chiralProjectPlus(psi[s]);
412  }
413 
414  //Now apply the inverse of U^{T}. Forward elimination
415  //
416  // U^T z' = z . This is the same as L z' = z above but with
417  // a P+ => b P-
418 
419  z_prime[0][rb[cb]] = z[0];
420  for(int s = 0; s < N5-1; s++) {
421 
422  // OLD CODE
423  // The factor of 1/2 is for the projection part
424  // fact = Real(0.5)*b[s];
425  // z_prime[s+1][rb[cb]] = z[s+1] - fact*(z_prime[s] - GammaConst<Ns,Ns*Ns-1>()*z_prime[s]);
426 
427  // Recoded using Chiral Projectors
428  z_prime[s+1][rb[cb]] = z[s+1] - b[s]*chiralProjectMinus(z_prime[s]);
429  }
430 
431  // z = D^{-1} z'
432  for(int s=0; s < N5; s++) {
433  fact = Real(1)/d[s];
434  z[s][rb[cb]] = fact*z_prime[s];
435  }
436 
437  // z' = (L^T)^{-1} z ie L^T z' = z
438  //
439  // This is the same as the U z' =z case above
440  // but with b P_ => a P+
441 
442  z_prime[N5-1][rb[cb]] = z[N5-1];
443  for(int s=N5-2; s >=0; s-- ) {
444 
445  // OLD CODE
446  // fact = Real(0.5)*a[s];
447  //z_prime[s][rb[cb]] = z[s] - fact*(z_prime[s+1] + GammaConst<Ns,Ns*Ns-1>()*z_prime[s+1]);
448 
449  // Recoded using Chiral Projectors
450  z_prime[s][rb[cb]] = z[s] - a[s]*chiralProjectPlus(z_prime[s+1]);
451  }
452 
453  //Finally the inverse of Lm^T
454  //Same as the inverse of R above but with r P+ => l P-
455  chi[N5-1][rb[cb]] = z_prime[N5-1];
456  for(int s=0; s < N5-1; s++) {
457 
458  // OLD CODE
459  // fact = Real(0.5)*l[s];
460  // chi[s][rb[cb]] = z_prime[s] - fact*(chi[N5-1] - GammaConst<Ns,Ns*Ns-1>()*chi[N5-1]);
461 
462 
463  // Recoded using Chiral Projectors
464  chi[s][rb[cb]] = z_prime[s] - l[s]*chiralProjectMinus(chi[N5-1]);
465  }
466 
467  }
468  break ;
469  }
470 
471  END_CODE();
472  }
473 
474  //! Apply the even-odd (odd-even) coupling piece of the NEF operator
475  /*!
476  * \ingroup linop
477  *
478  * The operator acts on the entire lattice
479  *
480  * \param psi Pseudofermion field (Read)
481  * \param isign Flag ( PLUS | MINUS ) (Read)
482  * \param cb checkerboard ( 0 | 1 ) (Read)
483  */
484  void
485  EvenOddPrecGenNEFDWLinOpArray::applyOffDiag(multi1d<LatticeFermion>& chi,
486  const multi1d<LatticeFermion>& psi,
487  enum PlusMinus isign,
488  int cb) const
489  {
490  START_CODE();
491 
492  if( chi.size() != N5 ) chi.resize(N5);
493 
494  switch ( isign )
495  {
496  case PLUS:
497  {
498  multi1d<LatticeFermion> tmp(N5); moveToFastMemoryHint(tmp);
499  Real fact1;
500  Real fact2;
501  Real fact2mf;
502 
503  int otherCB = (cb + 1)%2 ;
504 
505  // OLD CODE
506  /*
507  LatticeFermion tmp1;
508  LatticeFermion tmp2;
509 
510 
511  fact1 = -Real(0.5)*b5[0];
512  fact2 = -Real(0.25)*c5[0];
513  tmp1[rb[otherCB]] = psi[1] - m_q*psi[N5-1];
514  tmp2[rb[otherCB]] = psi[1] + m_q*psi[N5-1];
515 
516  tmp[0][rb[otherCB]] = fact1*psi[0]
517  + fact2*(tmp1 - GammaConst<Ns, Ns*Ns-1>()*tmp2);
518  */
519 
520  // Recoded using Chiral Projectors
521  fact1 = -Real(0.5)*b5[0];
522  fact2 = -Real(0.5)*c5[0];
523  fact2mf = m_q*fact2;
524 
525  tmp[0][rb[otherCB]] = fact1*psi[0];
526  tmp[0][rb[otherCB]] += fact2*chiralProjectMinus(psi[1]);
527  tmp[0][rb[otherCB]] -= fact2mf*chiralProjectPlus(psi[N5-1]);
528 
529 
530  // OLD CODE
531  /*
532  fact1 = -Real(0.5)*b5[N5-1];
533  fact2 = -Real(0.25)*c5[N5-1];
534 
535  tmp1[rb[otherCB]]=psi[N5-2]-m_q*psi[0];
536  tmp2[rb[otherCB]]=psi[N5-2]+m_q*psi[0];
537 
538  tmp[N5-1][rb[otherCB]] = fact1*psi[N5-1]
539  + fact2*(tmp1 + GammaConst<Ns,Ns*Ns-1>()*tmp2 );
540  */
541 
542  // Recoded using Chiral Projectors
543  fact1 = -Real(0.5)*b5[N5-1];
544  fact2 = -Real(0.5)*c5[N5-1];
545  fact2mf = m_q*fact2;
546 
547  tmp[N5-1][rb[otherCB]] = fact1*psi[N5-1];
548  tmp[N5-1][rb[otherCB]] += fact2*chiralProjectPlus(psi[N5-2]);
549  tmp[N5-1][rb[otherCB]] -= fact2mf*chiralProjectMinus(psi[0]);
550 
551  for(int s=1; s < N5-1; s++) {
552 
553  // OLD CODE
554  /*
555  fact1 = -Real(0.5)*b5[s];
556  fact2 = -Real(0.25)*c5[s];
557 
558  tmp1[rb[otherCB]]=psi[s-1] + psi[s+1];
559  tmp2[rb[otherCB]]=psi[s-1] - psi[s+1];
560 
561  tmp[s][rb[otherCB]]= fact1*psi[s]
562  + fact2*(tmp1 + GammaConst<Ns,Ns*Ns-1>()*tmp2) ;
563  */
564 
565  // Recoded using Chiral Projectors
566  fact1 = -Real(0.5)*b5[s];
567  fact2 = -Real(0.5)*c5[s];
568  tmp[s][rb[otherCB]] = fact1*psi[s];
569  tmp[s][rb[otherCB]] += fact2*chiralProjectPlus(psi[s-1]);
570  tmp[s][rb[otherCB]] += fact2*chiralProjectMinus(psi[s+1]);
571  }
572 
573  // Replace with std::vector dslash later -- done
574  D.apply(chi, tmp, isign, cb);
575 
576 
577  }
578  break ;
579 
580  case MINUS:
581  {
582  multi1d<LatticeFermion> tmp_d(N5) ; moveToFastMemoryHint(tmp_d);
583 
584  // Replace with a std::vector dslash later -- done
585  D.apply(tmp_d, psi, isign, cb);
586 
587 
588  // OLD CODE
589  /*
590  LatticeFermion tmp1;
591  LatticeFermion tmp2;
592 
593 
594  Real one_quarter = Real(0.25);
595 
596  Real factb= -Real(0.5)*b5[0];
597  Real fact = m_q*c5[N5-1];
598  tmp1[rb[cb]] = c5[1]*tmp_d[1] - fact*tmp_d[N5-1];
599  tmp2[rb[cb]] = c5[1]*tmp_d[1] + fact*tmp_d[N5-1];
600 
601  chi[0][rb[cb]] = factb*tmp_d[0]
602  - one_quarter*( tmp1 + GammaConst<Ns,Ns*Ns-1>()*tmp2 );
603  */
604 
605  // Recoded with Chiral Projector
606  Real factb = -Real(0.5)*b5[0];
607  Real factc1 = -Real(0.5)*c5[1];
608  Real factc2 = -Real(0.5)*m_q*c5[N5-1];
609 
610  chi[0][rb[cb]] = factb*tmp_d[0];
611  chi[0][rb[cb]] += factc1*chiralProjectPlus(tmp_d[1]);
612  chi[0][rb[cb]] -= factc2*chiralProjectMinus(tmp_d[N5-1]);
613 
614 
615  // OLD CODE
616  /*
617  factb = -Real(0.5)*b5[N5-1];
618  fact = m_q * c5[0];
619 
620  tmp1[rb[cb]] = c5[N5-2]*tmp_d[N5-2] - fact * tmp_d[0];
621  tmp2[rb[cb]] = c5[N5-2]*tmp_d[N5-2] + fact * tmp_d[0];
622 
623  chi[N5-1][rb[cb]] = factb * tmp_d[N5-1]
624  - one_quarter * ( tmp1 - GammaConst<Ns, Ns*Ns-1>()*tmp2 );
625  */
626 
627  // Recoded with Chiral Projector
628  factb = -Real(0.5)*b5[N5-1];
629  factc1 = -Real(0.5)*c5[N5-2];
630  factc2 = -Real(0.5)*m_q*c5[0];
631 
632  chi[N5-1][rb[cb]] = factb*tmp_d[N5-1];
633  chi[N5-1][rb[cb]] += factc1*chiralProjectMinus(tmp_d[N5-2]);
634  chi[N5-1][rb[cb]] -= factc2*chiralProjectPlus(tmp_d[0]);
635 
636 
637  for(int s=1; s<N5-1; s++){
638  factb = -Real(0.5) * b5[s];
639  factc1 = -Real(0.5) * c5[s-1];
640  factc2 = -Real(0.5) * c5[s+1];
641 
642  // OLD CODE
643  // tmp1[rb[cb]] = c5[s-1]*tmp_d[s-1] + c5[s+1]*tmp_d[s+1];
644  // tmp2[rb[cb]] = c5[s-1]*tmp_d[s-1] - c5[s+1]*tmp_d[s+1];
645 
646  // chi[s][rb[cb]] = factb*tmp_d[s]
647  // - one_quarter*( tmp1 - GammaConst<Ns, Ns*Ns-1>()*tmp2 );
648 
649  // Recoded with Chiral Projector
650  chi[s][rb[cb]] = factb*tmp_d[s];
651  chi[s][rb[cb]] += factc1*chiralProjectMinus(tmp_d[s-1]);
652  chi[s][rb[cb]] += factc2*chiralProjectPlus(tmp_d[s+1]);
653 
654  }
655  }
656  break ;
657  }
658 
659 
660  END_CODE();
661  }
662 
663 
664  //! Apply the Dminus operator on a lattice fermion. See my notes ;-)
665  void
666  EvenOddPrecGenNEFDWLinOpArray::Dminus(LatticeFermion& chi,
667  const LatticeFermion& psi,
668  enum PlusMinus isign,
669  int s5) const
670  {
671  LatticeFermion tt ; moveToFastMemoryHint(tt);
672  D.apply(tt,psi,isign,0);
673  D.apply(tt,psi,isign,1);
674  Real fact = Real(0.5)*c5[s5];
675  chi = f_minus[s5]*psi - fact*tt ;
676  }
677 
678 
679 
680  // Derivative of even-odd linop component
681  /*
682  * This is a copy of the above applyOffDiag with the D.apply(...) replaced
683  * by D.deriv(ds_tmp,...) like calls.
684  */
685  void
686  EvenOddPrecGenNEFDWLinOpArray::applyDerivOffDiag(multi1d<LatticeColorMatrix>& ds_u,
687  const multi1d<LatticeFermion>& chi,
688  const multi1d<LatticeFermion>& psi,
689  enum PlusMinus isign,
690  int cb) const
691  {
692  START_CODE();
693 
694  ds_u.resize(Nd);
695  multi1d<LatticeColorMatrix> ds_tmp(Nd);
696 
697  ds_u = zero;
698 
699  switch ( isign ) {
700  case PLUS:
701  {
702  LatticeFermion tmp; moveToFastMemoryHint(tmp);
703  Real fact1;
704  Real fact2;
705  Real fact2mf;
706  int otherCB = (cb + 1)%2 ;
707 
708  // OLD CODE
709  /*
710  LatticeFermion tmp1;
711  LatticeFermion tmp2;
712 
713 
714  fact1 = -Real(0.5)*b5[0];
715  fact2 = -Real(0.25)*c5[0];
716  tmp1[rb[otherCB]] = psi[1] - m_q*psi[N5-1];
717  tmp2[rb[otherCB]] = psi[1] + m_q*psi[N5-1];
718 
719  tmp[rb[otherCB]] = fact1*psi[0]
720  + fact2*(tmp1 - GammaConst<Ns, Ns*Ns-1>()*tmp2);
721  */
722 
723  // Recoded with chiralProject
724  fact1 = -Real(0.5)*b5[0];
725  fact2 = -Real(0.5)*c5[0];
726  fact2mf = m_q*fact2;
727 
728  tmp[rb[otherCB]] = fact1*psi[0];
729  tmp[rb[otherCB]] += fact2*chiralProjectMinus(psi[1]);
730  tmp[rb[otherCB]] -= fact2mf*chiralProjectPlus(psi[N5-1]);
731 
732  D.deriv(ds_u, chi[0], tmp, isign, cb);
733 
734  /*
735  fact1 = -Real(0.5)*b5[N5-1];
736  fact2 = -Real(0.25)*c5[N5-1];
737 
738  tmp1[rb[otherCB]]=psi[N5-2]-m_q*psi[0];
739  tmp2[rb[otherCB]]=psi[N5-2]+m_q*psi[0];
740 
741  tmp[rb[otherCB]] = fact1*psi[N5-1]
742  + fact2*(tmp1 + GammaConst<Ns,Ns*Ns-1>()*tmp2 );
743  */
744 
745  // Recoded using Chiral Projectors
746  fact1 = -Real(0.5)*b5[N5-1];
747  fact2 = -Real(0.5)*c5[N5-1];
748  fact2mf = m_q*fact2;
749 
750  tmp[rb[otherCB]] = fact1*psi[N5-1];
751  tmp[rb[otherCB]] += fact2*chiralProjectPlus(psi[N5-2]);
752  tmp[rb[otherCB]] -= fact2mf*chiralProjectMinus(psi[0]);
753 
754  D.deriv(ds_tmp, chi[N5-1], tmp, isign, cb);
755  ds_u += ds_tmp;
756 
757 
758  for(int s=1; s < N5-1; s++) {
759 
760  // OLD CODE
761  /*
762  fact1 = -Real(0.5)*b5[s];
763  fact2 = -Real(0.25)*c5[s];
764 
765  tmp1[rb[otherCB]]=psi[s-1] + psi[s+1];
766  tmp2[rb[otherCB]]=psi[s-1] - psi[s+1];
767 
768  tmp[rb[otherCB]]= fact1*psi[s]
769  + fact2*(tmp1 + GammaConst<Ns,Ns*Ns-1>()*tmp2) ;
770  */
771 
772  // Recoded using Chiral Projectors
773  fact1 = -Real(0.5)*b5[s];
774  fact2 = -Real(0.5)*c5[s];
775  tmp[rb[otherCB]] = fact1*psi[s];
776  tmp[rb[otherCB]] += fact2*chiralProjectPlus(psi[s-1]);
777  tmp[rb[otherCB]] += fact2*chiralProjectMinus(psi[s+1]);
778  D.deriv(ds_tmp, chi[s], tmp, isign, cb);
779  ds_u += ds_tmp;
780  }
781  }
782  break ;
783 
784  case MINUS:
785  {
786 
787  // New Minus case by Balint
788  // There are 3 cases which are summed over.
789  // Best way to compute them is to go back to basics.
790 
791  // s=0:
792  //
793  // F = -(b_5[0]/2) chi(0,cb)^dag Ddot(dagger, cb) psi(0, 1-cb)
794  // -(c_5[1]/2) chi(0,cb)^dag P_{+} Ddot(dagger, cb) psi(1, 1-cb)
795  // +( m c_5[N-1]/2) chi(0,cb)^dag P_{-} Ddot(dagger, cb) psi(N-1, 1-cb)
796 
797  // General 0 < s < N-1:
798  //
799  // F = -(b_5[s]/2) chi(s,cb)^dag Ddot(dagger, cb) psi(s, 1-cb)
800  // -(c_5[s-1]/2) chi(s,cb)^dag P_{-} Ddot(dagger, cb) psi(s-1, 1-cb)
801  // -(c_5[s+1]/2) chi(s,cb)^dag P_{+} Ddot(dagger, cb) psi(s+1, 1-cb)
802 
803  // and for s = N-1
804  //
805  // F = -(b_5[N-1]/2) chi(N-1,cb)^dag Ddot(dagger,cb) psi(N-1, 1-cb)
806  // -(c_5[N-2]/2) chi(N-1,cb)^dag P_{-} Ddot(dagger,cb) psi(N-2, 1-cb)
807  // +( m c_5[0] /2) chi(N-1,cb)^dag P_{+} Ddot(dagger,cb) psi(0, 1-cb)
808  //
809  //
810  // If we compute X, P_{+}X, and P_{-}X up front then we need to
811  // evaluate 3 force terms for each i. If we multiply out the projections // and do gamma_5's etc it brings the number of terms up to 6.
812  //
813  // We also fold the coefficients into the right hand (chi) vectors
814  // as needed
815 
816  // Storage for P_{+} X
817  multi1d<LatticeFermion> chi_plus(N5); moveToFastMemoryHint(chi_plus);
818  multi1d<LatticeFermion> chi_minus(N5); moveToFastMemoryHint(chi_minus);
819 
820  for(int s=0; s < N5; s++) {
821  // OLD CODE
822  /*
823  chi_plus[s][rb[cb]] = chi[s] + GammaConst<Ns,Ns*Ns-1>()*chi[s];
824  chi_plus[s][rb[cb]] *= Real(0.5);
825 
826  chi_minus[s][rb[cb]] = chi[s] - GammaConst<Ns,Ns*Ns-1>()*chi[s];
827  chi_minus[s][rb[cb]]*= Real(0.5);
828  */
829 
830  // Recoded using chiralProject
831  chi_plus[s][rb[cb]] = chiralProjectPlus(chi[s]);
832  chi_minus[s][rb[cb]] = chiralProjectMinus(chi[s]);
833 
834  // Zero out unwanted checkerboard... May remove this later
835  chi_plus[s][rb[1-cb]] = zero;
836  chi_minus[s][rb[1-cb]] = zero;
837  }
838 
839 
840  ds_u = zero;
841  LatticeFermion chi_tmp = zero; moveToFastMemoryHint(chi_tmp);
842 
843 
844  Real ftmp;
845 
846  // First term:
847  // -(b5[0]/2) chi[0]^dag Ddot^dag psi[0]
848  ftmp=Real(0.5)*b5[0];
849  chi_tmp[rb[cb]] = ftmp*chi[0];
850  D.deriv(ds_tmp, chi_tmp, psi[0], isign, cb);
851  ds_u -= ds_tmp;
852 
853 
854  // -(c5[1]/2) chi[0]^dag P+ Ddot^dag psi[1]
855  ftmp=Real(0.5)*c5[1];
856  chi_tmp[rb[cb]] = ftmp*chi_plus[0];
857  D.deriv(ds_tmp, chi_tmp, psi[1], isign, cb);
858  ds_u -= ds_tmp;
859 
860  // +(m c5[N-1]) chi[0]^dag P_- Ddot^dag psi[N-1]
861  ftmp=Real(0.5)*c5[N5-1]*m_q;
862  chi_tmp[rb[cb]] = ftmp*chi_minus[0];
863  D.deriv(ds_tmp, chi_tmp, psi[N5-1], isign, cb);
864  ds_u += ds_tmp;
865 
866 
867  // Middle bits
868  for(int s=1; s < N5-1; s++) {
869 
870  // -(b5[s]/2) chi[s] Ddot^dag psi[s]
871  ftmp = Real(0.5)*b5[s];
872  chi_tmp[rb[cb]] = ftmp*chi[s];
873  D.deriv(ds_tmp, chi_tmp, psi[s], isign, cb);
874  ds_u -= ds_tmp;
875 
876  // -(c5[s-1]/2) chi[s]^dag P- Ddot^dag psi[s-1]
877  ftmp = Real(0.5)*c5[s-1];
878  chi_tmp[rb[cb]] = ftmp*chi_minus[s];
879  D.deriv(ds_tmp, chi_tmp, psi[s-1], isign,cb);
880  ds_u -= ds_tmp;
881 
882  // -(c5[s+1]/2) chi[s]^dag P+ Ddot^dag psi[s+1]
883  ftmp = Real(0.5)*c5[s+1];
884  chi_tmp[rb[cb]] = ftmp*chi_plus[s];
885  D.deriv(ds_tmp, chi_tmp, psi[s+1], isign, cb);
886  ds_u -= ds_tmp;
887 
888  }
889 
890 
891  // Lat bit
892  // - (b5[N5-1]/2) chi[N5-1] Ddot^dagger psi[N5-1]
893  ftmp = Real(0.5)*b5[N5-1];
894  chi_tmp[rb[cb]] = ftmp * chi[N5-1];
895  D.deriv(ds_tmp, chi_tmp, psi[N5-1], isign, cb);
896  ds_u -= ds_tmp;
897 
898  // -(c5[N5-2]/2) chi[N5-1] P_ Ddot^dagger psi[N5-2]
899  ftmp = Real(0.5)*c5[N5-2];
900  chi_tmp[rb[cb]] = ftmp*chi_minus[N5-1];
901  D.deriv(ds_tmp, chi_tmp, psi[N5-2], isign, cb);
902  ds_u -= ds_tmp;
903 
904  // +(m c5[0]/2) chi[N5-1] P+ Ddot^dagger psi[0]
905  ftmp = Real(0.5)*c5[0]*m_q;
906  chi_tmp[rb[cb]] = ftmp * chi_plus[N5-1];
907  D.deriv(ds_tmp, chi_tmp, psi[0], isign, cb);
908  ds_u += ds_tmp;
909 
910  }
911  break ;
912  }
913 
914 
915  END_CODE();
916  }
917 
918 } // End Namespace Chroma
919 
chromabase.h
Primary include file for CHROMA library code.
END_CODE
#define END_CODE()
Definition: chromabase.h:65
START_CODE
#define START_CODE()
Definition: chromabase.h:64
Chroma::FermState
Support class for fermion actions and linear operators.
Definition: state.h:94
Chroma::Handle
Class for counted reference semantics.
Definition: handle.h:33
N5
int N5
Definition: eoprec_dwf_fermact_array_w.cc:78
eoprec_nef_general_linop_array_w.h
4D Even Odd preconditioned NEF domain-wall fermion linear operator generalised to take array of b_5 a...
s
unsigned s
Definition: ldumul_w.cc:37
i
unsigned i
Definition: ldumul_w.cc:34
tmp
Double tmp
Definition: meslate.cc:60
z
int z
Definition: meslate.cc:36
Nd
Nd
Definition: meslate.cc:74
tmp2
Double tmp2
Definition: mesq.cc:30
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::MINUS
@ MINUS
Definition: chromabase.h:45
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::a
Complex a
Definition: invbicg.cc:95
Chroma::d
DComplex d
Definition: invbicg.cc:99
Chroma::b
Complex b
Definition: invbicg.cc:96
Chroma::cb
int cb
Definition: invbicg.cc:120
Chroma::zero
Double zero
Definition: invbicg.cc:106
r
multi1d< LatticeFermion > r(Ncb)
chi
chi
Definition: pade_trln_w.cc:24
psi
psi
Definition: pade_trln_w.cc:191
l
int l
Definition: pade_trln_w.cc:111
Generated by doxygen 1.9.1