CHROMA
lwldslash_qdpopt_w.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! \file
3  * \brief Wilson Dslash linear operator
4  */
5 
6 #ifndef __lwldslash_qdpopt_h__
7 #define __lwldslash_qdpopt_h__
8 
9 #include "state.h"
10 #include "io/aniso_io.h"
12 
13 
14 namespace Chroma
15 {
16  //! General Wilson-Dirac dslash
17  /*!
18  * \ingroup linop
19  *
20  * DSLASH
21  *
22  * This routine is specific to Wilson fermions!
23  *
24  * Description:
25  *
26  * This routine applies the operator D' to Psi, putting the result in Chi.
27  *
28  * Nd-1
29  * ---
30  * \
31  * chi(x) := > U (x) (1 - isign gamma ) psi(x+mu)
32  * / mu mu
33  * ---
34  * mu=0
35  *
36  * Nd-1
37  * ---
38  * \ +
39  * + > U (x-mu) (1 + isign gamma ) psi(x-mu)
40  * / mu mu
41  * ---
42  * mu=0
43  *
44  */
45 
46  template<typename T, typename P, typename Q>
47  class QDPWilsonDslashOptT : public WilsonDslashBase<T,P,Q>
48  {
49  public:
50 
51  //! Empty constructor. Must use create later
53 
54  //! Full constructor
56 
57  //! Full constructor with anisotropy
59  const AnisoParam_t& aniso_);
60 
61  //! Full constructor with general coefficients
63  const multi1d<Real>& coeffs_);
64 
65  //! Creation routine
67 
68  //! Creation routine with anisotropy
70  const AnisoParam_t& aniso_);
71 
72  //! Full constructor with general coefficients
74  const multi1d<Real>& coeffs_);
75 
76  //! No real need for cleanup here
78 
79  /**
80  * Apply a dslash
81  *
82  * \param chi result (Write)
83  * \param psi source (Read)
84  * \param isign D'^dag or D' ( MINUS | PLUS ) resp. (Read)
85  * \param cb Checkerboard of OUTPUT std::vector (Read)
86  *
87  * \return The output of applying dslash on psi
88  */
89  void apply (T& chi, const T& psi, enum PlusMinus isign, int cb) const;
90 
91  //! Return the fermion BC object for this linear operator
92  const FermBC<T,P,Q>& getFermBC() const {return *fbc;}
93 
94  protected:
95  //! Get the anisotropy parameters
96  const multi1d<Real>& getCoeffs() const {return coeffs;}
97 
98  private:
99  multi1d<Real> coeffs; /*!< Nd array of coefficients of terms in the action */
101  Q u;
102  };
103 
104 
105 #if 0
106  template<typename T>
107  class HalfFermionType{};
108 
109  template<>
110  class HalfFermionType<LatticeFermionF> {
111  public:
112  typedef LatticeHalfFermionF Type_t;
113  };
114 
115  template<>
116  class HalfFermionType<LatticeFermionD> {
117  public:
118  typedef LatticeHalfFermionD Type_t;
119  };
120 #endif
121 
122  //! General Wilson-Dirac dslash
123  /*! \ingroup linop
124  * DSLASH
125  *
126  * This routine is specific to Wilson fermions!
127  *
128  * Description:
129  *
130  * This routine applies the operator D' to Psi, putting the result in Chi.
131  *
132  * Nd-1
133  * ---
134  * \
135  * chi(x) := > U (x) (1 - isign gamma ) psi(x+mu)
136  * / mu mu
137  * ---
138  * mu=0
139  *
140  * Nd-1
141  * ---
142  * \ +
143  * + > U (x-mu) (1 + isign gamma ) psi(x-mu)
144  * / mu mu
145  * ---
146  * mu=0
147  *
148  */
149 
150  //! Empty constructor
151  template<typename T, typename P, typename Q>
153 
154  //! Full constructor
155  template<typename T, typename P, typename Q>
157  {
158  create(state);
159  }
160 
161  //! Full constructor with anisotropy
162  template<typename T, typename P, typename Q>
164  const AnisoParam_t& aniso_)
165  {
166  create(state, aniso_);
167  }
168 
169  //! Full constructor with general coefficients
170  template<typename T, typename P, typename Q>
172  const multi1d<Real>& coeffs_)
173  {
174  create(state, coeffs_);
175  }
176 
177  //! Creation routine
178  template<typename T, typename P, typename Q>
180  {
181  multi1d<Real> cf(Nd);
182  cf = 1.0;
183  create(state, cf);
184  }
185 
186  //! Creation routine with anisotropy
187  template<typename T, typename P, typename Q>
189  const AnisoParam_t& anisoParam)
190  {
191  START_CODE();
192 
193  create(state, makeFermCoeffs(anisoParam));
194 
195  END_CODE();
196  }
197 
198  //! Full constructor with general coefficients
199  template<typename T, typename P, typename Q>
201  const multi1d<Real>& coeffs_)
202  {
203  START_CODE();
204 
205  // Save a copy of the aniso params original fields and with aniso folded in
206  coeffs = coeffs_;
207 
208  // Save a copy of the fermbc
209  fbc = state->getFermBC();
210 
211  // Sanity check
212  if (fbc.operator->() == 0)
213  {
214  QDPIO::cerr << "QDPWilsonDslashOpt: error: fbc is null" << std::endl;
215  QDP_abort(1);
216  }
217 
218  // Fold in anisotropy
219  u.resize(Nd);
220  for(int mu=0; mu < u.size(); ++mu) {
221  u[mu] = (state->getLinks())[mu];
222  }
223  // Rescale the u fields by the anisotropy
224  for(int mu=0; mu < u.size(); ++mu)
225  {
226  u[mu] *= coeffs[mu];
227  }
228  }
229 
230 
231  //! General Wilson-Dirac dslash
232  /*! \ingroup linop
233  * Wilson dslash
234  *
235  * Arguments:
236  *
237  * \param chi Result (Write)
238  * \param psi Pseudofermion field (Read)
239  * \param isign D'^dag or D' ( MINUS | PLUS ) resp. (Read)
240  * \param cb Checkerboard of OUTPUT std::vector (Read)
241  */
242  template<typename T, typename P, typename Q>
243  void
245  enum PlusMinus isign, int cb) const
246  {
247  START_CODE();
248 #if QDP_NC == 3
249  /* F
250  * a2 (x) := U (x) (1 - isign gamma ) psi(x)
251  * mu mu mu
252  */
253  /* B +
254  * a2 (x) := U (x-mu) (1 + isign gamma ) psi(x-mu)
255  * mu mu mu
256  */
257  // Recontruct the bottom two spinor components from the top two
258  /* F B
259  * chi(x) := sum_mu a2 (x) + a2 (x)
260  * mu mu
261  */
262 
263  /* Why are these lines split? An array syntax would help, but the problem is deeper.
264  * The expression templates require NO variable args (like int's) to a function
265  * and all args must be known at compile time. Hence, the function names carry
266  * (as functions usually do) the meaning (and implicit args) to a function.
267  */
268 
269  int otherCB = (cb == 0 ? 1 : 0);
270 
271  switch (isign)
272  {
273  case PLUS:
274  {
275 
276 
279 
280  tmp[rb[otherCB]] = spinProjectDir0Minus(psi);
281  tmp2[rb[cb]] = shift(tmp, FORWARD, 0);
282  chi[rb[cb]] = spinReconstructDir0Minus(u[0]*tmp2);
283 
284 
285  tmp[rb[otherCB]] = spinProjectDir1Minus(psi);
286  tmp2[rb[cb]] = shift(tmp, FORWARD, 1);
287  chi[rb[cb]] += spinReconstructDir1Minus(u[1]*tmp2);
288 
289  tmp[rb[otherCB]] = spinProjectDir2Minus(psi);
290  tmp2[rb[cb]] = shift(tmp, FORWARD, 2);
291  chi[rb[cb]] += spinReconstructDir2Minus(u[2]*tmp2);
292 
293  tmp[rb[otherCB]] = spinProjectDir3Minus(psi);
294  tmp2[rb[cb]] = shift(tmp, FORWARD, 3);
295  chi[rb[cb]] += spinReconstructDir3Minus(u[3]*tmp2);
296 
297 
298  tmp[rb[otherCB]] = adj(u[0])*spinProjectDir0Plus(psi);
299  tmp2[rb[cb]] = shift(tmp, BACKWARD, 0);
300  chi[rb[cb]] += spinReconstructDir0Plus(tmp2);
301 
302  tmp[rb[otherCB]] = adj(u[1])*spinProjectDir1Plus(psi);
303  tmp2[rb[cb]] = shift(tmp, BACKWARD, 1);
304  chi[rb[cb]] += spinReconstructDir1Plus(tmp2);
305 
306  tmp[rb[otherCB]] = adj(u[2])*spinProjectDir2Plus(psi);
307  tmp2[rb[cb]] = shift(tmp, BACKWARD, 2);
308  chi[rb[cb]] += spinReconstructDir2Plus(tmp2);
309 
310  tmp[rb[otherCB]] = adj(u[3])*spinProjectDir3Plus(psi);
311  tmp2[rb[cb]] = shift(tmp, BACKWARD, 3);
312  chi[rb[cb]] += spinReconstructDir3Plus(tmp2);
313 
314  }
315 
316 
317  break;
318 
319  case MINUS:
320  {
321 
324 
325  tmp[rb[otherCB]] = spinProjectDir0Plus(psi);
326  tmp2[rb[cb]] = shift(tmp, FORWARD, 0);
327  chi[rb[cb]] = spinReconstructDir0Plus(u[0]*tmp2);
328 
329  tmp[rb[otherCB]] = spinProjectDir1Plus(psi);
330  tmp2[rb[cb]] = shift(tmp, FORWARD, 1);
331  chi[rb[cb]] += spinReconstructDir1Plus(u[1]*tmp2);
332 
333  tmp[rb[otherCB]] = spinProjectDir2Plus(psi);
334  tmp2[rb[cb]] = shift(tmp, FORWARD, 2);
335  chi[rb[cb]] += spinReconstructDir2Plus(u[2]*tmp2);
336 
337  tmp[rb[otherCB]] = spinProjectDir3Plus(psi);
338  tmp2[rb[cb]] = shift(tmp, FORWARD, 3);
339  chi[rb[cb]] += spinReconstructDir3Plus(u[3]*tmp2);
340 
341 
342  tmp[rb[otherCB]] = adj(u[0])*spinProjectDir0Minus(psi);
343  tmp2[rb[cb]] = shift(tmp, BACKWARD, 0);
344  chi[rb[cb]] += spinReconstructDir0Minus(tmp2);
345 
346  tmp[rb[otherCB]] = adj(u[1])*spinProjectDir1Minus(psi);
347  tmp2[rb[cb]] = shift(tmp, BACKWARD, 1);
348  chi[rb[cb]] += spinReconstructDir1Minus(tmp2);
349 
350  tmp[rb[otherCB]] = adj(u[2])*spinProjectDir2Minus(psi);
351  tmp2[rb[cb]] = shift(tmp, BACKWARD, 2);
352  chi[rb[cb]] += spinReconstructDir2Minus(tmp2);
353 
354  tmp[rb[otherCB]] = adj(u[3])*spinProjectDir3Minus(psi);
355  tmp2 = shift(tmp, BACKWARD, 3);
356  chi[rb[cb]] += spinReconstructDir3Minus(tmp2);
357 
358  }
359 
360  break;
361  }
362 
363  getFermBC().modifyF(chi, QDP::rb[cb]);
364 #else
365  QDPIO::cerr<<"lwldshash_qdpopt_w: not implemented for NC!=3\n";
366  QDP_abort(13);
367 #endif
368 
369  END_CODE();
370  }
371 
372  typedef QDPWilsonDslashOptT< LatticeFermion,
373  multi1d<LatticeColorMatrix>,
374  multi1d<LatticeColorMatrix> > QDPWilsonDslashOpt;
375 
376  typedef QDPWilsonDslashOptT< LatticeFermionF,
377  multi1d<LatticeColorMatrixF>,
378  multi1d<LatticeColorMatrixF> > QDPWilsonDslashOptF;
379 
380  typedef QDPWilsonDslashOptT< LatticeFermionD,
381  multi1d<LatticeColorMatrixD>,
382  multi1d<LatticeColorMatrixD> > QDPWilsonDslashOptD;
383 
384 
385 } // End Namespace Chroma
386 
387 
388 #endif
Anisotropy parameters.
Base class for all fermion action boundary conditions.
Definition: fermbc.h:20
Support class for fermion actions and linear operators.
Definition: state.h:94
Class for counted reference semantics.
Definition: handle.h:33
General Wilson-Dirac dslash.
Handle< FermBC< T, P, Q > > fbc
const multi1d< Real > & getCoeffs() const
Get the anisotropy parameters.
~QDPWilsonDslashOptT()
No real need for cleanup here.
void create(Handle< FermState< T, P, Q > > state)
Creation routine.
const FermBC< T, P, Q > & getFermBC() const
Return the fermion BC object for this linear operator.
General Wilson-Dirac dslash.
int mu
Definition: cool.cc:24
void apply(T &chi, const T &psi, enum PlusMinus isign, int cb) const
General Wilson-Dirac dslash.
QDPWilsonDslashOptT()
Empty constructor. Must use create later.
Wilson Dslash linear operator.
Nd
Definition: meslate.cc:74
Double tmp2
Definition: mesq.cc:30
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
QDPWilsonDslashOptT< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > QDPWilsonDslashOpt
LinOpSysSolverMGProtoClover::Q Q
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
LinOpSysSolverMGProtoClover::T T
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
QDPWilsonDslashOptT< LatticeFermionF, multi1d< LatticeColorMatrixF >, multi1d< LatticeColorMatrixF > > QDPWilsonDslashOptF
LatticeFermion psi
Definition: mespbg5p_w.cc:35
QDPWilsonDslashOptT< LatticeFermionD, multi1d< LatticeColorMatrixD >, multi1d< LatticeColorMatrixD > > QDPWilsonDslashOptD
START_CODE()
int cb
Definition: invbicg.cc:120
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
multi1d< Real > makeFermCoeffs(const AnisoParam_t &aniso)
Make fermion coefficients.
Definition: aniso_io.cc:63
#define FORWARD
Definition: primitives.h:82
#define BACKWARD
Definition: primitives.h:83
Support class for fermion actions and linear operators.
Parameters for anisotropy.
Definition: aniso_io.h:24