CHROMA
eoprec_ovlap_contfrac5d_linop_base_array_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Base class for even-odd prec. 5D continued fraction linop
3  */
4 
7 
8 using namespace QDP::Hints;
9 
10 namespace Chroma
11 {
12  EvenOddPrecOvlapContFrac5DLinOpBaseArray::EvenOddPrecOvlapContFrac5DLinOpBaseArray(
14  const Real& _m_q,
15  const Real& _OverMass,
16  int _N5,
17  const Real& _scale_fac,
18  const multi1d<Real>& _alpha,
19  const multi1d<Real>& _beta,
20  const bool _isLastZeroP ) :
21  m_q(_m_q), OverMass(_OverMass), N5(_N5), scale_fac(_scale_fac),
22  alpha(_alpha), beta(_beta), isLastZeroP(_isLastZeroP)
23  {
24  START_CODE();
25 
26  int dslash_length;
27  if( isLastZeroP ) {
28  dslash_length = N5-1;
29  }
30  else {
31  dslash_length = N5;
32  }
33 
34  Dslash.create(state,dslash_length);
35 
36  // The mass ratio
37  Real mass = ( Real(1) + m_q ) / (Real(1) - m_q);
38 
39  // Now compute some coefficients.
40  // First the beta_tilde_i
41  // Basically this is beta[n]*Hsign*scale_fac
42  // Now N5 is always odd. So the first Hsign is +
43  // and the last one should also be
44  // Hence at the end of this loop Hsign should be flipped from +->-
45  beta_tilde.resize(N5);
46  int Hsign = 1;
47  for(int i=0; i < N5; i++)
48  {
49  // Flip Hsign
50  beta_tilde[i] = beta[i]*Hsign*scale_fac;
51 
52  /*
53  QDPIO::cout << "beta["<<i<<"]=" << beta[i]
54  << " Hsign=" << Hsign
55  << " scale_fac=" << scale_fac
56  << " beta_tilde["<<i<<"]=" << beta_tilde[i] << std::endl;
57 
58  */
59  Hsign = -Hsign;
60  }
61 
62  // Sanity Check
63  if ( Hsign > 0 ) {
64  QDPIO::cerr << "Something is wrong. At the end of this loop"
65  << " Hsign should be -ve" << std::endl;
66  }
67 
68  // Now the a_i's and b_i's
69  a.resize(N5);
70  for(int i=0; i < N5-1; i++) {
71  a[i] = beta_tilde[i]*(Nd - OverMass);
72  }
73  a[N5-1] = mass + (beta_tilde[N5-1]*(Nd - OverMass));
74 
75  /*
76  QDPIO::cout << "Nd - OverMass = " << Nd- OverMass << std::endl;
77  for(int i=0; i < N5; i++) {
78  QDPIO::cout << "a["<<i<<"]= " << a[i] << std::endl;
79  }
80  */
81  // Now the d-s
82  multi1d<Real> d(N5);
83  invd.resize(N5);
84 
85  d[0] = a[0];
86  invd[0] = Real(1)/d[0];
87 
88  for(int i=1; i < N5; i++) {
89  d[i] = a[i] - (alpha[i-1]*alpha[i-1])/d[i-1];
90  invd[i] = Real(1)/d[i];
91  }
92 
93 
94  /*
95  for(int i=0; i < N5; i++) {
96  QDPIO::cout << "d["<<i<<"]=" << d[i] << std::endl;
97  }
98  */
99 
100  // Now the u-s
101  u.resize(N5-1);
102  for(int i=0; i < N5-1; i++) {
103  u[i] = alpha[i]/d[i];
104  }
105 
106  off_diag_coeff.resize(N5);
107  for(int i=0; i < N5; i++) {
108  off_diag_coeff[i] = -Real(0.5)*beta_tilde[i];
109  }
110  /*
111  for(int i=0; i < N5-1; i++) {
112  QDPIO::cout << "u["<<i<<"] = " << u[i] << std::endl;
113  }
114  */
115  END_CODE();
116  }
117 
118 
119 
120 
121  //! Apply the even-even (odd-odd) coupling piece of the domain-wall fermion operator
122  /*!
123  * \ingroup linop
124  *
125  * The operator acts on the entire lattice
126  *
127  * \param psi Pseudofermion field (Read)
128  * \param isign Flag ( PLUS | MINUS ) (Read)
129  * \param cb checkerboard ( 0 | 1 ) (Read)
130  *
131  *
132  * Flopcount: N5*6NcNs + (N5-2)*4NcNs = NcNs( 6N5 +4(N5-2) ) = (10N5-8) Nc Ns / cb_site
133  */
134  void
136  multi1d<LatticeFermion>& chi,
137  const multi1d<LatticeFermion>& psi,
138  enum PlusMinus isign,
139  const int cb) const
140  {
141  START_CODE();
142 
143  if( chi.size() != N5 ) chi.resize(N5);
144 
145  // We don't care about the isign because our operator is Hermitian
146  // Apply matrix
147  // [ A_0 B_0 0 ... ] [ psi_0 ]
148  // [ B_0 A_1 B_1 ... ... ] [ psi_1 ]
149  // [ 0 ... ... ... ... ] [ psi_2 ]
150  // [ ... ... 0 B_N5-3 A_N5-2 B_N5-2 ] [ psi_N5-2 ]
151  // [ ... ... ... 0 B_N5-2 A_N5-1 ] [ psi_N5-1 ]
152 
153  // With A_i = gamma_5 a_i = a_i gamma_5
154  // and B_i = b_i I = alpha_i I
155 
156  LatticeFermion tmp; moveToFastMemoryHint(tmp);
157  int G5=Ns*Ns-1;
158 
159  // First 0ne
160  // Operation: chi[0][rb[cb]] = a[0] G5 psi[0] + alpha[0]*psi[1]
161  //
162  // Useful flops: 6Nc Ns / site
163  // tmp[rb[cb]] = Gamma(G5)*psi[0];
164  // chi[0][rb[cb]] = a[0]*tmp;
165 
166  if( N5 > 1 ) {
167  chi[0][rb[cb]] = alpha[0]*psi[1] + a[0]*(GammaConst<Ns,Ns*Ns-1>()*psi[0]);
168  }
169  else {
170  chi[0][rb[cb]] = a[0]*(GammaConst<Ns,Ns*Ns-1>()*psi[0]);
171  }
172 
173  // All the rest
174  for(int i=1; i < N5; i++) {
175 
176  // Operation:
177  // N5 - 1 times:
178  // chi[i] = alpha[i-1]*psi[i-1] + a[i] Gamma_5 *psi[i]
179  // N5 - 2 times:
180  // chi[i] += alpha[i]*psi[i+1];
181  // Useful flops = (N5-1) * 6NcNs + (N5-2)*4Nc*Ns
182 
183  /*
184  // B_{i-1} psi_[i-1]
185  chi[i][rb[cb]] = alpha[i-1]*psi[i-1];
186 
187  // A_{i} psi[i] = a_{i} g_5 psi[i]
188  tmp[rb[cb]] = Gamma(G5)*psi[i];
189  chi[i][rb[cb]] += a[i]*tmp;
190  */
191  chi[i][rb[cb]] = alpha[i-1]*psi[i-1] + a[i]*(GammaConst<Ns,Ns*Ns-1>()*psi[i]);
192 
193  // When i hits N5-1, we don't have the B_N5-1 term
194  if(i < N5-1) {
195  chi[i][rb[cb]] += alpha[i]*psi[i+1];
196  }
197  }
198 
199  END_CODE();
200  }
201 
202 
203  //! Apply the inverse even-even (odd-odd)
204  /*!
205  * \ingroup linop
206  *
207  * Here we apply the LDU decomposition
208  *
209  * \param psi Pseudofermion field (Read)
210  * \param isign Flag ( PLUS | MINUS ) (Read)
211  * \param cb checkerboard ( 0 | 1 ) (Read)
212 
213  *
214  * Total flopcount: (N5-1)*4NcNs + 2NcNs + (N5-1)*6NcNs
215  * = (N5-1)*10NcNs + 2NcNs
216  * = (10N5-8) Nc Ns per_cb_site
217  */
218  void
220  multi1d<LatticeFermion>& chi,
221  const multi1d<LatticeFermion>& psi,
222  enum PlusMinus isign,
223  const int cb) const
224  {
225  START_CODE();
226 
227  if( chi.size() != N5 ) chi.resize(N5);
228 
229  multi1d<LatticeFermion> y(N5); moveToFastMemoryHint(y);
230 
231  Real coeff;
232 
233  const int G5 = Ns*Ns-1;
234 
235  // Solve L y = psi
236  y[0][rb[cb]] = psi[0];
237 
238  // (N5-1)*4NcNs
239  for(int i = 1; i < N5; i++) {
240  // tmp[rb[cb]] = Gamma(G5)*y[i-1];
241  // y[i][rb[cb]] = psi[i] - u[i-1]*tmp;
242  y[i][rb[cb]] = psi[i] - u[i-1]*(GammaConst<Ns,Ns*Ns-1>()*y[i-1]);
243  }
244 
245  // Invert diagonal piece y <- D^{-1} y
246  // N5 times: y = (1/d_i) gamma_5 y[i]
247  // Useful flops: N5 * 2NcNs flops / site
248  // Rolled this into the bottom loop
249 
250  // Backsubstitute U chi = y
251 
252  // 2NcNs
253 
254  chi[N5-1][rb[cb]] = invd[N5-1]*(GammaConst<Ns,Ns*Ns-1>()*y[N5-1]);
255 
256  // N5-1 * 6NcNs
257  for(int i = N5-2; i >= 0; i--) {
258  // tmp[rb[cb]] = Gamma(G5)*chi[i+1]
259  // chi[i][rb[cb]] = y[i] - u[i]*tmp;
260  // y[i][rb[cb]] = invd[i]*(GammaConst<Ns,Ns*Ns-1>()*y[i]);
261  // chi[i][rb[cb]] = y[i] - u[i]*(GammaConst<Ns,Ns*Ns-1>()*chi[i+1]);
262  chi[i][rb[cb]] = GammaConst<Ns,Ns*Ns-1>()*(invd[i]*y[i]-u[i]*chi[i+1]);
263  }
264 
265  //Done! That was not that bad after all....
266  //See, I told you so...
267  END_CODE();
268  }
269 
270  //! Apply the off diagonal block
271  /*!
272  * \param chi result (Modify)
273  * \param psi source (Read)
274  * \param isign Flag ( PLUS | MINUS ) (Read)
275  * \param cb checkerboard ( 0 | 1 ) (Read)
276  *
277  * Total flopcount: (N5-1)*(Dslash_cb_flops + 2NcNs) per cb_site (isLastZeroP==true)
278  * N5*(Dslash_cb_flops+2NcNs) per cb_site (isLastZeroP==false)
279  *
280  *
281  */
283  multi1d<LatticeFermion>& chi,
284  const multi1d<LatticeFermion>& psi,
285  enum PlusMinus isign,
286  const int cb) const
287  {
288  START_CODE();
289 
290  if( chi.size() != N5 ) chi.resize(N5);
291 
292  multi1d<LatticeFermion> tmp(N5); moveToFastMemoryHint(tmp);
293 
294  int G5 = Ns*Ns-1;
295 
296  // Optimisation... do up to the penultimate block...
297 
298  // (N5-1)( Dslash + 2NcNs) flops/site
299  Dslash.apply(tmp, psi, PLUS, cb);
300 
301  for(int i=0; i < N5-1; i++) {
302  /*
303  // CB is CB of TARGET
304  // gamma_5 Dslash is hermitian so I can ignore isign
305 
306  // Apply g5 Dslash
307  Dslash.apply(tmp, psi[i], PLUS, cb);
308 
309  // chi[i][rb[cb]] = Gamma(G5)*tmp;
310  */
311  // Chi_i is now -(1/2) beta_tilde_i Dslash
312  chi[i][rb[cb]] = off_diag_coeff[i]*(GammaConst<Ns,Ns*Ns-1>()*tmp[i]);
313  }
314 
315 
316  // Only do last block if beta_tilde[i] is not zero
317  // ( Dslash + 2NcNs) flops per site if done.
318  if( !isLastZeroP ) {
319 
320  // Vector Dslash gets this appropriately
321  //Dslash.apply(tmp, psi[N5-1], PLUS, cb);
322  // chi[N5-1][rb[cb]] = Gamma(G5)*tmp;
323 
324  // Chi_i is now -(1/2) beta_tilde_i Dslash
325  chi[N5-1][rb[cb]] = off_diag_coeff[N5-1]*(GammaConst<Ns,Ns*Ns-1>()*tmp[N5-1]);
326  }
327  else {
328  chi[N5-1][rb[cb]] = zero;
329  }
330 
331  END_CODE();
332  }
333 
334 
335  // Derivative of even-odd linop component
336  /*
337  * This is a copy of the above applyOffDiag with the D.apply(...) replaced
338  * by D.deriv(ds_tmp,...) like calls.
339  */
340  void
342  multi1d<LatticeColorMatrix>& ds_u,
343  const multi1d<LatticeFermion>& chi,
344  const multi1d<LatticeFermion>& psi,
345  enum PlusMinus isign,
346  int cb) const
347  {
348  START_CODE();
349 
350  ds_u.resize(Nd);
351  ds_u = zero;
352 
353  multi1d<LatticeColorMatrix> ds_tmp(Nd);
354 
355  LatticeFermion tmp; moveToFastMemoryHint(tmp);
356 
357  Real coeff;
358  int G5 = Ns*Ns-1;
359 
360  switch (isign)
361  {
362  case PLUS:
363  // Optimisation... do up to the penultimate block...
364  for(int i=0; i < N5; i++)
365  {
366  if (i == N5-1 && isLastZeroP) continue;
367 
368  // CB is CB of TARGET
369  // consider case of gamma_5 Dslash
370  tmp[rb[cb]] = Gamma(G5)*chi[i];
371 
372  // Multiply coefficient
373  coeff = -Real(0.5)*beta_tilde[i];
374 
375  // Chi_i is now -(1/2) beta_tilde_i Dslash
376  tmp[rb[cb]] *= coeff;
377 
378  // Apply g5 Dslash
379  Dslash.deriv(ds_tmp, tmp, psi[i], PLUS, cb);
380  ds_u += ds_tmp;
381  }
382  break;
383 
384  case MINUS:
385  // Optimisation... do up to the penultimate block...
386  for(int i=0; i < N5; i++)
387  {
388  if (i == N5-1 && isLastZeroP) continue;
389 
390  // CB is CB of TARGET
391  // consider case of Dslash^dag gamma_5
392  tmp[rb[1-cb]] = Gamma(G5)*psi[i];
393 
394  // Multiply coefficient
395  coeff = -Real(0.5)*beta_tilde[i];
396 
397  // Chi_i is now -(1/2) beta_tilde_i Dslash
398  tmp[rb[1-cb]] *= coeff;
399 
400  // Apply g5 Dslash
401  Dslash.deriv(ds_tmp, chi[i], tmp, MINUS, cb);
402  ds_u += ds_tmp;
403  }
404  break;
405 
406  default:
407  QDP_error_exit("unknown case");
408  }
409 
410  END_CODE();
411  }
412 
413 } // End Namespace Chroma
virtual void applyOffDiag(multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign, const int cb) const
Apply the off diagonal block.
virtual void applyDerivOffDiag(multi1d< LatticeColorMatrix > &ds_u, const multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign, int cb) const
Apply the even-odd (odd-even) coupling piece of the NEF operator.
Support class for fermion actions and linear operators.
Definition: state.h:94
Class for counted reference semantics.
Definition: handle.h:33
virtual void deriv(multi1d< LatticeColorMatrix > &ds_u, const multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign) const
Take deriv of D.
Include possibly optimized Wilson dslash.
Base class for Even-odd prec. 5D continued fraction linop.
virtual 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.
void apply(multi1d< LatticeFermion > &chi, const multi1d< LatticeFermion > &psi, enum PlusMinus isign, int cb) const
General Wilson-Dirac dslash.
void create(Handle< FermState< T, P, Q > > state, int N5_)
Creation routine.
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.
int y
Definition: meslate.cc:35
Nd
Definition: meslate.cc:74
static const LatticeInteger & beta(const int dim)
Definition: stag_phases_s.h:47
static const LatticeInteger & alpha(const int dim)
Definition: stag_phases_s.h:43
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
int G5
Definition: pbg5p_w.cc:57
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
int i
Definition: pbg5p_w.cc:55
Double mass
Definition: pbg5p_w.cc:54
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
DComplex d
Definition: invbicg.cc:99
START_CODE()
int cb
Definition: invbicg.cc:120
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
Double zero
Definition: invbicg.cc:106