CHROMA
eoprec_ovdwf_linop_array_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief 4D Even Odd preconditioned Overlap-DWF (Borici) linear operator
3  */
4 
5 #include "chromabase.h"
6 #include "actions/ferm/linop/eoprec_ovdwf_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 m_q_ quark mass (Read)
19  * \param N5_ extent of 5D (Read)
20  */
21  void
22  EvenOddPrecOvDWLinOpArray::create(Handle< FermState<T,P,Q> > fs,
23  const Real& WilsonMass_, const Real& m_q_, int N5_)
24  {
25  WilsonMass = WilsonMass_;
26  m_q = m_q_;
27  a5 = 1.0;
28  N5 = N5_;
29 
30 
31  // Extract a 4D fermstate to use for the Wilson dslash
32  D.create(fs);
33 
34  c5InvTwoKappa = 1.0 - (Nd-WilsonMass) ;
35  c5TwoKappa = 1.0 / c5InvTwoKappa ;
36 
37  b5InvTwoKappa = 1.0 + (Nd-WilsonMass) ;
38  b5TwoKappa = 1.0 / b5InvTwoKappa ;
39 
40  TwoKappa = c5InvTwoKappa/b5InvTwoKappa ;
41  Kappa = TwoKappa/2.0 ;
42 
43  invDfactor =1.0/(1.0 + m_q*pow(TwoKappa,N5)) ;
44 
45  }
46 
47 
48  //! Apply the even-even (odd-odd) coupling piece of the Borici operator
49  /*!
50  * \ingroup linop
51  *
52  * The operator acts on the entire lattice
53  *
54  * \param psi Pseudofermion field (Read)
55  * \param isign Flag ( PLUS | MINUS ) (Read)
56  * \param cb checkerboard ( 0 | 1 ) (Read)
57  */
58  void
59  EvenOddPrecOvDWLinOpArray::applyDiag(multi1d<LatticeFermion>& chi,
60  const multi1d<LatticeFermion>& psi,
61  enum PlusMinus isign,
62  const int cb) const
63  {
64  START_CODE();
65 
66  if( chi.size() != N5 ) chi.resize(N5);
67 
68  switch ( isign ) {
69 
70  case PLUS:
71  {
72  for(int s(1);s<N5-1;s++) // 1/2k psi[s] + P_- * psi[s+1] + P_+ * psi[s-1]
73  chi[s][rb[cb]] = b5InvTwoKappa*psi[s] -
74  (0.5*c5InvTwoKappa)*( psi[s+1] + psi[s-1] + GammaConst<Ns,Ns*Ns-1>()*(psi[s-1] - psi[s+1]) ) ;
75 
76  int N5m1(N5-1) ;
77  //s=0 -- 1/2k psi[0] - P_- * psi[1] + mf* P_+ * psi[N5-1]
78  chi[0][rb[cb]] = b5InvTwoKappa*psi[0] -
79  (0.5*c5InvTwoKappa)*( psi[1] - m_q*psi[N5m1] - GammaConst<Ns,Ns*Ns-1>()*(m_q*psi[N5m1] + psi[1]) ) ;
80 
81  int N5m2(N5-2);
82  //s=N5-1 -- 1/2k psi[N5-1] +mf* P_- * psi[0] - P_+ * psi[N5-2]
83  chi[N5m1][rb[cb]] = b5InvTwoKappa*psi[N5m1] -
84  (0.5*c5InvTwoKappa)*( psi[N5m2] - m_q *psi[0] + GammaConst<Ns,Ns*Ns-1>()*(psi[N5m2] + m_q * psi[0]) );
85  }
86  break ;
87 
88  case MINUS:
89  {
90  for(int s(1);s<N5-1;s++) // 1/2k psi[s] - P_+ * psi[s+1] - P_- * psi[s-1]
91  chi[s][rb[cb]] = b5InvTwoKappa*psi[s] -
92  (0.5*c5InvTwoKappa)*( psi[s+1] + psi[s-1] + GammaConst<Ns,Ns*Ns-1>()*(psi[s+1] - psi[s-1]) ) ;
93 
94  int N5m1(N5-1) ;
95  //s=0 -- 1/2k psi[0] - P_+ * psi[1] + mf* P_- * psi[N5-1]
96  chi[0][rb[cb]] = b5InvTwoKappa*psi[0] -
97  (0.5*c5InvTwoKappa)*( psi[1] - m_q*psi[N5m1] + GammaConst<Ns,Ns*Ns-1>()*( psi[1]+m_q*psi[N5m1]) ) ;
98 
99  int N5m2(N5-2);
100  //s=N5-1 -- 1/2k psi[N5-1] + mf* P_+ * psi[0] - P_- * psi[N5-2]
101  chi[N5m1][rb[cb]] = b5InvTwoKappa*psi[N5m1] -
102  (0.5*c5InvTwoKappa)*( psi[N5m2] - m_q *psi[0] - GammaConst<Ns,Ns*Ns-1>()*(psi[N5m2] + m_q * psi[0]) );
103  }
104  break ;
105  }
106 
107  END_CODE();
108  }
109 
110 
111  //! Apply the inverse even-even (odd-odd) coupling piece of the Borici fermion operator
112  /*!
113  * \ingroup linop
114  *
115  * The operator acts on the entire lattice
116  *
117  * \param psi Pseudofermion field (Read)
118  * \param isign Flag ( PLUS | MINUS ) (Read)
119  * \param cb checkerboard ( 0 | 1 ) (Read)
120  */
121  void
122  EvenOddPrecOvDWLinOpArray::applyDiagInv(multi1d<LatticeFermion>& chi,
123  const multi1d<LatticeFermion>& psi,
124  enum PlusMinus isign,
125  const int cb) const
126  {
127  START_CODE();
128 
129  if( chi.size() != N5 ) chi.resize(N5);
130 
131  // Copy and scale by TwoKappa (1/M0)
132  for(int s(0);s<N5;s++)
133  chi[s][rb[cb]] = b5TwoKappa * psi[s] ;
134 
135 
136  switch ( isign ) {
137 
138  case PLUS:
139  {
140 
141  // First apply the inverse of Lm
142  Real fact(0.5*m_q*TwoKappa) ;
143  for(int s(0);s<N5-1;s++){
144  chi[N5-1][rb[cb]] -= fact * (chi[s] - GammaConst<Ns,Ns*Ns-1>()*chi[s]) ;
145  fact *= TwoKappa ;
146  }
147 
148  //Now apply the inverse of L. Forward elimination
149  for(int s(1);s<N5;s++)
150  chi[s][rb[cb]] += Kappa*(chi[s-1] + GammaConst<Ns,Ns*Ns-1>()*chi[s-1]) ;
151 
152  //The inverse of D now
153  chi[N5-1][rb[cb]] *= invDfactor ;
154  // That was easy....
155 
156  //The inverse of R. Back substitution...... Getting there!
157  for(int s(N5-2);s>-1;s--)
158  chi[s][rb[cb]] += Kappa*(chi[s+1] - GammaConst<Ns,Ns*Ns-1>()*chi[s+1]) ;
159 
160  //Finally the inverse of Rm
161  LatticeFermion tt; moveToFastMemoryHint(tt);
162  fact = 0.5*m_q*TwoKappa;
163  tt[rb[cb]] = fact*(chi[N5-1] + GammaConst<Ns,Ns*Ns-1>()*chi[N5-1]);
164  for(int s(0);s<N5-1;s++){
165  chi[s][rb[cb]] -= tt ;
166  tt[rb[cb]] *= TwoKappa ;
167  }
168  }
169  break ;
170 
171  case MINUS:
172  {
173 
174  // First apply the inverse of Lm
175  Real fact(0.5*m_q*TwoKappa) ;
176  for(int s(0);s<N5-1;s++){
177  chi[N5-1][rb[cb]] -= fact * (chi[s] + GammaConst<Ns,Ns*Ns-1>()*chi[s]) ;
178  fact *= TwoKappa ;
179  }
180 
181  //Now apply the inverse of L. Forward elimination
182  for(int s(1);s<N5;s++)
183  chi[s][rb[cb]] += Kappa*(chi[s-1] - GammaConst<Ns,Ns*Ns-1>()*chi[s-1]) ;
184 
185  //The inverse of D now
186  chi[N5-1][rb[cb]] *= invDfactor ;
187  // That was easy....
188 
189  //The inverse of R. Back substitution...... Getting there!
190  for(int s(N5-2);s>-1;s--)
191  chi[s][rb[cb]] += Kappa*(chi[s+1] + GammaConst<Ns,Ns*Ns-1>()*chi[s+1]) ;
192 
193  //Finally the inverse of Rm
194  LatticeFermion tt; moveToFastMemoryHint(tt);
195  tt[rb[cb]] = (0.5*m_q*TwoKappa)*(chi[N5-1] - GammaConst<Ns,Ns*Ns-1>()*chi[N5-1]);
196  for(int s(0);s<N5-1;s++){
197  chi[s][rb[cb]] -= tt ;
198  tt[rb[cb]] *= TwoKappa ;
199  }
200  }
201  break ;
202  }
203 
204  END_CODE();
205  }
206 
207 
208  //! Apply the even-odd (odd-even) coupling piece of the Borici operator
209  /*!
210  * \ingroup linop
211  *
212  * The operator acts on the entire lattice
213  *
214  * \param psi Pseudofermion field (Read)
215  * \param isign Flag ( PLUS | MINUS ) (Read)
216  * \param cb checkerboard ( 0 | 1 ) (Read)
217  */
218  void
219  EvenOddPrecOvDWLinOpArray::applyOffDiag(multi1d<LatticeFermion>& chi,
220  const multi1d<LatticeFermion>& psi,
221  enum PlusMinus isign,
222  const int cb) const
223  {
224  START_CODE();
225 
226  if( chi.size() != N5 ) chi.resize(N5);
227 
228 
229  switch ( isign )
230  {
231  case PLUS:
232  {
233  LatticeFermion tmp; moveToFastMemoryHint(tmp);
234  int otherCB = (cb + 1)%2 ;
235  for(int s(1);s<N5-1;s++)
236  {
237  tmp[rb[otherCB]] = psi[s] +
238  0.5*( psi[s+1] + psi[s-1] +
239  GammaConst<Ns,Ns*Ns-1>()*(psi[s-1] - psi[s+1]));
240  D.apply(chi[s],tmp,isign,cb);
241  chi[s][rb[cb]] *= (-0.5);
242  }
243 
244  int N5m1(N5-1);
245  tmp[rb[otherCB]] = psi[0] +
246  0.5*(psi[1] - m_q*psi[N5m1]
247  - GammaConst<Ns,Ns*Ns-1>()*(m_q*psi[N5m1] + psi[1]));
248  D.apply(chi[0],tmp,isign,cb);
249  chi[0][rb[cb]] *= (-0.5);
250 
251  int N5m2(N5-2);
252  tmp[rb[otherCB]] = psi[N5m1] +
253  0.5*(psi[N5m2] - m_q *psi[0]
254  + GammaConst<Ns,Ns*Ns-1>()*(psi[N5m2] + m_q * psi[0]));
255 
256  D.apply(chi[N5m1],tmp,isign,cb);
257  chi[N5m1][rb[cb]] *= (-0.5);
258 
259  }
260  break ;
261 
262  case MINUS:
263  {
264  multi1d<LatticeFermion> tmp(N5) ; moveToFastMemoryHint(tmp);
265  for(int s(0);s<N5;s++){
266  D.apply(tmp[s],psi[s],isign,cb);
267  tmp[s][rb[cb]] *= (-0.5) ;
268  }
269  for(int s(1);s<N5-1;s++){
270  chi[s][rb[cb]] = tmp[s] +
271  0.5*(tmp[s+1] + tmp[s-1] -
272  GammaConst<Ns,Ns*Ns-1>()*(tmp[s-1]-tmp[s+1]));
273  }
274  int N5m1(N5-1);
275  chi[0][rb[cb]] = tmp[0] +
276  0.5*(tmp[1] - m_q*tmp[N5m1] +
277  GammaConst<Ns,Ns*Ns-1>()*(m_q*tmp[N5m1] + tmp[1]));
278 
279 
280  int N5m2(N5-2);
281  chi[N5m1][rb[cb]] = tmp[N5m1] +
282  0.5*(tmp[N5m2] - m_q *tmp[0] -
283  GammaConst<Ns,Ns*Ns-1>()*(tmp[N5m2] + m_q * tmp[0]));
284 
285  }
286  break ;
287  }
288 
289  //Done! That was not that bad after all....
290  //See, I told you so...
291  END_CODE();
292  }
293 
294  //! Apply the Dminus operator on a lattice fermion. See my notes ;-)
295  void
296  EvenOddPrecOvDWLinOpArray::Dminus(LatticeFermion& chi,
297  const LatticeFermion& psi,
298  enum PlusMinus isign,
299  int s5) const
300  {
301  LatticeFermion tt ; moveToFastMemoryHint(tt);
302  D.apply(tt,psi,isign,0);
303  D.apply(tt,psi,isign,1);
304  chi = c5InvTwoKappa*psi +0.5*tt ; //really -(-.5)D
305  }
306 
307 } // End Namespace Chroma
308 
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
a5
Real a5
Definition: eoprec_dwf_fermact_array_w.cc:79
N5
int N5
Definition: eoprec_dwf_fermact_array_w.cc:78
eoprec_ovdwf_linop_array_w.h
4D Even Odd preconditioned Overlap-DWF (Borici) linear operator
s
unsigned s
Definition: ldumul_w.cc:37
tmp
Double tmp
Definition: meslate.cc:60
Nd
Nd
Definition: meslate.cc:74
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::cb
int cb
Definition: invbicg.cc:120
chi
chi
Definition: pade_trln_w.cc:24
psi
psi
Definition: pade_trln_w.cc:191
Generated by doxygen 1.9.1