QUDA  v1.1.0
A library for QCD on GPUs
dirac_twisted_clover.cpp
Go to the documentation of this file.
1 #include <dirac_quda.h>
2 #include <blas_quda.h>
3 #include <iostream>
4 #include <multigrid.h>
5 
6 #define NEW_DSLASH
7 
8 namespace quda {
9 
11  DiracWilson(param, nDim),
12  mu(param.mu),
14  clover(param.clover)
15  {
16  }
17 
20  mu(dirac.mu),
22  clover(dirac.clover)
23  {
24  }
25 
27 
29  {
30  if (&dirac != this)
31  {
33  clover = dirac.clover;
34  }
35 
36  return *this;
37  }
38 
40  {
41  Dirac::checkParitySpinor(out, in);
42 
43  if (out.Volume() != clover->VolumeCB())
44  errorQuda("Parity spinor volume %lu doesn't match clover checkboard volume %lu", out.Volume(), clover->VolumeCB());
45  }
46 
47  // Protected method for applying twist
49  {
50  checkParitySpinor(out, in);
51  ApplyTwistClover(out, in, *clover, kappa, mu, 0.0, parity, dagger, twistType);
52 
53  if (twistType == QUDA_TWIST_GAMMA5_INVERSE) flops += 1056ll*in.Volume();
54  else flops += 552ll*in.Volume();
55  }
56 
57 
58  // Public method to apply the twist
60  {
62  }
63 
65  {
66  // this would really just be a Wilson dslash (not actually instantiated at present)
67  errorQuda("Not implemented");
68  }
69 
71  const ColorSpinorField &x, const double &k) const
72  {
73 
74  if (in.TwistFlavor() == QUDA_TWIST_SINGLET) {
75  // k * D * in + (1 + i*2*mu*kappa*gamma_5) *x
76  ApplyTwistedClover(out, in, *gauge, *clover, k, 2 * mu * kappa, x, parity, dagger, commDim, profile);
77  flops += (1320ll + 552ll) * in.Volume();
78 
79  } else {
80  errorQuda("Non-degenerate operator is not implemented");
81  }
82  }
83 
85  {
86  checkFullSpinor(out, in);
87  if (in.TwistFlavor() != out.TwistFlavor())
88  errorQuda("Twist flavors %d %d don't match", in.TwistFlavor(), out.TwistFlavor());
89 
91  errorQuda("Twist flavor not set %d", in.TwistFlavor());
92  }
93 
95  profile);
96  flops += (1320ll + 552ll) * in.Volume();
97  }
98 
100  {
101  checkFullSpinor(out, in);
102  bool reset = newTmp(&tmp1, in);
103 
104  M(*tmp1, in);
105  Mdag(out, *tmp1);
106 
107  deleteTmp(&tmp1, reset);
108  }
109 
111  ColorSpinorField &b, const QudaSolutionType solType) const
112  {
113  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
114  errorQuda("Preconditioned solution requires a preconditioned solve_type");
115  }
116 
117  src = &b;
118  sol = &x;
119  }
120 
122  const QudaSolutionType solType) const
123  {
124  // do nothing
125  }
126 
128  {
129  Dirac::prefetch(mem_space, stream);
131  }
132 
134  double kappa, double mass, double mu, double mu_factor) const {
136  errorQuda("Wilson-type operators only support aggregation coarsening");
137 
138  double a = 2.0 * kappa * mu * T.Vectors().TwistFlavor();
140  }
141 
144  reverse(false)
145  {
146  }
147 
149  DiracTwistedClover(param, nDim),
150  reverse(false)
151  {
152  }
153 
155  {
156 
157  }
158 
160  {
161  if (&dirac != this) {
163  }
164  return *this;
165  }
166 
167  // Public method to apply the inverse twist
169  {
171  }
172 
173  // apply hopping term, then inverse twist: (A_ee^-1 D_eo) or (A_oo^-1 D_oe),
174  // and likewise for dagger: (D^dagger_eo D_ee^-1) or (D^dagger_oe A_oo^-1)
176  {
177  checkParitySpinor(in, out);
178  checkSpinorAlias(in, out);
179 
180  if (in.TwistFlavor() != out.TwistFlavor())
181  errorQuda("Twist flavors %d %d don't match", in.TwistFlavor(), out.TwistFlavor());
183  errorQuda("Twist flavor not set %d", in.TwistFlavor());
184 
185  bool symmetric = (matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
186  if (dagger && symmetric && !reverse) {
187  bool reset = newTmp(&tmp2, in);
188  TwistCloverInv(*tmp2, in, 1 - parity);
190  deleteTmp(&tmp2, reset);
191  } else {
192  ApplyTwistedCloverPreconditioned(out, in, *gauge, *clover, 1.0, -2.0 * kappa * mu, false, in, parity, dagger,
193  commDim, profile);
194  flops += (1320ll + 552ll) * in.Volume();
195  }
196  }
197 
198  // xpay version of the above
200  const ColorSpinorField &x, const double &k) const
201  {
202  checkParitySpinor(in, out);
203  checkSpinorAlias(in, out);
204  if (in.TwistFlavor() != out.TwistFlavor())
205  errorQuda("Twist flavors %d %d don't match", in.TwistFlavor(), out.TwistFlavor());
207  errorQuda("Twist flavor not set %d", in.TwistFlavor());
208 
209  bool symmetric = (matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
210  if (dagger && symmetric && !reverse) {
211  bool reset = newTmp(&tmp2, in);
212  TwistCloverInv(*tmp2, in, 1 - parity);
213  DiracWilson::DslashXpay(out, *tmp2, parity, x, k);
214  deleteTmp(&tmp2, reset);
215  } else {
216  ApplyTwistedCloverPreconditioned(out, in, *gauge, *clover, k, -2.0 * kappa * mu, true, x, parity, dagger, commDim,
217  profile);
218  flops += (1320ll + 552ll) * in.Volume();
219  }
220  }
221 
223  {
224  double kappa2 = -kappa*kappa;
225  bool reset = newTmp(&tmp1, in);
226 
227  bool symmetric =(matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
228  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
229  QudaParity parity[2] = {static_cast<QudaParity>((1 + odd_bit) % 2), static_cast<QudaParity>((0 + odd_bit) % 2)};
230 
231  if (in.TwistFlavor() == QUDA_TWIST_SINGLET) {
232  if (!symmetric) { // asymmetric preconditioning
233  Dslash(*tmp1, in, parity[0]);
234  DiracTwistedClover::DslashXpay(out, *tmp1, parity[1], in, kappa2);
235  } else if (!dagger) { // symmetric preconditioning
236  Dslash(*tmp1, in, parity[0]);
237  DslashXpay(out, *tmp1, parity[1], in, kappa2);
238  } else { // symmetric preconditioning, dagger
239  TwistCloverInv(out, in, parity[1]);
240  reverse = true;
241  Dslash(*tmp1, out, parity[0]);
242  reverse = false;
243  DiracWilson::DslashXpay(out, *tmp1, parity[1], in, kappa2);
244  }
245  } else { //Twist doublet
246  errorQuda("Non-degenerate operator is not implemented");
247  }
248 
249  deleteTmp(&tmp1, reset);
250  }
251 
253  {
254  // need extra temporary because of symmetric preconditioning dagger
255  bool reset = newTmp(&tmp2, in);
256  M(*tmp2, in);
257  Mdag(out, *tmp2);
258  deleteTmp(&tmp2, reset);
259  }
260 
262  ColorSpinorField &b, const QudaSolutionType solType) const
263  {
264  // we desire solution to preconditioned system
265  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
266  src = &b;
267  sol = &x;
268  return;
269  }
270 
271  bool reset = newTmp(&tmp1, b.Even());
272 
273  bool symmetric = (matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
274  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
275 
276  src = odd_bit ? &(x.Even()) : &(x.Odd());
277  sol = odd_bit ? &(x.Odd()) : &(x.Even());
278 
279  TwistCloverInv(symmetric ? *src : *tmp1, odd_bit ? b.Even() : b.Odd(), odd_bit ? QUDA_EVEN_PARITY : QUDA_ODD_PARITY);
280 
281  // we desire solution to full system
282  if (b.TwistFlavor() == QUDA_TWIST_SINGLET) {
283 
285  // src = A_ee^-1 (b_e + k D_eo A_oo^-1 b_o)
287  } else if (matpcType == QUDA_MATPC_ODD_ODD) {
288  // src = A_oo^-1 (b_o + k D_oe A_ee^-1 b_e)
291  // src = b_e + k D_eo A_oo^-1 b_o
294  // src = b_o + k D_oe A_ee^-1 b_e
296  } else {
297  errorQuda("MatPCType %d not valid for DiracTwistedCloverPC", matpcType);
298  }
299 
300  } else { // doublet:
301  errorQuda("Non-degenerate operator is not implemented");
302  } // end of doublet
303 
304  if (symmetric) TwistCloverInv(*src, *tmp1, odd_bit ? QUDA_ODD_PARITY : QUDA_EVEN_PARITY);
305 
306  // here we use final solution to store parity solution and parity source
307  // b is now up for grabs if we want
308 
309  deleteTmp(&tmp1, reset);
310  }
311 
313  const QudaSolutionType solType) const
314  {
315  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) { return; }
316 
317  checkFullSpinor(x, b);
318  bool reset = newTmp(&tmp1, b.Even());
319  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
320 
321  // create full solution
322  if (b.TwistFlavor() == QUDA_TWIST_SINGLET) {
324  // x_o = A_oo^-1 (b_o + k D_oe x_e)
327  // x_e = A_ee^-1 (b_e + k D_eo x_o)
329  } else {
330  errorQuda("MatPCType %d not valid for DiracTwistedCloverPC", matpcType);
331  }
332  } else { // twist doublet:
333  errorQuda("Non-degenerate operator is not implemented");
334  } // end of twist doublet...
335 
336  TwistCloverInv(odd_bit ? x.Even() : x.Odd(), *tmp1, odd_bit ? QUDA_EVEN_PARITY : QUDA_ODD_PARITY);
337  deleteTmp(&tmp1, reset);
338  }
339 
341  double kappa, double mass, double mu, double mu_factor) const {
343  errorQuda("Wilson-type operators only support aggregation coarsening");
344 
345  double a = -2.0 * kappa * mu * T.Vectors().TwistFlavor();
347  }
348 
350  {
351  Dirac::prefetch(mem_space, stream);
352 
353  bool symmetric = (matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
354  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
355  QudaParity parity[2] = {static_cast<QudaParity>((1 + odd_bit) % 2), static_cast<QudaParity>((0 + odd_bit) % 2)};
356 
357  if (symmetric) {
359  } else {
362  }
363  }
364 } // namespace quda
const ColorSpinorField & Odd() const
QudaTwistFlavorType TwistFlavor() const
const ColorSpinorField & Even() const
unsigned long long flops
Definition: dirac_quda.h:150
bool newTmp(ColorSpinorField **, const ColorSpinorField &) const
Definition: dirac.cpp:72
double kappa
Definition: dirac_quda.h:145
virtual void prefetch(QudaFieldLocation mem_space, qudaStream_t stream=0) const
If managed memory and prefetch is enabled, prefetch the gauge field and temporary spinors to the CPU ...
Definition: dirac.cpp:305
QudaMatPCType matpcType
Definition: dirac_quda.h:148
cudaGaugeField * gauge
Definition: dirac_quda.h:144
double mass
Definition: dirac_quda.h:146
void deleteTmp(ColorSpinorField **, const bool &reset) const
Definition: dirac.cpp:83
virtual void checkParitySpinor(const ColorSpinorField &, const ColorSpinorField &) const
Check parity spinors are usable (check geometry ?)
Definition: dirac.cpp:108
ColorSpinorField * tmp1
Definition: dirac_quda.h:151
TimeProfile profile
Definition: dirac_quda.h:161
QudaDagType dagger
Definition: dirac_quda.h:149
void checkSpinorAlias(const ColorSpinorField &, const ColorSpinorField &) const
check spinors do not alias
Definition: dirac.cpp:146
virtual void checkFullSpinor(const ColorSpinorField &, const ColorSpinorField &) const
check full spinors are compatible (check geometry ?)
Definition: dirac.cpp:138
ColorSpinorField * tmp2
Definition: dirac_quda.h:152
int commDim[QUDA_MAX_DIM]
Definition: dirac_quda.h:159
void Mdag(ColorSpinorField &out, const ColorSpinorField &in) const
Apply Mdag (daggered operator of M.
Definition: dirac.cpp:92
cudaCloverField * clover
Definition: dirac_quda.h:1049
void TwistClover(ColorSpinorField &out, const ColorSpinorField &in, const int parity) const
virtual void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
Xpay version of Dslash.
virtual void reconstruct(ColorSpinorField &x, const ColorSpinorField &b, const QudaSolutionType) const
virtual void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const
Apply MdagM operator which may be optimized.
DiracTwistedClover(const DiracTwistedClover &dirac)
void twistedCloverApply(ColorSpinorField &out, const ColorSpinorField &in, const QudaTwistGamma5Type twistType, const int parity) const
virtual void Dslash(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
apply 'dslash' operator for the DiracOp. This may be e.g. AD
void checkParitySpinor(const ColorSpinorField &, const ColorSpinorField &) const
Check parity spinors are usable (check geometry ?)
virtual void prepare(ColorSpinorField *&src, ColorSpinorField *&sol, ColorSpinorField &x, ColorSpinorField &b, const QudaSolutionType) const
void createCoarseOp(GaugeField &Y, GaugeField &X, const Transfer &T, double kappa, double mass, double mu, double mu_factor=0.) const
Create the coarse twisted-clover operator.
virtual void prefetch(QudaFieldLocation mem_space, qudaStream_t stream=0) const
If managed memory and prefetch is enabled, prefetch all relevant memory fields (gauge,...
virtual void M(ColorSpinorField &out, const ColorSpinorField &in) const
Apply M for the dirac op. E.g. the Schur Complement operator.
DiracTwistedClover & operator=(const DiracTwistedClover &dirac)
void createCoarseOp(GaugeField &Y, GaugeField &X, const Transfer &T, double kappa, double mass, double mu, double mu_factor=0.) const
Create the coarse even-odd preconditioned twisted-clover operator. Unlike the Wilson operator,...
void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const
Apply MdagM operator which may be optimized.
void M(ColorSpinorField &out, const ColorSpinorField &in) const
Apply M for the dirac op. E.g. the Schur Complement operator.
virtual void Dslash(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
apply 'dslash' operator for the DiracOp. This may be e.g. AD
void prepare(ColorSpinorField *&src, ColorSpinorField *&sol, ColorSpinorField &x, ColorSpinorField &b, const QudaSolutionType) const
void TwistCloverInv(ColorSpinorField &out, const ColorSpinorField &in, const int parity) const
virtual void prefetch(QudaFieldLocation mem_space, qudaStream_t stream=0) const
If managed memory and prefetch is enabled, prefetch all relevant memory fields (gauge,...
DiracTwistedCloverPC & operator=(const DiracTwistedCloverPC &dirac)
void reconstruct(ColorSpinorField &x, const ColorSpinorField &b, const QudaSolutionType) const
DiracTwistedCloverPC(const DiracTwistedCloverPC &dirac)
virtual void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
Xpay version of Dslash.
DiracWilson & operator=(const DiracWilson &dirac)
virtual void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
Xpay version of Dslash.
virtual void Dslash(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
apply 'dslash' operator for the DiracOp. This may be e.g. AD
size_t VolumeCB() const
QudaTransferType getTransferType() const
Definition: transfer.h:240
const ColorSpinorField & Vectors(QudaFieldLocation location=QUDA_INVALID_FIELD_LOCATION) const
Definition: transfer.h:209
void prefetch(QudaFieldLocation mem_space, qudaStream_t stream=0) const
If managed memory and prefetch is enabled, prefetch the clover, the norm field (as appropriate),...
double kappa
double mass
double epsilon
double mu
quda::mgarray< double > mu_factor
GaugeCovDev * dirac
Definition: covdev_test.cpp:42
QudaParity parity
Definition: covdev_test.cpp:40
@ QUDA_TWISTED_CLOVERPC_DIRAC
Definition: enum_quda.h:314
@ QUDA_TWISTED_CLOVER_DIRAC
Definition: enum_quda.h:313
enum QudaTwistGamma5Type_s QudaTwistGamma5Type
@ QUDA_EVEN_PARITY
Definition: enum_quda.h:284
@ QUDA_ODD_PARITY
Definition: enum_quda.h:284
@ QUDA_INVALID_PARITY
Definition: enum_quda.h:284
@ QUDA_TRANSFER_AGGREGATE
Definition: enum_quda.h:453
enum QudaSolutionType_s QudaSolutionType
enum QudaFieldLocation_s QudaFieldLocation
@ QUDA_MATPC_ODD_ODD_ASYMMETRIC
Definition: enum_quda.h:219
@ QUDA_MATPC_EVEN_EVEN_ASYMMETRIC
Definition: enum_quda.h:218
@ QUDA_MATPC_ODD_ODD
Definition: enum_quda.h:217
@ QUDA_MATPC_EVEN_EVEN
Definition: enum_quda.h:216
@ QUDA_MATPC_INVALID
Definition: enum_quda.h:220
@ QUDA_TWIST_GAMMA5_INVERSE
Definition: enum_quda.h:424
@ QUDA_TWIST_GAMMA5_DIRECT
Definition: enum_quda.h:423
@ QUDA_MATPC_SOLUTION
Definition: enum_quda.h:159
@ QUDA_MATPCDAG_MATPC_SOLUTION
Definition: enum_quda.h:161
@ QUDA_TWIST_SINGLET
Definition: enum_quda.h:400
@ QUDA_TWIST_NO
Definition: enum_quda.h:403
@ QUDA_TWIST_INVALID
Definition: enum_quda.h:404
enum QudaParity_s QudaParity
void ApplyTwistedCloverPreconditioned(ColorSpinorField &out, const ColorSpinorField &in, const GaugeField &U, const CloverField &C, double a, double b, bool xpay, const ColorSpinorField &x, int parity, bool dagger, const int *comm_override, TimeProfile &profile)
Driver for applying the preconditioned twisted-clover stencil.
void ApplyTwistClover(ColorSpinorField &out, const ColorSpinorField &in, const CloverField &clover, double kappa, double mu, double epsilon, int parity, int dagger, QudaTwistGamma5Type twist)
Apply twisted clover-matrix field to a color-spinor field.
qudaStream_t * stream
void ApplyTwistedClover(ColorSpinorField &out, const ColorSpinorField &in, const GaugeField &U, const CloverField &C, double a, double b, const ColorSpinorField &x, int parity, bool dagger, const int *comm_override, TimeProfile &profile)
Driver for applying the twisted-clover stencil.
void CoarseOp(GaugeField &Y, GaugeField &X, const Transfer &T, const cudaGaugeField &gauge, const cudaCloverField *clover, double kappa, double mass, double mu, double mu_factor, QudaDiracType dirac, QudaMatPCType matpc)
Coarse operator construction from a fine-grid operator (Wilson / Clover)
QudaGaugeParam param
Definition: pack_test.cpp:18
cudaStream_t qudaStream_t
Definition: quda_api.h:9
#define errorQuda(...)
Definition: util_quda.h:120