CHROMA
central_tprec_linop.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! @file
3  * @brief Time-preconditioned Linear Operators
4  */
5 
6 #ifndef central_tprec_linop_h
7 #define central_tprec_linop_h
8 
9 #include "qdp_config.h"
10 
11 #if QDP_NS == 4
12 #if QDP_NC == 3
13 #if QDP_ND == 4
14 
15 #include "linearop.h"
16 #include "actions/ferm/fermbcs/schroedinger_fermbc_w.h"
17 #include <typeinfo>
18 
19 namespace Chroma
20 {
21 
22  //! Central Preconditioned Linear Operators
23 
24  //! For now, no forces just yet - come later...
25  template<typename T, typename P, typename Q>
26  class CentralTimePrecLinearOperator : public DiffLinearOperator<T,P,Q>
27  {
28  public:
29  //! Virtual destructor to help with cleanup;
30  virtual ~CentralTimePrecLinearOperator() {}
31 
32  //! Defined on the entire lattice
33  const Subset& subset() const = 0;
34 
35  //! Return the fermion BC object for this linear operator
36  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
37 
38 
39  //! Do we have SchroedingerBCs in time?
40  //! Yucky but doable as a default
41  virtual bool schroedingerTP() const
42  {
43 
44  // Get the BCs
45  const FermBC<T,P,Q>& fbc=getFermBC();
46 
47  // If they are nontrivial
48  if( fbc.nontrivialP() ) {
49 
50  // Try and cast to a SchrFermBC base class
51  try {
52  const SchrFermBC& schrReference = dynamic_cast<const SchrFermBC&>(fbc);
53  // Success -- check whether its dir is the same as my tDir
54  return ( schrReference.getDir() == tDir() );
55  }
56  catch( std::bad_cast ) {
57  // Cast failed - so not Schroedinger(?)
58  return false;
59  }
60  }
61 
62  // Wasn't nontrivial to start with.
63  return false;
64  }
65 
66 
67  //! The time direction
68  virtual int tDir() const = 0;
69 
70  //! Apply inv (C_L)^{-1}
71  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
72 
73  //! Apply inv (C_R)^{-1}
74  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
75 
76  //! Apply C_L
77  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
78 
79  //! Apply C_R
80  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
81 
82 
83  //! Apply the operator onto a source std::vector
84  // This varies a lot amongst the families
85  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const = 0;
86 
87 
88  //! Apply the UNPRECONDITIONED operator onto a source std::vector
89  /*! Mainly intended for debugging */
90  /*! This by the way defines the essence of the preconditioning ie:
91  * M_unprec = C^{-1}_L M_prec C^{-1}_R
92  *
93  * This applies whether we use
94  * no spatial preconditioning: M_prec = 1 + C_L D_s C_R
95  * ILU preconditioning M_prec = L + U + L (A-1) U
96  * Here L & U involve both D_s and C_L & C_R, with
97  * A being a left over diagonal piece.
98  * However det(L) = det(C_R), det(U) = det(C_L)
99  * EO3DPrec M_prec = L D U
100  * here LDU is a Schur Decomposition of M_prec from no
101  * spatial decomposition case.
102  */
103  virtual void unprecLinOp(T& chi, const T& psi,
104  enum PlusMinus isign) const
105  {
106 
107  switch (isign) {
108  case PLUS:
109  {
110  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
111 
112  invCRightLinOp(tmp1, psi, isign);
113  (*this)(tmp2, tmp1, isign);
114  invCLeftLinOp(chi, tmp2, isign);
115  }
116  break;
117  case MINUS:
118  {
119  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
120 
121  invCLeftLinOp(tmp1, psi, isign);
122  (*this)(tmp2, tmp1, isign);
123  invCRightLinOp(chi, tmp2, isign);
124  }
125  break;
126  default:
127  QDPIO::cerr << "unknown sign" << std::endl;
128  QDP_abort(1);
129  }
130 
131 
132  getFermBC().modifyF(chi);
133  }
134 
135 
136  //! Apply the d/dt of the preconditioned linop
137  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
138 
139  //! Get log det ( T^\dag T )
140  virtual Double logDetTDagT(void) const = 0;
141 
142  //! Get the force due to the det T^\dag T bit
143  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const = 0;
144  };
145 
146 
147  template<typename T, typename P, typename Q>
148  class Central2TimePrecLinearOperator : public DiffLinearOperator<T,P,Q>
149  {
150  public:
151  //! Virtual destructor to help with cleanup;
152  virtual ~Central2TimePrecLinearOperator() {}
153 
154  //! Defined on the entire lattice
155  const Subset& subset() const = 0;
156 
157  //! Return the fermion BC object for this linear operator
158  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
159 
160 
161  //! Do we have SchroedingerBCs in time?
162  //! Yucky but doable as a default
163  virtual bool schroedingerTP() const
164  {
165 
166  // Get the BCs
167  const FermBC<T,P,Q>& fbc=getFermBC();
168 
169  // If they are nontrivial
170  if( fbc.nontrivialP() ) {
171 
172  // Try and cast to a SchrFermBC base class
173  try {
174  const SchrFermBC& schrReference = dynamic_cast<const SchrFermBC&>(fbc);
175  // Success -- check whether its dir is the same as my tDir
176  return ( schrReference.getDir() == tDir() );
177  }
178  catch( std::bad_cast ) {
179  // Cast failed - so not Schroedinger(?)
180  return false;
181  }
182  }
183 
184  // Wasn't nontrivial to start with.
185  return false;
186  }
187 
188 
189  //! The time direction
190  virtual int tDir() const = 0;
191 
192  //! Apply inv (C_L)^{-1}
193  virtual void invLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
194 
195  //! Apply inv (C_R)^{-1}
196  virtual void invRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
197 
198  //! Apply C_L
199  virtual void leftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
200 
201  //! Apply C_R
202  virtual void rightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
203 
204 
205  //! Apply the operator onto a source std::vector
206  // This varies a lot amongst the families
207  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const = 0;
208 
209 
210  //! Apply the UNPRECONDITIONED operator onto a source std::vector
211  /*! Mainly intended for debugging */
212  /*! This by the way defines the essence of the preconditioning ie:
213  * M_unprec = C^{-1}_L M_prec C^{-1}_R
214  *
215  * This applies whether we use
216  * no spatial preconditioning: M_prec = 1 + C_L D_s C_R
217  * ILU preconditioning M_prec = L + U + L (A-1) U
218  * Here L & U involve both D_s and C_L & C_R, with
219  * A being a left over diagonal piece.
220  */
221  virtual void unprecLinOp(T& chi, const T& psi,
222  enum PlusMinus isign) const
223  {
224 
225  switch (isign) {
226  case PLUS:
227  {
228  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
229 
230  invRightLinOp(tmp1, psi, isign);
231  (*this)(tmp2, tmp1, isign);
232  invLeftLinOp(chi, tmp2, isign);
233  }
234  break;
235  case MINUS:
236  {
237  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
238 
239  invLeftLinOp(tmp1, psi, isign);
240  (*this)(tmp2, tmp1, isign);
241  invRightLinOp(chi, tmp2, isign);
242  }
243  break;
244  default:
245  QDPIO::cerr << "unknown sign" << std::endl;
246  QDP_abort(1);
247  }
248 
249 
250  getFermBC().modifyF(chi);
251  }
252 
253  protected:
254  //! Apply inv (C_L)^{-1}
255  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
256 
257  //! Apply inv (C_R)^{-1}
258  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
259 
260  //! Apply C_L
261  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
262 
263  //! Apply C_R
264  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
265 
266  //! Apply the d/dt of the preconditioned linop
267  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
268 
269 
270  };
271 
272 
273  //-----------------------------------------------------------------------------------
274  //! Time preconditioned linear operator
275  /*! @ingroup linop
276  *
277  * Support for time preconditioned linear operators
278  * Given a matrix M written in block form:
279  *
280  * M = D_t + D_s
281  *
282  * The preconditioning consists of defining matrices C_l and C_r
283  *
284  * so that C_l^{-1} C_r^{-1} = D_t
285  *
286  * and \gamma_5 C_l^{-1} \gamma_5 = C_r^\dagger
287  * and \gamma_5 C_r^{-1} \gamma_5 = C_l^\dagger
288  *
289  * The preconditioned op is then
290  *
291  * M = 1 + C_l D_s C_r
292  *
293  * and
294  *
295  * M^\dag = 1 + C_r^\dag D^\dag_s C_l^\dag
296  *
297  * The simplest way of performing the preconditioning is
298  * To have no other preconditioning than the time preconditioning
299  */
300 
301  //! For now, no forces just yet - come later...
302  template<typename T, typename P, typename Q>
303  class UnprecSpaceCentralPrecTimeLinearOperator : public CentralTimePrecLinearOperator<T,P,Q>
304  {
305  public:
306  //! Virtual destructor to help with cleanup;
307  virtual ~UnprecSpaceCentralPrecTimeLinearOperator() {}
308 
309  //! Defined on the entire lattice
310  const Subset& subset() const {return all;}
311 
312  //! Return the fermion BC object for this linear operator
313  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
314 
315 
316  //! The time direction
317  virtual int tDir() const = 0;
318 
319  //! Apply inv (C_L)^{-1}
320  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
321 
322  //! Apply inv (C_R)^{-1}
323  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
324 
325  //! Apply C_L
326  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
327 
328  //! Apply C_R
329  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
330 
331 
332  //! Apply the the space block onto a source std::vector
333  virtual void spaceLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
334 
335  //! Override operator() for particular preconditioning
336  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const
337  {
338  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
339 
340  switch (isign)
341  {
342  case PLUS:
343  // chi = ( 1 + C_L D_s C_R ) psi
344  cRightLinOp(tmp1, psi, isign);
345  spaceLinOp(tmp2, tmp1, isign);
346  cLeftLinOp(tmp1, tmp2, isign);
347 
348  chi = psi + tmp1;
349 
350 
351  break;
352 
353  case MINUS:
354  // chi = (1 + C_R^\dag D_s C_L^\dag ) psi
355  cLeftLinOp(tmp1, psi, isign);
356  spaceLinOp(tmp2, tmp1, isign);
357  cRightLinOp(tmp1, tmp2, isign);
358 
359  chi = psi + tmp1;
360 
361  break;
362 
363  default:
364  QDPIO::cerr << "unknown sign" << std::endl;
365  QDP_abort(1);
366  }
367 
368  getFermBC().modifyF(chi);
369 
370  }
371 
372 
373 
374  //! Apply the d/dt of the preconditioned linop
375  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
376 
377  //! Get log det ( T^\dag T )
378  virtual Double logDetTDagT(void) const = 0;
379 
380  //! Get the force due to the det T^\dag T bit
381  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const = 0;
382 
383 
384  };
385 
386 
387 
388  //-----------------------------------------------------------------------------------
389  //! Time preconditioned linear operator
390  /*! @ingroup linop
391  *
392  * Support for time preconditioned linear operators with ILU spatial preconditioning
393  * Given a matrix M written in block form:
394  *
395  * M = D_t + D_s + A
396  *
397  * Preconditioning consists of writing:
398  *
399  * D_t = P_{-} T + P_{+} T^\dagger
400  *
401  * Then define T+ = P_{+} + P_{-} T
402  * T- = P_{-} + P_{+} T^\dagger
403  *
404  * T is same as before - essentailly the T+ and T- are like the
405  * C_L^{-1} and C_R^{-1} of the unpreconditioned case.
406  *
407  *
408  * We now write C_L^{-1} = [ T+_{ee} Ds_{eo} ]
409  * [ 0 T+_{oo} ]
410  *
411  * and C_R^{-1} = [ T-_{ee} 0 ]
412  * [ Ds_{oe} T-_{oo} ]
413  *
414  * and
415  * M = C_L^{-1} + C_R^{-1} + tilda{A}
416  *
417  * where tilde{A} = A - 1
418  *
419  * So trivial wilson like case: A = 0, clover case A = -1/4 sigma_munu F_munu
420  *
421  * The preconditioning is: M = C_L^{-1} C_L M C_R C_R^{-1} = C_L^{-1} ( \tilde{M} ) C_R^{-1}
422  * with the preconditioned matrix:
423  *
424  * \tilde{M} = C_R + C_L + C_L (A - 1) C_R
425  *
426  * C_R = [ ( T-_{ee} )^{-1} 0 ]
427  * [-( T-_{oo} )^{-1} Ds_{oe} ( T-_{ee} )^{-1} (T-_{oo})^{-1} ]
428  *
429  * and
430  *
431  * C_L = [ ( T+_{ee} )^{-1} -( T+_{ee} )^{-1} Ds_{eo} ( T+_{oo} )^{-1} ]
432  * [ 0 T+_{oo}^{-1} ]
433  *
434  *
435  * why is this a preconditioning? Certainly it is an ILU form C_R is lower, C_L is upper and there is a residual term... C_L C_R
436  *
437  */
438 
439  //! For now, no forces just yet - come later...
440  template<typename T, typename P, typename Q>
441  class ILUPrecSpaceCentralPrecTimeLinearOperator : public CentralTimePrecLinearOperator<T,P,Q>
442  {
443  public:
444  //! Virtual destructor to help with cleanup;
445  virtual ~ILUPrecSpaceCentralPrecTimeLinearOperator() {}
446 
447  //! Defined on the entire lattice
448  const Subset& subset() const {return all;}
449 
450  //! Return the fermion BC object for this linear operator
451  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
452 
453  //! The time direction
454  virtual int tDir() const = 0;
455 
456  //! Apply inv (C_L)^{-1}
457  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
458 
459  //! Apply inv (C_R)^{-1}
460  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
461 
462  //! Apply C_L
463  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
464 
465  //! Apply C_R
466  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
467 
468  //! Apply A - 1
469  virtual void AMinusOneOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
470 
471  //! Override operator() For this preconditioning
472  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const
473  {
474  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
475 
476  switch (isign)
477  {
478  case PLUS:
479  // Eisenstat's Trick:
480  // chi = ( C_L + C_R + C_L (A - 1 ) C_R ) psi
481  // = C_L ( 1 + (A - 1 ) C_R) psi + C_R psi
482  //
483  // eval as tmp1 = C_R psi
484  // tmp2 = psi + (A-1)tmp1
485  // chi = C_L tmp2
486  // chi += tmp1
487  cRightLinOp(tmp1, psi, isign);
488  AMinusOneOp(tmp2, tmp1,isign);
489  tmp2 +=psi;
490  cLeftLinOp(chi, tmp2, isign);
491  chi += tmp1;
492 
493  break;
494 
495  case MINUS:
496  // Eisenstat's Trick:
497 
498  // chi = ( C_R^\dag + C_L^\dag + C_R^\dag (A-1)^\dag C_L^\dag ) \psi
499  // = C^R\dag ( 1 + (A-1)^\dag C_L^\dag ) psi + C_L^\dag \psi
500 
501  cLeftLinOp(tmp1, psi, isign);
502  AMinusOneOp(tmp2, tmp1, isign);
503  tmp2 += psi;
504  cRightLinOp(chi, tmp2, isign);
505  chi += tmp1;
506 
507  break;
508 
509  default:
510  QDPIO::cerr << "unknown sign" << std::endl;
511  QDP_abort(1);
512  }
513  getFermBC().modifyF(chi);
514 
515  }
516 
517  //! Apply the UNPRECONDITIONED operator onto a source std::vector
518  /*! Mainly intended for debugging */
519  virtual void unprecLinOp(T& chi, const T& psi,
520  enum PlusMinus isign) const
521  {
522 
523  switch (isign) {
524  case PLUS:
525  {
526  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
527 
528  invCRightLinOp(tmp1, psi, isign);
529  (*this)(tmp2, tmp1, isign);
530  invCLeftLinOp(chi, tmp2, isign);
531  }
532  break;
533  case MINUS:
534  {
535  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
536 
537  invCLeftLinOp(tmp1, psi, isign);
538  (*this)(tmp2, tmp1, isign);
539  invCRightLinOp(chi, tmp2, isign);
540  }
541  break;
542  default:
543  QDPIO::cerr << "unknown sign" << std::endl;
544  QDP_abort(1);
545  }
546 
547  getFermBC().modifyF(chi);
548  }
549 
550  virtual void derivCLeft(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
551 
552  virtual void derivCRight(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
553 
554  virtual void derivAMinusOne(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
555 
556  //! Apply the d/dt of the preconditioned linop
557  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const {
558  P ds_tmp;
559  T T_1, T_2, T_3, T_4;
560 
561  switch( isign ) {
562  case PLUS :
563  {
564  // T_1 = C_R Y
565  cRightLinOp(T_1, Y, PLUS);
566 
567  // T_2 = C_L^\dag X => T_2^\dag = X^\dag C_L
568  cLeftLinOp(T_2, X, MINUS);
569 
570 
571  // T_3 = (A-1)^\dagger T_2
572  AMinusOneOp(T_3, T_2, MINUS);
573 
574  // T_3 = (A-1)^\dagger T_2 + X => T_3^\dagger = X^\dagger
575  // + T^2\dagger (A - 1)
576  // = X^\dagger [ 1 + C_L (A-1)]
577  T_3 += X;
578 
579  // T_4 = (A-1) T_1
580  AMinusOneOp(T_4, T_1, PLUS);
581  // T_4 = Y + (A-1) T_1 = Y + (A-1) C_R Y = [ 1 + (A-1)C_R ] Y
582  T_4 += Y;
583 
584 
585  // T_2^\dagger \dot( A-1 ) T_1
586  derivAMinusOne(ds_u, T_2, T_1, PLUS);
587 
588  // T_3^\dagger dot(C_R) Y
589  derivCRight(ds_tmp, T_3, Y, PLUS);
590  ds_u += ds_tmp;
591 
592  derivCLeft(ds_tmp, X, T_4, PLUS);
593  ds_u += ds_tmp;
594  }
595  break;
596  case MINUS :
597  {
598  // T_1 = C_L^\dag Y
599  cLeftLinOp(T_1, Y, MINUS);
600 
601  // T_2 = C_R X
602  cRightLinOp(T_2, X, PLUS);
603 
604  // T_3 = Y + (A-1)^\dag T_1
605  AMinusOneOp(T_3, T_1, MINUS);
606  T_3 += Y;
607 
608  // T_4 = X + (A-1) T_2
609  AMinusOneOp(T_4, T_2, PLUS);
610  T_4 += X;
611 
612  // T_2^\dagger dot(A)^\dagger T_1
613  derivAMinusOne(ds_u, T_2, T_1, MINUS);
614 
615  // T_4^\dagger dot(C_L)^\dagger Y
616  derivCLeft(ds_tmp, T_4, Y, MINUS);
617  ds_u += ds_tmp;
618 
619  // X^\dagger dot(C_R)^\dagger T_3
620  derivCRight(ds_tmp, X, T_3, MINUS);
621  ds_u += ds_tmp;
622 
623  }
624  break;
625  default:
626  QDPIO::cerr << "Bad Case: Should never get here" << std::endl;
627  QDP_abort(1);
628  }
629 
630  getFermBC().zero(ds_u);
631  }
632 
633  //! Get log det ( T^\dag T )
634  virtual Double logDetTDagT(void) const = 0;
635 
636  //! Get the force due to the det T^\dag T bit
637  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const = 0;
638 
639  };
640 
641  //-----------------------------------------------------------------------------------
642  //! Time preconditioned linear operator
643  /*! @ingroup linop
644  *
645  * Support for time preconditioned linear operators with ILU spatial preconditioning
646  * Given a matrix M written in block form:
647  *
648  * M = D_t + D_s
649  *
650  * Preconditioning consists of writing:
651  *
652  * D_t = P_{-} T + P_{+} T^\dagger
653  *
654  * Then define C_L^{-1} = P_{+} + P_{-} T
655  * C_R^{-1} = P_{-} + P_{+} T^\dagger
656  *
657  *
658  * Now we can write
659  *
660  * C_L M C_R = ( 1 + C_L D_s C_R )
661  *
662  *
663  * And now we write D_s in a 3d even odd preconditioned form:
664  *
665  * and M = [ 1 C_L D_s^{eo} C_R ]
666  * [ C_L D^{oe}_s C_R 1 ]
667  *
668  *
669  * and this can be Schur decomposed as
670  * [ 1 0 ] [ 1 0 ] [ 1 C_L D^{oe} C_R ]
671  * [C_L D^{oe} C_R 1 ] [ 0 1 - C_L D^{oe}_s C_R C_L D_s^{eo} C_R ] [ 0 1 ]
672  *
673  */
674 
675 
676  //! For now, no forces just yet - come later...
677  template<typename T, typename P, typename Q>
678  class EO3DPrecSpaceCentralPrecTimeLinearOperator : public CentralTimePrecLinearOperator<T,P,Q>
679  {
680  public:
681  //! Virtual destructor to help with cleanup;
682  virtual ~EO3DPrecSpaceCentralPrecTimeLinearOperator() {}
683 
684  //! Defined on the cb 1
685  const Subset& subset() const {return rb3[1];}
686 
687  //! Return the fermion BC object for this linear operator
688  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
689 
690  //! The time direction
691  virtual int tDir() const = 0;
692 
693 
694 
695 
696  //! Apply inv (C_L)^{-1}
697  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3) const = 0;
698 
699  //! Apply inv (C_R)^{-1}
700  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3) const = 0;
701 
702  //! Apply C_L
703  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3) const = 0;
704 
705  //! Apply C_R
706  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3) const = 0;
707 
708  //! Apply inv (C_L)^{-1}
709  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
710  invCLeftLinOp(chi, psi, isign, 0);
711  invCLeftLinOp(chi, psi, isign, 1);
712  }
713 
714  //! Apply inv (C_R)^{-1}
715  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
716  invCRightLinOp(chi, psi, isign, 0);
717  invCRightLinOp(chi, psi, isign, 1);
718  }
719 
720  //! Apply C_L
721  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
722  cLeftLinOp(chi, psi, isign, 0);
723  cLeftLinOp(chi, psi, isign, 1);
724  }
725 
726  //! Apply C_R
727  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
728  cRightLinOp(chi, psi, isign, 0);
729  cRightLinOp(chi, psi, isign, 1);
730  }
731 
732 
733  //! 3D preconditioning
734 
735  //! M_{ee} where the checkerboarding is in 3d
736  virtual void evenEvenLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
737 
738  //! M_{eo} where the checkerboarding is in 3d
739  virtual void evenOddLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
740 
741  //! M_{oe} where the checkerboarding is in 3d
742  virtual void oddEvenLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
743 
744  //! M_{oo} where the checkerboarding is in 3d
745  virtual void oddOddLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
746 
747  //! M_{ee}^{-1} where the checkerboarding is in 3d - this is awkward for clover
748  // I think. But let's check it for wilson first
749  virtual void evenEvenInvLinOp(T& chi, const T& psi, enum PlusMinus isign) const = 0;
750 
751 
752 
753  //! Apply the operator onto a source std::vector
754  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const
755  {
756  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
757 
758  switch (isign)
759  {
760  case PLUS:
761  // \chi = ( M_oo - M_oe M^{-1}_ee M_eo ) \psi
762 
763  // chi = M_oo \psi
764  oddOddLinOp(chi, psi, PLUS);
765 
766  // tmp2 = M_oe M^{-1}_ee M_eo \psi
767  evenOddLinOp(tmp2, psi, PLUS);
768  evenEvenInvLinOp(tmp1, tmp2, PLUS);
769  oddEvenLinOp(tmp2, tmp1, PLUS);
770 
771  chi[rb3[1]] -= tmp2;
772  break;
773 
774  case MINUS:
775  // \chi = ( M_oo - M_oe M^{-1}_ee M_eo )^\dagger \psi
776  // = M^\dagger_oo \psi - M^\dagger_{oe} ( M^{-1}_ee )^\dagger M^\dagger{eo}
777  //
778  // NB: Daggering acts on checkerboarding to achieve the result above.
779 
780  oddOddLinOp(chi, psi, MINUS);
781 
782  // tmp2 = M_oe M^{-1}_ee M_eo \psi
783  evenOddLinOp(tmp2, psi, MINUS);
784  evenEvenInvLinOp(tmp1, tmp2, MINUS);
785  oddEvenLinOp(tmp2, tmp1, MINUS);
786 
787  chi[rb3[1]] -= tmp2;
788  break;
789 
790  default:
791  QDPIO::cerr << "unknown sign" << std::endl;
792  QDP_abort(1);
793  }
794 
795  getFermBC().modifyF(chi, rb3[1]);
796  }
797 
798  //! Apply the UNPRECONDITIONED operator onto a source std::vector
799  /*! Mainly intended for debugging */
800  virtual void unprecLinOp(T& chi, const T& psi,
801  enum PlusMinus isign) const
802  {
803 
804  T tmp1, tmp2,tmp3;
805  moveToFastMemoryHint(tmp1);
806  moveToFastMemoryHint(tmp2);
807  moveToFastMemoryHint(tmp3);
808 
809  switch (isign) {
810  case PLUS:
811  {
812 
813  // tmp 1 = C_R^{-1} \psi
814  invCRightLinOp(tmp1, psi, isign);
815 
816  // Now apply ( 1 M_ee^{-1} M_eo )
817  // ( 0 1 )
818  evenOddLinOp(tmp2, tmp1, PLUS);
819  evenEvenInvLinOp(tmp3,tmp2, PLUS);
820  tmp1[rb3[0]] += tmp3;
821 
822  // rb[1] part
823  (*this)(tmp2, tmp1, PLUS);
824 
825  // rb[0] part
826  evenEvenLinOp(tmp2, tmp1, PLUS);
827 
828  // Now apply ( 1 0 )
829  // (M_oe M_ee^{-1} 1 )
830  evenEvenInvLinOp(tmp1, tmp2, PLUS);
831  oddEvenLinOp(tmp3, tmp1, PLUS);
832  tmp2[rb3[1]] += tmp3;
833 
834 
835  // chi = C_L^{-1} tmp2
836  invCLeftLinOp(chi, tmp2, isign);
837  }
838  break;
839  case MINUS:
840  {
841 
842  moveToFastMemoryHint(tmp1);
843  moveToFastMemoryHint(tmp2);
844  moveToFastMemoryHint(tmp3);
845 
846 
847  // tmp 1 = C_R^{-\dagger} \psi
848  invCLeftLinOp(tmp1, psi, isign);
849 
850  // Now apply ( 1 M_ee^{-\dagger} M_eo^\dagger )
851  // ( 0 1 )
852  evenOddLinOp(tmp2, tmp1, isign);
853  evenEvenInvLinOp(tmp3,tmp2, isign);
854  tmp1[rb3[0]] += tmp3;
855 
856  // rb[1] part
857  (*this)(tmp2, tmp1, isign);
858 
859  // rb[0] part
860  evenEvenLinOp(tmp2, tmp1, isign);
861 
862  // Now apply ( 1 0 )
863  // (M_oe^\dagger M_ee^{-dagger} 1 )
864 
865  evenEvenInvLinOp(tmp1, tmp2, isign);
866  oddEvenLinOp(tmp3, tmp1, isign);
867  tmp2[rb3[1]] += tmp3;
868 
869 
870  // chi = C_L^{-\dagger} tmp2
871  invCRightLinOp(chi, tmp2, isign);
872 
873  }
874  break;
875  default:
876  QDPIO::cerr << "unknown sign" << std::endl;
877  QDP_abort(1);
878  }
879 
880  getFermBC().modifyF(chi);
881  }
882 
883  //! Apply the d/dt of the preconditioned linop
884  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const = 0;
885 
886  //! Get log det ( T^\dag T )
887  virtual Double logDetTDagT(void) const = 0;
888 
889  //! Get the force due to the det T^\dag T bit
890  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const = 0;
891 
892  };
893 
894 
895 
896 
897 
898 
899  //-----------------------------------------------------------------------------------
900  //! Time preconditioned linear operator
901  /*! @ingroup linop
902  *
903  * Support for time preconditioned linear operators with Mike's style ILU spatial preconditioning
904 
905  */
906 
907  //! For now, no forces just yet - come later...
908  template<typename T, typename P, typename Q>
909  class ILU2PrecSpaceCentralPrecTimeLinearOperator : public Central2TimePrecLinearOperator<T,P,Q>
910  {
911  public:
912  //! Virtual destructor to help with cleanup;
913  virtual ~ILU2PrecSpaceCentralPrecTimeLinearOperator() {}
914 
915  //! Defined on the entire lattice
916  const Subset& subset() const {return all;}
917 
918  //! Return the fermion BC object for this linear operator
919  virtual const FermBC<T,P,Q>& getFermBC() const = 0;
920 
921  //! The time direction
922  virtual int tDir() const = 0;
923 
924 
925  //! Apply S_L
926  virtual void leftLinOp(T& chi, const T& psi, enum PlusMinus isign) const {
927  if( isign == PLUS ) {
928  T tmp1;
929  cLeftLinOp(chi, psi, PLUS,0); // C^e_L psi_e // done
930  cLeftLinOp(chi, psi, PLUS,1); // C^o_L psi_o
931  DBar(tmp1, chi, PLUS, 1); // CB of TARGET: tmp1 = D^{oe} C^e_L
932  chi[rb3[1]] += tmp1; // chi[1] = C^{o}_L psi_o D^{oe}C^e_L psi_e
933 
934  }
935  else {
936  T tmp1;
937  DBar(tmp1, psi, MINUS, 0); // CB of TARGET: tmp1 = D^\dag^eo psi_o
938  tmp1[rb3[0]] += psi; // tmp1 = psi_e + D^\dag^eo psi_o
939  tmp1[rb3[1]] = psi; // tmp1 = psi_o
940  cLeftLinOp(chi, tmp1, MINUS, 0); // chi = C^e_L [psi_e + D^\dag^eo psi_o]
941  cLeftLinOp(chi, tmp1, MINUS, 1); // chi = C^o_L psi_o
942  }
943  }
944 
945  //! Apply S_R
946  virtual void rightLinOp(T& chi, const T& psi, enum PlusMinus isign) const
947  {
948  if ( isign == PLUS ) {
949  T tmp1;
950  DBar(tmp1, psi, PLUS, 0); // CB of TARGET: tmp1 = D^eo psi_o
951  tmp1[rb3[0]] += psi; // tmp1 = psi_e + D psi_o
952  tmp1[rb3[1]] = psi; // tmp1 = psi_o
953  cRightLinOp(chi, tmp1, PLUS, 0); // chi = C^e_R [psi_e + D psi_o]
954  cRightLinOp(chi, tmp1, PLUS, 1); // chi = C^o_R psi_o
955 
956  }
957  else {
958  T tmp1;
959  cRightLinOp(chi, psi, MINUS,0); // (C^e_R)^\dag psi_e // done
960  cRightLinOp(chi, psi, MINUS,1); // (C^o_R)^\dag psi_o
961 
962  DBar(tmp1, chi, MINUS, 1); // CB of TARGET: tmp1 = D^{oe}\dag C^e_R\dag
963  chi[rb3[1]] += tmp1; // chi[1] = C^{o}_L psi_o D^{oe}C^e_L psi_e
964  }
965  }
966 
967 
968  //! Apply S_L^{-1}
969  virtual void invLeftLinOp(T& chi, const T& psi, enum PlusMinus isign) const
970  {
971  T tmp1, tmp2;
972  if(isign == PLUS ) {
973 
974  DBar(tmp1, psi, PLUS, 1); // tmp1= D^{oe} psi_e
975  tmp2[rb3[1]] = psi - tmp1; // tmp2= -D^{oe} psi_e + psi_o
976 
977  invCLeftLinOp(chi, psi, PLUS, 0);
978  invCLeftLinOp(chi, tmp2, PLUS, 1);
979  }
980  else {
981  invCLeftLinOp(chi, psi, MINUS, 0);
982  invCLeftLinOp(chi, psi, MINUS, 1);
983  DBar(tmp1, chi, MINUS, 0); // tmp1 = D^{eo}^\dagger (C^o_L)^{-dagger} psi_o
984  chi[rb3[0]] -= tmp1;
985 
986  }
987  }
988 
989  //! Apply S_R^{-1}
990  virtual void invRightLinOp(T& chi, const T& psi, enum PlusMinus isign) const
991  {
992  T tmp1, tmp2;
993  if( isign == PLUS) {
994  invCRightLinOp(chi, psi, PLUS, 0);
995  invCRightLinOp(chi, psi, PLUS, 1);
996  DBar(tmp1, chi, PLUS, 0); // tmp1 = D^{eo}^\dagger (C^o_L)^{-dagger} psi_o
997  chi[rb3[0]] -= tmp1;
998  }
999  else {
1000  DBar(tmp1, psi, MINUS, 1); // tmp1= D^{oe}^dagger psi_e
1001  tmp2[rb3[1]] = psi - tmp1; // tmp2= -D^{oe}^dagger psi_e + psi_o
1002 
1003  invCRightLinOp(chi, psi, MINUS, 0);
1004  invCRightLinOp(chi, tmp2, MINUS, 1);
1005  }
1006  }
1007 
1008 
1009 
1010  //! Override operator() For this preconditioning
1011  virtual void operator() (T& chi, const T& psi, enum PlusMinus isign) const
1012  {
1013  T tmp1, tmp2,tmp3; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
1014 
1015  switch (isign)
1016  {
1017  case PLUS:
1018 
1019 
1020  chi = psi;
1021  DBar(tmp1, psi, PLUS, 0); // tmp1 = Dbar^{eo} psi_o
1022 
1023  tmp3[rb3[0]] = tmp1 + psi; // tmp1 = psi_e + Dbar^{eo} psi_o
1024  ABar(tmp2, tmp3, PLUS, 0); // tmp2 = ABar [ psi_e + Dbar^{eo} psi_o ]
1025  chi[rb3[0]] += tmp2;
1026 
1027  ABar(tmp3, psi, PLUS, 1); // t3 = \Bar(A)_o psi_o
1028 
1029  chi[rb3[1]] += tmp3;
1030  tmp3[rb3[0]] = tmp2 - tmp1;
1031  DBar(tmp1, tmp3, PLUS, 1);
1032 
1033  chi[rb3[1]] += tmp1;
1034 
1035 
1036  break;
1037 
1038 
1039  case MINUS:
1040  chi = psi;
1041 
1042  DBar(tmp1, psi, MINUS, 0); // tmp1 = D^\dag_eo psi_o
1043 
1044  tmp2[rb3[0]] = psi + tmp1; // tmp2 = psi_e + D^\dag_eo psi_o
1045  ABar(tmp3, tmp2, MINUS, 0); // t3 = A^\dag{ = psi_e + D^\dag_eo psi_o }
1046 
1047 
1048  chi[rb3[0]] += tmp3;
1049 
1050 
1051  tmp2[rb3[0]] = tmp3 - tmp1; // t2= A^\dag [ psi_e + D^\dag psi_o ] - D^\dag psi_p
1052  DBar(tmp3, tmp2, MINUS, 1);
1053  ABar(tmp1, psi, MINUS, 1);
1054  chi[rb3[1]] += tmp3;
1055  chi[rb3[1]] += tmp1;
1056 
1057  break;
1058 
1059  default:
1060  QDPIO::cerr << "unknown sign" << std::endl;
1061  QDP_abort(1);
1062  }
1063  getFermBC().modifyF(chi);
1064 
1065  }
1066 
1067  //! Apply the UNPRECONDITIONED operator onto a source std::vector
1068  /*! Mainly intended for debugging */
1069  virtual void unprecLinOp(T& chi, const T& psi,
1070  enum PlusMinus isign) const
1071  {
1072 
1073  switch (isign) {
1074  case PLUS:
1075  {
1076  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
1077  invRightLinOp(tmp1, psi, isign);
1078  (*this)(tmp2, tmp1, isign);
1079  invLeftLinOp(chi, tmp2, isign);
1080  }
1081  break;
1082  case MINUS:
1083  {
1084  T tmp1, tmp2; moveToFastMemoryHint(tmp1); moveToFastMemoryHint(tmp2);
1085 
1086  invLeftLinOp(tmp1, psi, isign);
1087  (*this)(tmp2, tmp1, isign);
1088  invRightLinOp(chi, tmp2, isign);
1089  }
1090  break;
1091  default:
1092  QDPIO::cerr << "unknown sign" << std::endl;
1093  QDP_abort(1);
1094  }
1095 
1096  getFermBC().modifyF(chi);
1097  }
1098 
1099  protected:
1100  //! Apply inv (C_L)^{-1}
1101  virtual void invCLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1102 
1103  //! Apply inv (C_R)^{-1}
1104  virtual void invCRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1105 
1106  //! Apply C_L
1107  virtual void cLeftLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1108 
1109  //! Apply C_R
1110  virtual void cRightLinOp(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1111 
1112  virtual void DBar(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
1113  {
1114  T tmp1, tmp2;
1115  if (isign == PLUS) {
1116  cRightLinOp(tmp1, psi, PLUS, (1-cb3d));
1117  Dslash3D(tmp2, tmp1, PLUS, cb3d);
1118  cLeftLinOp(chi, tmp2, PLUS, cb3d);
1119  }
1120  else {
1121  cLeftLinOp(tmp1, psi, MINUS, 1-cb3d);
1122  Dslash3D(tmp2, tmp1, MINUS, cb3d);
1123  cRightLinOp(chi, tmp2, MINUS, cb3d);
1124  }
1125  }
1126 
1127 
1128  virtual void ABar(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const
1129  {
1130  T tmp1, tmp2;
1131  if( isign == PLUS ) {
1132  cRightLinOp(tmp1, psi, PLUS, cb3d);
1133  AH(tmp2, tmp1, PLUS, cb3d);
1134  cLeftLinOp(chi, tmp2, PLUS, cb3d);
1135  }
1136  else {
1137  cLeftLinOp(tmp1, psi, MINUS, cb3d);
1138  AH(tmp2, tmp1, MINUS, cb3d);
1139  cRightLinOp(chi, tmp2, MINUS, cb3d);
1140  }
1141  }
1142 
1143  virtual void Dslash3D(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1144  virtual void AH(T& chi, const T& psi, enum PlusMinus isign, int cb3d) const = 0;
1145 
1146  //! Apply the d/dt of the preconditioned linop
1147  virtual void deriv(P& ds_u, const T& X, const T& Y, enum PlusMinus isign) const
1148  {
1149  QDPIO::cerr << "Not Yet Implemented " << std::endl;
1150  QDP_abort(1);
1151  }
1152 
1153 
1154  //! Get log det ( T^\dag T )
1155  virtual Double logDetTDagT(void) const
1156  {
1157  QDPIO::cerr << "Not Yet Implemented " << std::endl;
1158  QDP_abort(1);
1159  Double tmp; // Make compiler happy
1160  return tmp;
1161  }
1162 
1163  //! Get the force due to the det T^\dag T bit
1164  virtual void derivLogDetTDagT(P& ds_u, enum PlusMinus isign) const
1165  {
1166  QDPIO::cerr<< "Not Yet Implemented" << std::endl;
1167  QDP_abort(1);
1168  }
1169 
1170  };
1171 
1172 }
1173 
1174 #endif
1175 #endif
1176 #endif
1177 
1178 #endif
linearop.h
Linear Operators.
tmp2
Double tmp2
Definition: mesq.cc:30
tmp3
Double tmp3
Definition: mesq.cc:31
Chroma::SimpleBaryonSeqSourceEnv::operator()
multi1d< Hadron2PtContraction_t > operator()(const multi1d< LatticeColorMatrix > &u)
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::tDir
static int tDir()
Definition: coulgauge.cc:14
Chroma::tmp
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
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
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
schroedinger_fermbc_w.h
Fermion action boundary conditions.
T
LatticeFermion T
Definition: t_clover.cc:11
P
multi1d< LatticeColorMatrix > P
Definition: t_clover.cc:13
deriv
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.
Definition: t_precact_4d.cc:141
Generated by doxygen 1.9.1