CHROMA
iluprec_s_cprec_t_wilsonlike_linop_w.h
Go to the documentation of this file.
1 #ifndef ILUPREC_S_CPREC_T_WILSONLIKE_LINOP_W_H
2 #define ILUPREC_S_CPREC_T_WILSONLIKE_LINOP_W_H
3 
4 #include "qdp_config.h"
5 #if QDP_NS == 4
6 #if QDP_ND == 4
7 #if QDP_NC == 3
8 
9 #include "linearop.h"
10 #include "central_tprec_linop.h"
11 
12 #include "actions/ferm/linop/dslash_w.h"
13 #include "actions/ferm/linop/central_tprec_nospin_utils.h"
14 
15 namespace Chroma
16 {
17  //! Wilson Dirac Operator - Unpreconditioned in Space, Centrally Preconditioned in time
18  /*!
19  * \ingroup linop
20  *
21  * This routine is specific to Wilson fermions!
22  *
23  * ~ ~+
24  */
25 
26  class ILUPrecSCprecTWilsonLikeLinOp :
27  public ILUPrecSpaceCentralPrecTimeLinearOperator<LatticeFermion,
28  multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
29  {
30  public:
31  // Typedefs to save typing
32  typedef LatticeFermion T;
33  typedef multi1d<LatticeColorMatrix> P;
34  typedef multi1d<LatticeColorMatrix> Q;
35  typedef PScalar<PColorMatrix<RComplex<REAL>, 3> > CMat; // Useful type: ColorMat with no Outer<>
36  typedef PSpinVector<PColorVector<RComplex<REAL>, 3>, 2> HVec_site; // Useful type: Half Vec with no Outer<>
37 
38  //! Destructor is automatic
39  virtual ~ILUPrecSCprecTWilsonLikeLinOp() {}
40 
41 
42  protected:
43 
44 
45  virtual void spaceLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
46  virtual void derivSpaceOp(P& ds_u, const T& X, const T& Y, enum PlusMinus isign, int cb3d) const = 0;
47 
48  virtual void AMinusOneOp(T& chi, const T& psi, enum PlusMinus isign) const =0 ;
49  virtual void derivAMinusOne(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
50 
51  //! Get log det ( T^\dag T )
52  virtual Double logDetTDagT(void) const = 0;
53 
54  // Blanks to fill out...
55  virtual const multi3d< int >& getTsiteArray(void) const = 0;
56  virtual const multi3d< CMat >& getPMatrixArray(void) const = 0;
57  virtual const multi3d< CMat >& getPMatrixDaggerArray(void) const = 0;
58  virtual const multi2d< CMat >& getQMatrixInvArray(void) const = 0;
59  virtual const multi2d< CMat >& getQMatrixDaggerInvArray(void) const = 0;
60  virtual const Real& getFactor() const = 0;
61  virtual const Real& getInvFactor() const = 0;
62  virtual const multi1d< LatticeColorMatrix>& getLinks(void) const = 0;
63  virtual int getTMax(void) const = 0;
64  public:
65 
66  //! Return the fermion BC object for this linear operator
67  const FermBC<T,P,Q>& getFermBC() const = 0;
68 
69  //! The time direction
70  int tDir() const {
71  return 3;
72  }
73 
74 
75  //! Apply inv (C_L)^{-1}
76  virtual
77  void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const
78  {
79 
80  switch(isign) {
81  case PLUS:
82  {
83  //[ T+ D_s ] =
84  //[ 0 T+ ]
85  T tmp1;
86  TPlusOp(chi, psi, PLUS, 0);
87  TPlusOp(chi, psi, PLUS, 1);
88  spaceLinOp(tmp1, psi, PLUS, 0);
89  chi[rb3[0]] += tmp1;
90  }
91  break;
92  case MINUS:
93  {
94  //
95  // Remember to dagger in Even Odd Space too.
96  // [ (T+)^\dag 0 ]
97  // [ D_s^\dag (T+)^\dag ]
98  T tmp1;
99  TPlusOp(chi, psi, MINUS, 0);
100  TPlusOp(chi, psi, MINUS, 1);
101  spaceLinOp(tmp1, psi, MINUS, 1);
102  chi[rb3[1]] += tmp1;
103  }
104  break;
105  default:
106  QDPIO::cout << "Unknown sign " << std::endl;
107  QDP_abort(1);
108  }
109 
110  getFermBC().modifyF(chi);
111 
112  }
113 
114  //! Apply inv (C_R)^{-1}
115  virtual
116  void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
117  switch(isign) {
118  case PLUS:
119  {
120  //[ T- 0 ]
121  //[ D_s T- ]
122  T tmp1;
123  TMinusOp(chi, psi, PLUS, 0);
124  TMinusOp(chi, psi, PLUS, 1);
125  spaceLinOp(tmp1, psi, PLUS, 1);
126  chi[rb3[1]] += tmp1;
127  }
128  break;
129  case MINUS:
130  {
131  //
132  // Remember to dagger in Even Odd Space too.
133  // [ (T-)^\dag D_s^\dag ]
134  // [ 0 (T-)^\dag ]
135  T tmp1;
136  TMinusOp(chi, psi, MINUS, 0);
137  TMinusOp(chi, psi, MINUS, 1);
138  spaceLinOp(tmp1, psi, MINUS, 0);
139  chi[rb3[0]] += tmp1;
140  }
141  break;
142  default:
143  QDPIO::cout << "Unknown sign " << std::endl;
144  QDP_abort(1);
145  };
146  getFermBC().modifyF(chi);
147  }
148 
149  //! Apply C_L
150  virtual
151  void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
152 
153  switch(isign) {
154  case PLUS:
155  {
156  //[ (T+)^{-1} -(T+)^{-1}_{e} D_s_eo (T+)^{-1}_o ]
157  //[ 0 (T+)^{-1}_o ]
158  T tmp1, tmp2;
159 
160 
161  // tmp1_0 = (T+)^{-1}_{o} \psi_o
162  invTPlusOp(tmp1, psi, PLUS, 1);
163 
164  // chi_o = (T+)^{-1}_{o} \psi_o
165  chi[rb3[1]] = tmp1;
166 
167  // tmp2_e = D_s_eo (T+)^{-1}_{o} \psi_o
168  spaceLinOp(tmp2, tmp1, PLUS, 0);
169 
170  // tmp1_e = \psi_e - D_s_eo (T+)^{-1}_{o} \psi_o
171  tmp1[rb3[0]] = psi - tmp2;
172 
173 
174  // chi_e = (T+)^{-1}_{e} ( \psi_e - D_s_eo (T+)^{-1}_{o} \psi_o )
175  invTPlusOp(chi, tmp1, PLUS, 0);
176 
177  }
178  break;
179  case MINUS:
180  {
181  //[ (T+)^{-\dag 0 ]
182  //[ -(T+)^{-\dag}_{o} D_s_oe^\dag (T+)^{-\dag}_e (T+)^{-\dag}_o ]
183 
184  // = [ 1 0 ] [ T+^{-\dag} 0 ]
185  // [ 0 T+_o^{-\dag} ] [ -D_{oe}^\dag T+^{-\dag} 1 ]
186  T tmp1, tmp2;
187 
188  // tmp1 = T+^{-dag}_{e} \psi_{e}
189  invTPlusOp(tmp1, psi, MINUS, 0);
190 
191  // chi_e = T+^{-dag}_{e} \psi_{e}
192  chi[rb3[0]] = tmp1;
193 
194  // tmp2_o = D_{oe} T+^{-\dag} psi_e
195  spaceLinOp(tmp2, tmp1, MINUS, 1);
196 
197  // tmp1_o = psi_o - D_{oe} T+^{-\dag} psi_e
198  tmp1[rb3[1]] = psi - tmp2;
199 
200  invTPlusOp(chi, tmp1, MINUS, 1);
201  }
202  break;
203  default:
204  QDPIO::cout << "Unknown sign " << std::endl;
205  QDP_abort(1);
206  };
207  getFermBC().modifyF(chi);
208  }
209 
210  //! Apply C_R
211  virtual
212  void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
213 
214  switch(isign) {
215  case PLUS:
216  {
217  //[ (T-)^{-1}_e 0 ]
218  //[ -(T-)^{-1}_{o} D_s_oe^{-1} (T-)^{-1}_e (T-)^{-1}_o ]
219 
220  // = [ 1 0 ] [ T-^{-1}_e 0 ]
221  // [ 0 T-_o^{-1} ] [ -D_{oe}^\dag T-^{-1}_e 1 ]
222  T tmp1, tmp2;
223 
224  // tmp1 = T+^{-dag}_{e} \psi_{e}
225  invTMinusOp(tmp1, psi, PLUS, 0);
226 
227  // chi_e = T+^{-dag}_{e} \psi_{e}
228  chi[rb3[0]] = tmp1;
229 
230  // tmp2_o = D_{oe} T+^{-\dag} psi_e
231  spaceLinOp(tmp2, tmp1, PLUS, 1);
232 
233  // tmp1_o = psi_o - D_{oe} T+^{-\dag} psi_e
234  tmp1[rb3[1]] = psi - tmp2;
235 
236  invTMinusOp(chi, tmp1, PLUS, 1);
237  }
238 
239 
240  break;
241  case MINUS:
242  {
243  //[ (T-)^{-dag} -(T-)^{-\dag}_{e} D_s_eo^\dag (T-)^{-\dag}_o ]
244  //[ 0 (T-)^{-\dag}_o ]
245  T tmp1, tmp2;
246 
247 
248  // tmp1_0 = (T-)^{-\dag}_{o} \psi_o
249  invTMinusOp(tmp1, psi, MINUS, 1);
250 
251  // chi_o = (T-)^{-\dag}_{o} \psi_o
252  chi[rb3[1]] = tmp1;
253 
254  // tmp2_e = D_s_eo (T-)^{-1\dag_{o} \psi_o
255  spaceLinOp(tmp2, tmp1, MINUS, 0);
256 
257  // tmp1_e = \psi_e - D_s_eo (T-)^{-\dag}_{o} \psi_o
258  tmp1[rb3[0]] = psi - tmp2;
259 
260 
261  // chi_e = (T-)^{-\dag}_{e} ( \psi_e - D_s_eo (T-)^{-\dag}_{o} \psi_o )
262  invTMinusOp(chi, tmp1, MINUS, 0);
263 
264  }
265  break;
266  default:
267  QDPIO::cout << "Unknown sign " << std::endl;
268  QDP_abort(1);
269  };
270 
271  getFermBC().modifyF(chi);
272  }
273 
274  // A = 0 so A-1 = -1 => (A-1) psi = -psi and also for the dagger.
275 
276 
277  //! Flopcounter
278  virtual unsigned long nFlops() const
279  {
280  // QDPIO::cout << "Flopcount Not Yet Implemented " << std::endl;
281  return 0;
282  }
283 
284  protected:
285  virtual
286  void TPlusOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
287  {
288  LatticeHalfFermion tmp_minus;
289  LatticeHalfFermion tmp_plus;
290  LatticeHalfFermion tmp_T;
291 
292  // This is just a way of doing P+ psi - decompose and reconstruct.
293  // It may be more straightforward to just write a projector ... -- later
294  tmp_plus[ rb3[ cb3d ] ] = spinProjectDir3Plus(psi);
295  chi[ rb3[ cb3d ] ] = spinReconstructDir3Plus(tmp_plus);
296 
297  // Here I project - apply TOp to only the halfector - then reco nstruct
298  tmp_minus[ rb3[ cb3d ] ] = spinProjectDir3Minus(psi);
299 
300 
301  // Use shared routine to apply T or T^\dagger.
302  // Pass in u, and tsite for the subset
303  CentralTPrecNoSpinUtils::TOp(tmp_T,
304  tmp_minus,
305  getLinks(),
306  getTsiteArray()[cb3d],
307  getFactor(),
308  isign,
309  schroedingerTP());
310 
311  chi[rb3[ cb3d ]] += spinReconstructDir3Minus(tmp_T);
312  chi[rb3[ cb3d ]] *= Real(0.5);
313  // getFermBC().modifyF(chi);
314 
315 
316  }
317 
318  virtual
319  void TMinusOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
320  {
321  enum PlusMinus other_sign = (isign == PLUS ? MINUS : PLUS) ;
322  LatticeHalfFermion tmp_minus;
323  LatticeHalfFermion tmp_plus;
324  LatticeHalfFermion tmp_T;
325 
326  // This does the P- modulo a factor of 1/2
327  // Rather than spooling through 2 half vectors I could just do a
328  // ProjectDir3Minus rather than going through half vectirs
329  tmp_minus[rb3[cb3d]] = spinProjectDir3Minus(psi);
330  chi[rb3[cb3d]] = spinReconstructDir3Minus(tmp_minus);
331 
332  // This does the P+ T^\dagger modulo a factor of 1/2
333  tmp_plus[rb3[cb3d]] = spinProjectDir3Plus(psi);
334 
335 
336  // Use shared routine to apply T or T^\dagger.
337  // Pass in u, and tsite for the subset
338  CentralTPrecNoSpinUtils::TOp(tmp_T,
339  tmp_plus,
340  getLinks(),
341  getTsiteArray()[cb3d],
342  getFactor(),
343  other_sign,
344  schroedingerTP());
345 
346  chi[rb3[cb3d]] += spinReconstructDir3Plus(tmp_T);
347 
348  chi[rb3[cb3d]] *= Real(0.5); //The overall factor of 1/2
349 
350  // getFermBC().modifyF(chi);
351 
352 
353  }
354 
355  virtual
356  void invTPlusOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
357  {
358  LatticeHalfFermion tmp_minus;
359  LatticeHalfFermion tmp_plus;
360  LatticeHalfFermion tmp_T;
361 
362  tmp_plus[rb3[cb3d]] = spinProjectDir3Plus(psi);
363  chi[rb3[cb3d]] = spinReconstructDir3Plus(tmp_plus);
364 
365  tmp_minus[rb3[cb3d]] = spinProjectDir3Minus(psi);
366 
367  // Use shared routine to apply (T)^{-1} or (T^\dagger)^{-1}.
368  // Pass in u, and tsite, P, P^\dag, Q and Q^\dag for the subset.
369  CentralTPrecNoSpinUtils::invTOp( tmp_T,
370  tmp_minus,
371  getLinks(),
372  getTsiteArray()[cb3d],
373  getPMatrixArray()[cb3d],
374  getPMatrixDaggerArray()[cb3d],
375  getQMatrixInvArray()[cb3d],
376  getQMatrixDaggerInvArray()[cb3d],
377  getInvFactor(),
378  isign,
379  getTMax(),
380  schroedingerTP());
381 
382 
383  chi[rb3[cb3d]] += spinReconstructDir3Minus(tmp_T);
384  chi[rb3[cb3d]] *= Real(0.5);
385  // getFermBC().modifyF(chi);
386 
387 
388  }
389 
390  virtual
391  void invTMinusOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
392  {
393  enum PlusMinus other_sign = isign == PLUS ? MINUS : PLUS ;
394  LatticeHalfFermion tmp_minus;
395  LatticeHalfFermion tmp_plus;
396  LatticeHalfFermion tmp_T;
397 
398  tmp_minus[rb3[cb3d]] = spinProjectDir3Minus(psi);
399  chi[rb3[cb3d]] = spinReconstructDir3Minus(tmp_minus);
400 
401  tmp_plus[rb3[cb3d]] = spinProjectDir3Plus(psi);
402 
403  // Use shared routine to apply (T)^{-1} or (T^\dagger)^{-1}.
404  // Pass in u, and tsite, P, P^\dag, Q and Q^\dag for the subset.
405  CentralTPrecNoSpinUtils::invTOp( tmp_T,
406  tmp_plus,
407  getLinks(),
408  getTsiteArray()[cb3d],
409  getPMatrixArray()[cb3d],
410  getPMatrixDaggerArray()[cb3d],
411  getQMatrixInvArray()[cb3d],
412  getQMatrixDaggerInvArray()[cb3d],
413  getInvFactor(),
414  other_sign,
415  getTMax(),
416  schroedingerTP());
417 
418 
419  chi[rb3[cb3d]] += spinReconstructDir3Plus(tmp_T);
420  chi[rb3[cb3d]] *= Real(0.5);
421  // getFermBC().modifyF(chi);
422 
423 
424  }
425 
426  virtual
427  void derivInvTPlusOp(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const
428  {
429 
430  // Same code as C_L in unprec case
431  ds_u.resize(Nd);
432  for(int mu=0; mu < 3; mu++) {
433  ds_u[mu] = zero;
434  }
435 
436  if (isign == PLUS) {
437 
438  LatticeHalfFermion tmp1, tmp2;
439  T T1, T2;
440 
441  tmp1=spinProjectDir3Minus(Y);
442  for(int cb3d=0; cb3d < 2; cb3d++) {
443 
444  CentralTPrecNoSpinUtils::invTOp( tmp2,
445  tmp1,
446  getLinks(),
447  getTsiteArray()[cb3d],
448  getPMatrixArray()[cb3d],
449  getPMatrixDaggerArray()[cb3d],
450  getQMatrixInvArray()[cb3d],
451  getQMatrixDaggerInvArray()[cb3d],
452  getInvFactor(),
453  PLUS,
454  getTMax(),
455  schroedingerTP());
456  }
457 
458  T1 = spinReconstructDir3Minus(tmp2);
459 
460  tmp1=spinProjectDir3Minus(X);
461 
462  for(int cb3d=0; cb3d < 2; cb3d++) {
463 
464  CentralTPrecNoSpinUtils::invTOp( tmp2,
465  tmp1,
466  getLinks(),
467  getTsiteArray()[cb3d],
468  getPMatrixArray()[cb3d],
469  getPMatrixDaggerArray()[cb3d],
470  getQMatrixInvArray()[cb3d],
471  getQMatrixDaggerInvArray()[cb3d],
472  getInvFactor(),
473  MINUS,
474  getTMax(),
475  schroedingerTP());
476  }
477 
478  // Two factors of 0.5 from the projectors.
479  // Most cost efficient to apply them together to the half std::vector...
480  tmp2 *= Real(0.25);
481 
482  T2 = spinReconstructDir3Minus(tmp2);
483 
484 
485  LatticeFermion T3 = shift(T1, FORWARD, 3);
486 
487  // A minus from the fact that its dT^{-1}
488  // A minus from the fact that the shift has a -
489  // Makes a +
490  ds_u[3] = traceSpin(outerProduct(T3, T2));
491 
492  }
493  else {
494  ds_u[3] = zero;
495  }
496 
497  getFermBC().zero(ds_u);
498  }
499 
500  virtual
501  void derivInvTMinusOp(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const
502  {
503  // Same code as derivC_R
504  ds_u.resize(Nd);
505  for(int mu=0; mu < 3; mu++) {
506  ds_u[mu] = zero;
507  }
508 
509  if (isign == MINUS) {
510 
511  LatticeHalfFermion tmp1, tmp2;
512  T T1, T2;
513 
514  tmp1=spinProjectDir3Plus(Y);
515  for(int cb3d=0; cb3d < 2; cb3d++) {
516  CentralTPrecNoSpinUtils::invTOp( tmp2,
517  tmp1,
518  getLinks(),
519  getTsiteArray()[cb3d],
520  getPMatrixArray()[cb3d],
521  getPMatrixDaggerArray()[cb3d],
522  getQMatrixInvArray()[cb3d],
523  getQMatrixDaggerInvArray()[cb3d],
524  getInvFactor(),
525  PLUS,
526  getTMax(),
527  schroedingerTP());
528  }
529  T1 = spinReconstructDir3Plus(tmp2);
530 
531  tmp1=spinProjectDir3Plus(X);
532  for(int cb3d=0; cb3d < 2; cb3d++) {
533 
534  CentralTPrecNoSpinUtils::invTOp( tmp2,
535  tmp1,
536  getLinks(),
537  getTsiteArray()[cb3d],
538  getPMatrixArray()[cb3d],
539  getPMatrixDaggerArray()[cb3d],
540  getQMatrixInvArray()[cb3d],
541  getQMatrixDaggerInvArray()[cb3d],
542  getInvFactor(),
543  MINUS,
544  getTMax(),
545  schroedingerTP());
546  }
547  // Two factors of 0.5 from the projectors.
548  // Most cost efficient to apply them together to the half std::vector...
549  tmp2 *= Real(0.25);
550 
551  T2 = spinReconstructDir3Plus(tmp2);
552 
553  LatticeFermion T3 = shift(T1, FORWARD, 3);
554 
555  // A minus from the fact that its dT^{-1}
556  // A minus from the fact that the shift has a -
557  // Makes a +
558  ds_u[3] = traceSpin(outerProduct(T3, T2));
559 
560  }
561  else {
562  ds_u[3]= zero;
563  }
564 
565  getFermBC().zero(ds_u);
566  }
567 
568  virtual
569  void derivCRight(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const
570  {
571  T T_1, T_2, T_3, T_4;
572  switch(isign) {
573  case PLUS:
574  {
575  // Simple Case: -T^o^\dagger_2 \dot(Dslash^oe_S) T^e_1
576 
577  // [ T^e_1 ] [ (T-)^{-1}e 0 ] [ Y_e ] [ (T-)^{-1}_e Y_e ]
578  // [ ] = [ ] [ ] = [ ]
579  // { T^o_1 ] [ 0 1 ] [ Y_o ] [ Y_o ]
580  invTMinusOp(T_1, Y, PLUS, 0);
581  T_1[ rb3[1] ] = Y;
582 
583  // [ T^e_2 ] [ 1 0 ] [ X_e ] [ X_e ]
584  // [ ] = [ ] [ ] =[ ]
585  // [ T^o_2 ] [ 0 (T-)^{-dag}o ] [ X_o ] [ (T-)^{-dag}_o X_o ]
586  invTMinusOp(T_2, X, MINUS, 1);
587  T_2[ rb3[0] ] = X;
588 
589  derivSpaceOp(ds_u, T_2, T_1, PLUS, 1);
590  for(int mu=0; mu < Nd; mu++) {
591  ds_u[mu] *= Real(-1);
592  }
593  }
594  break;
595  case MINUS:
596  {
597  // [ T^e_1 ] [ 1 0 ] [ Y_e ] [ Y_e ]
598  // [ ] = [ ] [ ] =[ ]
599  // [ T^o_1 ] [ 0 (T-)^{-dag}o ] [ Y_o ] [ (T-)^{-dag}_o Y_o ]
600  invTMinusOp(T_1, Y, MINUS, 1);
601  T_1[ rb3[0] ] = Y;
602 
603 
604  // [ T^e_2 ] [ (T-)^{-1}e 0 ] [ X_e ] [ (T-)^{-1}_e X_e ]
605  // [ ] = [ ] [ ] = [ ]
606  // { T^o_2 ] [ 0 1 ] [ X_o ] [ X_o ]
607  invTMinusOp(T_2, X, PLUS, 0);
608  T_2[ rb3[1] ] = X;
609 
610  // [ T^e_3 ] [ 1 - D^{oe}^\dag_s ][ T^e_1 ]
611  // [ ] = [ ][ ]
612  // [ T^o_3 ] [ 0 1 ][ T^o_1 ]
613  // We only need the even part tho, and put Y_o in the bottom
614  T T_tmp;
615  spaceLinOp(T_tmp, T_1, MINUS, 0);
616  T_3[rb3[0]] = T_1 - T_tmp; // Mhalf from the (-1/2)dslash
617  T_3[rb3[1]] = Y;
618 
619 
620  // [ T^e_4 ] [ 1 0 ][ T^e_2 ]
621  // [ ] = [ ][ ]
622  // [ T^o_4 ] [ -D^{eo}_s 1 ][ T^o_2 ]
623  // We only need the odd part tho, and put X_o in the top
624  spaceLinOp(T_tmp, T_2, PLUS, 1);
625  T_4[rb3[1]] = T_2 - T_tmp; // Mhalf from the (-1/2)dslash
626  T_4[rb3[0]] = X;
627 
628 
629  derivInvTMinusOp(ds_u, T_4, T_3, MINUS);
630 
631  P ds_tmp;
632 
633  derivSpaceOp(ds_tmp, T_2, T_1, MINUS, 0);
634  for(int mu=0; mu < Nd; mu++) {
635  ds_u[mu] -= ds_tmp[mu];
636  }
637 
638  }
639  break;
640  default:
641  QDPIO::cerr << "Bad Case. Should never get here" << std::endl;
642  QDP_abort(1);
643  }
644  getFermBC().zero(ds_u);
645  }
646 
647  virtual
648  void derivCLeft(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const
649  {
650  T T_1, T_2, T_3, T_4;
651 
652  switch(isign) {
653  case PLUS:
654  {
655  // [ T^e_1 ] [ 1 0 ] [ Y_e ] [ Y_e ]
656  // [ ] = [ ] [ ] =[ ]
657  // [ T^o_1 ] [ 0 (T+)^{-1}_o ] [ Y_o ] [ (T+)^{-1}_o Y_o ]
658  invTPlusOp(T_1, Y, PLUS, 1);
659  T_1[ rb3[0] ] = Y;
660 
661 
662  // [ T^e_2 ] [ (T+)^{-dag}e 0 ] [ X_e ] [ (T+)^{-dag}_e X_e ]
663  // [ ] = [ ] [ ] = [ ]
664  // { T^o_2 ] [ 0 1 ] [ X_o ] [ X_o ]
665  invTPlusOp(T_2, X, MINUS, 0);
666  T_2[ rb3[1] ] = X;
667 
668  // [ T^e_3 ] [ 1 - D^{eo}_s ][ T^e_1 ]
669  // [ ] = [ ][ ]
670  // [ T^o_3 ] [ 0 1 ][ T^o_1 ]
671  // We only need the even part tho, and put Y_o in the bottom
672  T T_tmp;
673  spaceLinOp(T_tmp, T_1, PLUS, 0);
674  T_3[rb3[0]] = T_1 - T_tmp; // mhalf is from (-1/2)Dslash
675  T_3[rb3[1]] = Y;
676 
677 
678  // [ T^e_4 ] [ 1 0 ][ T^e_2 ]
679  // [ ] = [ ][ ]
680  // [ T^o_4 ] [ -D^{oe \dag}_s 1 ][ T^o_2 ]
681  // We only need the odd part tho, and put X_o in the top
682  spaceLinOp(T_tmp, T_2, MINUS, 1);
683  T_4[rb3[1]] = T_2 - T_tmp; // mhalf is from (-1/2) Dslash
684  T_4[rb3[0]] = X;
685 
686  P ds_tmp;
687  derivSpaceOp(ds_tmp, T_2, T_1, PLUS, 0);
688  derivInvTPlusOp(ds_u, T_4, T_3, PLUS);
689  for(int mu=0; mu < Nd; mu++) {
690  ds_u[mu] -= ds_tmp[mu];
691  }
692  }
693  break;
694 
695  case MINUS:
696  {
697  // Simple Case: -T^o^\dagger_2 \dot(Dslash^oe_S) T^e_1
698 
699  // [ T^e_1 ] [ (T+)^{-\dag}e 0 ] [ Y_e ] [ (T+)^{-\dag}_e Y_e ]
700  // [ ] = [ ] [ ] = [ ]
701  // { T^o_1 ] [ 0 1 ] [ Y_o ] [ Y_o ]
702  invTPlusOp(T_1, Y, MINUS, 0);
703  T_1[ rb3[1] ] = Y;
704 
705  // [ T^e_2 ] [ 1 0 ] [ X_e ] [ X_e ]
706  // [ ] = [ ] [ ] =[ ]
707  // [ T^o_2 ] [ 0 (T+)^{-1}o ] [ X_o ] [ (T+)^{-1}_o X_o ]
708  invTPlusOp(T_2, X, PLUS, 1);
709  T_2[ rb3[0] ] = X;
710 
711  derivSpaceOp(ds_u, T_2, T_1, MINUS, 1);
712  for(int mu=0; mu < Nd; mu++) {
713  ds_u[mu] *= Real(-1); // (Combinme factor of -1 and -(1/2)
714  }
715  }
716  break;
717  default:
718  QDPIO::cerr << "Bad Case. Should never get here" << std::endl;
719  QDP_abort(1);
720  }
721 
722  getFermBC().zero(ds_u);
723  }
724 
725  //! Get the force due to the det T^\dag T bit
726  // This is the same as for the unprec case as it turns
727  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const
728  {
729  // Derivative of a Hermitian quantity so ignore isign?
730  // Initial development -- set to 0
731  ds_u.resize(Nd);
732 
733 
734  for(int mu=0; mu < Nd; mu++) {
735  ds_u[mu] = zero;
736  }
737 
738  if( !schroedingerTP() ) {
739 
740  // Force from even sites
741  CentralTPrecNoSpinUtils::derivLogDet(ds_u,
742  getLinks(),
743  getQMatrixInvArray(),
744  getTsiteArray(),
745  tDir(),
746  getFactor(),
747  schroedingerTP());
748 
749 
750 
751  getFermBC().zero(ds_u);
752  }
753  else {
754  ds_u[tDir()] = zero;
755  }
756  }
757 
758  };
759 
760 
761 
762 } // End Namespace Chroma
763 
764 #endif
765 
766 #endif
767 #endif
768 #endif
central_tprec_linop.h
Time-preconditioned Linear Operators.
central_tprec_nospin_utils.h
Support for time preconditioning.
mu
int mu
Definition: cool.cc:24
dslash_w.h
Include possibly optimized Wilson dslash.
linearop.h
Linear Operators.
Nd
Nd
Definition: meslate.cc:74
tmp2
Double tmp2
Definition: mesq.cc:30
Chroma::InlineDiscoEigCGEnv::P
multi1d< LatticeColorMatrix > P
Definition: inline_disco_eigcg_w.cc:315
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::tDir
static int tDir()
Definition: coulgauge.cc:14
Chroma::Q
LinOpSysSolverMGProtoClover::Q Q
Definition: syssolver_linop_clover_mg_proto.cc:64
Chroma::T
LinOpSysSolverMGProtoClover::T T
Definition: syssolver_linop_clover_mg_proto.cc:63
Chroma::MINUS
@ MINUS
Definition: chromabase.h:45
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::chi
multi1d< LatticeFermion > chi(Ncb)
Chroma::psi
LatticeFermion psi
Definition: mespbg5p_w.cc:35
Chroma::zero
Double zero
Definition: invbicg.cc:106
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
FORWARD
#define FORWARD
Definition: primitives.h:82
T
LatticeFermion T
Definition: t_clover.cc:11
P
multi1d< LatticeColorMatrix > P
Definition: t_clover.cc:13
Generated by doxygen 1.9.1