CHROMA
lwldslash_array_qdpopt_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Wilson Dslash linear operator array
3  */
4 
5 #include "chromabase.h"
7 
8 
9 namespace Chroma
10 {
11  //! General Wilson-Dirac dslash of arrays
12  /*! \ingroup linop
13  * DSLASH
14  *
15  * This routine is specific to Wilson fermions!
16  *
17  * Description:
18  *
19  * This routine applies the operator D' to Psi, putting the result in Chi.
20  *
21  * Nd-1
22  * ---
23  * \
24  * chi(x) := > U (x) (1 - isign gamma ) psi(x+mu)
25  * / mu mu
26  * ---
27  * mu=0
28  *
29  * Nd-1
30  * ---
31  * \ +
32  * + > U (x-mu) (1 + isign gamma ) psi(x-mu)
33  * / mu mu
34  * ---
35  * mu=0
36  *
37  */
38 
39 
40  //! Creation routine
42  {
43  multi1d<Real> cf(Nd);
44  cf = 1.0;
45  create(state, N5_, cf);
46  }
47 
48 
49  //! Creation routine with anisotropy
51  const AnisoParam_t& anisoParam)
52  {
53  START_CODE();
54 
55  create(state, N5_, makeFermCoeffs(anisoParam));
56 
57  END_CODE();
58  }
59 
60  //! Creation routine
62  const multi1d<Real>& coeffs_)
63  {
64  START_CODE();
65 
66  N5 = N5_;
67  coeffs = coeffs_;
68 
69  // Save a copy of the fermbc
70  fbc = state->getFermBC();
71 
72  // Sanity check
73  if (fbc.operator->() == 0)
74  {
75  QDPIO::cerr << "QDPWilsonDslashArrayOpt: error: fbc is null" << std::endl;
76  QDP_abort(1);
77  }
78 
79  // Get links
80  u = state->getLinks();
81 
82  // Rescale the u fields by the anisotropy
83  for(int mu=0; mu < u.size(); ++mu)
84  {
85  u[mu] *= coeffs[mu];
86  }
87 
88  END_CODE();
89  }
90 
91 
92  //! General Wilson-Dirac dslash
93  /*! \ingroup linop
94  * Wilson dslash
95  *
96  * Arguments:
97  *
98  * \param chi Result (Write)
99  * \param psi Pseudofermion field (Read)
100  * \param isign D'^dag or D' ( MINUS | PLUS ) resp. (Read)
101  * \param cb Checkerboard of OUTPUT std::vector (Read)
102  */
103  void
104  QDPWilsonDslashArrayOpt::apply (multi1d<LatticeFermion>& chi,
105  const multi1d<LatticeFermion>& psi,
106  enum PlusMinus isign, int cb) const
107  {
108  START_CODE();
109 
110  if( chi.size() != N5 ) chi.resize(N5);
111 
112  for(int n=0; n < N5; ++n)
113  apply(chi[n], psi[n], isign, cb);
114 
115  END_CODE();
116  }
117 
118 
119 
120  //! General Wilson-Dirac dslash
121  /*! \ingroup linop
122  * Wilson dslash
123  *
124  * Arguments:
125  *
126  * \param chi Result (Write)
127  * \param psi Pseudofermion field (Read)
128  * \param isign D'^dag or D' ( MINUS | PLUS ) resp. (Read)
129  * \param cb Checkerboard of OUTPUT std::vector (Read)
130  */
131  void
132  QDPWilsonDslashArrayOpt::apply (LatticeFermion& chi, const LatticeFermion& psi,
133  enum PlusMinus isign, int cb) const
134  {
135  START_CODE();
136 #if (QDP_NC == 2) || (QDP_NC == 3)
137  /* F
138  * a2 (x) := U (x) (1 - isign gamma ) psi(x)
139  * mu mu mu
140  */
141  /* B +
142  * a2 (x) := U (x-mu) (1 + isign gamma ) psi(x-mu)
143  * mu mu mu
144  */
145  // Recontruct the bottom two spinor components from the top two
146  /* F B
147  * chi(x) := sum_mu a2 (x) + a2 (x)
148  * mu mu
149  */
150 
151  /* Why are these lines split? An array syntax would help, but the problem is deeper.
152  * The expression templates require NO variable args (like int's) to a function
153  * and all args must be known at compile time. Hence, the function names carry
154  * (as functions usually do) the meaning (and implicit args) to a function.
155  */
156  int otherCB = (cb == 0 ? 1 : 0);
157 
158  switch (isign)
159  {
160  case PLUS:
161  {
162 
163 
164  LatticeHalfFermion tmp;
165  LatticeHalfFermion tmp2;
166 
167  tmp[rb[otherCB]] = spinProjectDir0Minus(psi);
168  tmp2[rb[cb]] = shift(tmp, FORWARD, 0);
169  chi[rb[cb]] = spinReconstructDir0Minus(u[0]*tmp2);
170 
171 
172  tmp[rb[otherCB]] = spinProjectDir1Minus(psi);
173  tmp2[rb[cb]] = shift(tmp, FORWARD, 1);
174  chi[rb[cb]] += spinReconstructDir1Minus(u[1]*tmp2);
175 
176  tmp[rb[otherCB]] = spinProjectDir2Minus(psi);
177  tmp2[rb[cb]] = shift(tmp, FORWARD, 2);
178  chi[rb[cb]] += spinReconstructDir2Minus(u[2]*tmp2);
179 
180  tmp[rb[otherCB]] = spinProjectDir3Minus(psi);
181  tmp2[rb[cb]] = shift(tmp, FORWARD, 3);
182  chi[rb[cb]] += spinReconstructDir3Minus(u[3]*tmp2);
183 
184 
185  tmp[rb[otherCB]] = adj(u[0])*spinProjectDir0Plus(psi);
186  tmp2[rb[cb]] = shift(tmp, BACKWARD, 0);
187  chi[rb[cb]] += spinReconstructDir0Plus(tmp2);
188 
189  tmp[rb[otherCB]] = adj(u[1])*spinProjectDir1Plus(psi);
190  tmp2[rb[cb]] = shift(tmp, BACKWARD, 1);
191  chi[rb[cb]] += spinReconstructDir1Plus(tmp2);
192 
193  tmp[rb[otherCB]] = adj(u[2])*spinProjectDir2Plus(psi);
194  tmp2[rb[cb]] = shift(tmp, BACKWARD, 2);
195  chi[rb[cb]] += spinReconstructDir2Plus(tmp2);
196 
197  tmp[rb[otherCB]] = adj(u[3])*spinProjectDir3Plus(psi);
198  tmp2[rb[cb]] = shift(tmp, BACKWARD, 3);
199  chi[rb[cb]] += spinReconstructDir3Plus(tmp2);
200 
201  }
202 
203 
204  break;
205 
206  case MINUS:
207  {
208 
209 
210  LatticeHalfFermion tmp;
211  LatticeHalfFermion tmp2;
212 
213  tmp[rb[otherCB]] = spinProjectDir0Plus(psi);
214  tmp2[rb[cb]] = shift(tmp, FORWARD, 0);
215  chi[rb[cb]] = spinReconstructDir0Plus(u[0]*tmp2);
216 
217  tmp[rb[otherCB]] = spinProjectDir1Plus(psi);
218  tmp2[rb[cb]] = shift(tmp, FORWARD, 1);
219  chi[rb[cb]] += spinReconstructDir1Plus(u[1]*tmp2);
220 
221  tmp[rb[otherCB]] = spinProjectDir2Plus(psi);
222  tmp2[rb[cb]] = shift(tmp, FORWARD, 2);
223  chi[rb[cb]] += spinReconstructDir2Plus(u[2]*tmp2);
224 
225  tmp[rb[otherCB]] = spinProjectDir3Plus(psi);
226  tmp2[rb[cb]] = shift(tmp, FORWARD, 3);
227  chi[rb[cb]] += spinReconstructDir3Plus(u[3]*tmp2);
228 
229 
230  tmp[rb[otherCB]] = adj(u[0])*spinProjectDir0Minus(psi);
231  tmp2[rb[cb]] = shift(tmp, BACKWARD, 0);
232  chi[rb[cb]] += spinReconstructDir0Minus(tmp2);
233 
234  tmp[rb[otherCB]] = adj(u[1])*spinProjectDir1Minus(psi);
235  tmp2[rb[cb]] = shift(tmp, BACKWARD, 1);
236  chi[rb[cb]] += spinReconstructDir1Minus(tmp2);
237 
238  tmp[rb[otherCB]] = adj(u[2])*spinProjectDir2Minus(psi);
239  tmp2[rb[cb]] = shift(tmp, BACKWARD, 2);
240  chi[rb[cb]] += spinReconstructDir2Minus(tmp2);
241 
242  tmp[rb[otherCB]] = adj(u[3])*spinProjectDir3Minus(psi);
243  tmp2 = shift(tmp, BACKWARD, 3);
244  chi[rb[cb]] += spinReconstructDir3Minus(tmp2);
245 
246  }
247 
248  break;
249  }
250 
251  getFermBC().modifyF(chi, QDP::rb[cb]);
252 #else
253  QDPIO::cerr<<"lwldslash_array_qdpopt_w: not implemented for NC!=3\n";
254  QDP_abort(13) ;
255 #endif
256 
257  END_CODE();
258  }
259 
260 } // End Namespace Chroma
261 
Primary include file for CHROMA library code.
Support class for fermion actions and linear operators.
Definition: state.h:94
Class for counted reference semantics.
Definition: handle.h:33
multi1d< LatticeColorMatrix > u
const FermBC< T, P, Q > & getFermBC() const
Return the fermion BC object for this linear operator.
int mu
Definition: cool.cc:24
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.
unsigned n
Definition: ldumul_w.cc:36
Wilson Dslash linear operator over arrays.
Nd
Definition: meslate.cc:74
Double tmp2
Definition: mesq.cc:30
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
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
Parameters for anisotropy.
Definition: aniso_io.h:24