CHROMA
eoprec_parwilson_linop_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Even-odd preconditioned Wilson fermion linear operator with parity breaking term
3  */
4 
6 
7 namespace Chroma
8 {
9  //! Creation routine
10  /*!
11  * \param u_ gauge field (Read)
12  * \param Mass_ fermion mass (Read)
13  * \param H_ parity breaking term (Read)
14  */
16  const Real& Mass_, const Real& H_)
17  {
18  Mass = Mass_;
19  H = H_;
20 // u = u_;
21  D.create(fs);
22 
23  fact = Nd + Mass;
24  Real tmp = 1.0 / (fact*fact + H*H);
25  invfact1 = fact * tmp;
26  invfact2 = H * tmp;
27  }
28 
29 
30  //! Apply the the even-even block onto a source std::vector
31  void
32  EvenOddPrecParWilsonLinOp::evenEvenLinOp(LatticeFermion& chi, const LatticeFermion& psi,
33  enum PlusMinus isign) const
34  {
35  switch (isign)
36  {
37  case PLUS:
38  chi[rb[0]] = fact*psi + GammaConst<Ns,Ns*Ns-1>()*(H*timesI(psi));
39  break;
40 
41  case MINUS:
42  chi[rb[0]] = fact*psi - GammaConst<Ns,Ns*Ns-1>()*(H*timesI(psi));
43  break;
44  }
45  }
46 
47 
48  //! Return flops performed by the operator()
49  unsigned long EvenOddPrecParWilsonLinOp::nFlops() const
50  {
51  unsigned long cbsite_flops = 2*D.nFlops()+16*Nc*Ns;
52  return cbsite_flops*(Layout::sitesOnNode()/2);
53  }
54 
55 
56  //! Apply the inverse of the even-even block onto a source std::vector
57  void
58  EvenOddPrecParWilsonLinOp::evenEvenInvLinOp(LatticeFermion& chi, const LatticeFermion& psi,
59  enum PlusMinus isign) const
60  {
61  // tmp = D' (1-i isign H gamma_5) D' Psi
62  // O O,E E,O O
63  switch (isign)
64  {
65  case PLUS:
66  chi[rb[0]] = invfact1*psi - GammaConst<Ns,Ns*Ns-1>()*(invfact2*timesI(psi));
67  break;
68 
69  case MINUS:
70  chi[rb[0]] = invfact1*psi + GammaConst<Ns,Ns*Ns-1>()*(invfact2*timesI(psi));
71  break;
72  }
73  }
74 
75 
76  //! Apply the the odd-odd block onto a source std::vector
77  void
78  EvenOddPrecParWilsonLinOp::oddOddLinOp(LatticeFermion& chi, const LatticeFermion& psi,
79  enum PlusMinus isign) const
80  {
81  switch (isign)
82  {
83  case PLUS:
84  chi[rb[1]] = fact*psi + GammaConst<Ns,Ns*Ns-1>()*(H*timesI(psi));
85  break;
86 
87  case MINUS:
88  chi[rb[1]] = fact*psi - GammaConst<Ns,Ns*Ns-1>()*(H*timesI(psi));
89  break;
90  }
91  }
92 
93  //! Apply even-odd linop component
94  /*!
95  * The operator acts on the entire even sublattice
96  *
97  * \param chi Pseudofermion field (Write)
98  * \param psi Pseudofermion field (Read)
99  * \param isign Flag ( PLUS | MINUS ) (Read)
100  */
101  void
103  const LatticeFermion& psi,
104  enum PlusMinus isign) const
105  {
106  START_CODE();
107 
108  Real mhalf = -0.5;
109 
110  D.apply(chi, psi, isign, 0);
111  chi[rb[0]] *= mhalf;
112 
113  END_CODE();
114  }
115 
116 
117  //! Apply odd-even linop component
118  /*!
119  * The operator acts on the entire odd sublattice
120  *
121  * \param chi Pseudofermion field (Write)
122  * \param psi Pseudofermion field (Read)
123  * \param isign Flag ( PLUS | MINUS ) (Read)
124  */
125  void
127  const LatticeFermion& psi,
128  enum PlusMinus isign) const
129  {
130  START_CODE();
131 
132  Real mhalf = -0.5;
133 
134  D.apply(chi, psi, isign, 1);
135  chi[rb[1]] *= mhalf;
136 
137  END_CODE();
138  }
139 
140 
141  //! Override inherited one with a few more funkies
142  void
144  const LatticeFermion& psi,
145  enum PlusMinus isign) const
146  {
147  LatticeFermion tmp1, tmp2, tmp3; // if an array is used here,
148 
149  moveToFastMemoryHint(tmp1);
150  moveToFastMemoryHint(tmp2);
151  moveToFastMemoryHint(tmp3);
152 
153  Real mquarter = -0.25;
154 
155  // tmp = D' (1-i isign H gamma_5) D' Psi
156  // O O,E E,O O
157  switch (isign)
158  {
159  case PLUS:
160  // tmp1[0] = D_eo psi[1]
161  D.apply(tmp1, psi, isign, 0);
162 
163  tmp2[rb[0]] = invfact1*tmp1 - invfact2*(GammaConst<Ns,Ns*Ns-1>()*timesI(tmp1));
164 
165  // tmp2[1] = D_oe tmp2[0]
166  D.apply(tmp3, tmp2, isign, 1);
167 
168  chi[rb[1]] = fact*psi + mquarter*tmp3;
169  chi[rb[1]] += H*(GammaConst<Ns,Ns*Ns-1>()*timesI(psi));
170  break;
171 
172  case MINUS:
173  // tmp1[0] = D_eo psi[1]
174  D.apply(tmp1, psi, isign, 0);
175 
176  tmp2[rb[0]] = invfact1*tmp1 + invfact2*(GammaConst<Ns,Ns*Ns-1>()*timesI(tmp1));
177 
178  // tmp2[1] = D_oe tmp2[0]
179  D.apply(tmp3, tmp2, isign, 1);
180 
181  chi[rb[1]] = fact*psi + mquarter*tmp3;
182  chi[rb[1]] -= H*(GammaConst<Ns,Ns*Ns-1>()*timesI(psi));
183  break;
184  }
185 
186  getFermBC().modifyF(chi, rb[1]);
187  }
188 
189 
190  //! Derivative of even-odd linop component
191  void
192  EvenOddPrecParWilsonLinOp::derivEvenOddLinOp(multi1d<LatticeColorMatrix>& ds_u,
193  const LatticeFermion& chi,
194  const LatticeFermion& psi,
195  enum PlusMinus isign) const
196  {
197  START_CODE();
198 
199  ds_u.resize(Nd);
200 
201  D.deriv(ds_u, chi, psi, isign, 0);
202  for(int mu=0; mu < Nd; mu++) {
203  ds_u[mu] *= Real(-0.5);
204  }
205 
206  END_CODE();
207  }
208 
209 
210  //! Derivative of odd-even linop component
211  void
212  EvenOddPrecParWilsonLinOp::derivOddEvenLinOp(multi1d<LatticeColorMatrix>& ds_u,
213  const LatticeFermion& chi,
214  const LatticeFermion& psi,
215  enum PlusMinus isign) const
216  {
217  START_CODE();
218 
219  ds_u.resize(Nd);
220 
221  D.deriv(ds_u, chi, psi, isign, 1);
222  for(int mu=0; mu < Nd; mu++) {
223  ds_u[mu] *= Real(-0.5);
224  }
225  END_CODE();
226  }
227 
228 
229 
230 #if 0
231 // Code is here only as a reference. It should be deleted at some point.
232 // This is converted szin code.
233 // The above derivs will work, calling evenEvenInvLinOp to get that
234 // contribution.
235 // However, for more performance, one should copy the regular Wilson
236 // deriv() routine and modify it.
237 
238 #error "ANCIENT SZIN CODE"
239 #error "Not quite correct implementation"
240 
241 
242  //! Computes the derivative of the fermionic action respect to the link field
243  /*!
244  * | dS dS_f
245  * ds_u -- | ---- + ----- ( Write )
246  * | dU dU
247  *
248  * psi -- [1./(M_dag*M)]*chi_ ( read )
249  *
250  * \param ds_u result ( Write )
251  * \param state gauge field ( Read )
252  * \param psi solution to linear system ( Read )
253  */
254 
255  void
256  EvenOddPrecParWilsonFermAct::deriv(multi1d<LatticeColorMatrix>& ds_u,
257  const LatticeFermion& chi,
258  const LatticeFermion& psi,
259  enum PlusMinus isign) const
260  {
261  START_CODE();
262 
263  QDPIO::cerr << "EvenOddPrecParWilsonFermAct::dsdu not implemented" << std::endl;
264  QDP_abort(1);
265 
266  const multi1d<LatticeColorMatrix>& u = state->getLinks();
267 
268  LatticeColorMatrix utmp_1; moveToFastMemoryHint(utmp_1);
269  LatticeFermion phi; moveToFastMemoryHint(phi);
270  LatticeFermion rho; moveToFastMemoryHint(rho);
271  LatticeFermion sigma; moveToFastMemoryHint(sigma);
272  LatticeFermion ftmp_1; moveToFastMemoryHint(ftmp_1);
273  LatticeFermion ftmp_2; moveToFastMemoryHint(ftmp_2);
274  Double ddummy;
275  Real dummy;
276  int nu;
277  int cb;
278 
279  /* Do the Wilson fermion dS_f/dU with parity breaking term */
280 
281 // CONSTRUCT_LINEAR_OPERATOR(A, lwlhmpsi, u, KappaMD);
282 
283  /* phi = M(u)*psi */
284  A (A, psi, phi, 1, PLUS);
285 
286  /* rho = (1-i H gamma_5) * Dslash(0<-1) * psi */
287  dslash (u, psi, ftmp_1, PLUS, 1);
288  PARBREAK(ftmp_1, H_parity, rho, MINUS);
289 
290  /* phi = (KappaMD^2)*phi/(1+h^2) = -(KappaMD^2)*M*psi/(1+h^2) */
291  dummy = -(KappaMD*KappaMD) / (WORD_VALUE(WORD_H_parity, ONE) + H_parity*H_parity);
292  phi = phi * dummy;
293 
294  /* sigma = (1+i H gamma_5) * Dslash_dag(0<-1) * phi */
295  dslash (u, phi, ftmp_1, MINUS, 1);
296  PARBREAK(ftmp_1, H_parity, sigma, PLUS);
297 
298 
299  for(int mu = 0; mu < Nd; ++mu)
300  {
301  cb = 0;
302 
303  /* ftmp_2 = (gamma(mu))*psi */
304  SPIN_PRODUCT(psi,(INTEGER_LSHIFT_FUNCTION(1,mu)),ftmp_2);
305  /* ftmp_2(x) = -(psi(x) - ftmp_2(x)) = -(1 - gamma(mu))*psi( x ) */
306  ftmp_2 -= psi;
307  utmp_1 = -(shift(ftmp_2, cb, FORWARD, mu) * adj(sigma));
308 
309  /* ftmp_2 = (gamma(mu))*phi */
310  SPIN_PRODUCT(phi,(INTEGER_LSHIFT_FUNCTION(1,mu)),ftmp_2);
311  /* ftmp_2 = phi + ftmp_2 = (1 + gamma(mu))*phi( x) */
312  ftmp_2 += phi;
313  utmp_1 += shift(ftmp_2, cb, FORWARD, mu) * adj(rho);
314 
315  /* THIS NEEDS TO BE CHANGED (removed) */
316  ds_u[mu][cb] += u[mu][cb] * utmp_1;
317 
318  cb = 1;
319 
320  /* ftmp_2 = (gamma(mu))*ftmp_1 */
321  SPIN_PRODUCT(rho,(INTEGER_LSHIFT_FUNCTION(1,mu)),ftmp_2);
322  /* ftmp_2 = -(rho - ftmp_2) = -(1 - gamma(mu))*rho( x ) */
323  ftmp_2 -= rho;
324  utmp_1 = -(shift(ftmp_2, cb, FORWARD, mu) * adj(phi));
325 
326  /* ftmp_2 = (gamma(mu))*sigma */
327  SPIN_PRODUCT(sigma,(INTEGER_LSHIFT_FUNCTION(1,mu)),ftmp_2);
328  /* ftmp_1 = ftmp_1 + ftmp_2 = (1 + gamma(mu))*sigma( x + mu) */
329  ftmp_2 += sigma;
330  utmp_1 += shift(ftmp_2, cb, FORWARD, mu) * adj(psi);
331 
332  /* THIS NEEDS TO BE CHANGED (removed) */
333  ds_u[mu][cb] += u[mu][cb] * utmp_1;
334 
335  }
336 
337  getFermBC().zero(ds_u);
338 
339  END_CODE();
340  }
341 #endif
342 
343 } // End Namespace Chroma
344 
void evenEvenInvLinOp(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the inverse of the even-even block onto a source std::vector.
void oddEvenLinOp(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the odd-even block onto a source std::vector.
void evenOddLinOp(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the even-odd block onto a source std::vector.
void create(Handle< FermState< T, P, Q > > fs, const Real &Mass_, const Real &H_)
Creation routine.
unsigned long nFlops() const
Return flops performed by the operator()
void derivEvenOddLinOp(multi1d< LatticeColorMatrix > &ds_u, const LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the even-odd block onto a source std::vector.
void operator()(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Override inherited one with a few more funkies.
void evenEvenLinOp(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the even-even block onto a source std::vector.
const FermBC< T, P, Q > & getFermBC() const
Return the fermion BC object for this linear operator.
void derivOddEvenLinOp(multi1d< LatticeColorMatrix > &ds_u, const LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the odd-even block onto a source std::vector.
void oddOddLinOp(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the the odd-odd block onto a source std::vector.
Support class for fermion actions and linear operators.
Definition: state.h:94
virtual const FermBC< T, P, Q > & getFermBC() const
Return the fermion BC object for this action.
Definition: fermact.h:79
Class for counted reference semantics.
Definition: handle.h:33
void create(Handle< FermState< T, P, Q > > state)
Creation routine.
Definition: lwldslash_w.h:161
virtual void deriv(P &ds_u, const T &chi, const T &psi, enum PlusMinus isign) const
Take deriv of D.
unsigned long nFlops() const
Return flops performed by the operator()
Real H_parity
EXTERN Real KappaMD
int mu
Definition: cool.cc:24
int nu
Definition: cool.cc:25
Even-odd preconditioned Wilson fermion linear operator with parity breaking term.
void apply(T &chi, const T &psi, enum PlusMinus isign, int cb) const
General Wilson-Dirac dslash.
Definition: lwldslash_w.h:228
Nd
Definition: meslate.cc:74
Double tmp2
Definition: mesq.cc:30
Double tmp3
Definition: mesq.cc:31
LatticeReal phi
Definition: mesq.cc:27
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
dslash(u, eta, tmp, MINUS, 1)
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
START_CODE()
A(A, psi, r, Ncb, PLUS)
int cb
Definition: invbicg.cc:120
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
FloatingPoint< double > Double
Definition: gtest.h:7351
#define FORWARD
Definition: primitives.h:82
Real dummy
Definition: qtopcor.cc:36
multi1d< LatticeColorMatrix > deriv(const EvenOddPrecLinearOperator< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > &AP, const LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign)
Apply the operator onto a source std::vector.
#define ONE