QUDA  v1.1.0
A library for QCD on GPUs
dirac_clover.cpp
Go to the documentation of this file.
1 #include <dirac_quda.h>
2 #include <blas_quda.h>
3 #include <multigrid.h>
4 
5 namespace quda {
6 
8 
10 
12 
14  {
15  if (&dirac != this) {
17  clover = dirac.clover;
18  }
19  return *this;
20  }
21 
23  {
24  Dirac::checkParitySpinor(out, in);
25 
26  if (out.Volume() != clover->VolumeCB()) {
27  errorQuda("Parity spinor volume %lu doesn't match clover checkboard volume %lu", out.Volume(), clover->VolumeCB());
28  }
29  }
30 
33  const QudaParity parity, const ColorSpinorField &x,
34  const double &k) const
35  {
36  checkParitySpinor(in, out);
37  checkSpinorAlias(in, out);
38 
39  ApplyWilsonClover(out, in, *gauge, *clover, k, x, parity, dagger, commDim, profile);
40  flops += 1872ll*in.Volume();
41  }
42 
43  // Public method to apply the clover term only
45  {
46  checkParitySpinor(in, out);
47 
48  ApplyClover(out, in, *clover, false, parity);
49  flops += 504ll*in.Volume();
50  }
51 
52  void DiracClover::M(ColorSpinorField &out, const ColorSpinorField &in) const
53  {
55  flops += 1872ll * in.Volume();
56  }
57 
59  {
60  checkFullSpinor(out, in);
61 
62  bool reset = newTmp(&tmp1, in);
63  checkFullSpinor(*tmp1, in);
64 
65  M(*tmp1, in);
66  Mdag(out, *tmp1);
67 
68  deleteTmp(&tmp1, reset);
69  }
70 
73  const QudaSolutionType solType) const
74  {
75  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
76  errorQuda("Preconditioned solution requires a preconditioned solve_type");
77  }
78 
79  src = &b;
80  sol = &x;
81  }
82 
84  const QudaSolutionType solType) const
85  {
86  // do nothing
87  }
88 
90  double kappa, double mass, double mu, double mu_factor) const {
92  errorQuda("Wilson-type operators only support aggregation coarsening");
93 
94  double a = 2.0 * kappa * mu * T.Vectors().TwistFlavor();
96  }
97 
99  {
100  Dirac::prefetch(mem_space, stream);
102  }
103 
104  /*******
105  * DiracCloverPC Starts here
106  *******/
109  {
110  // For the preconditioned operator, we need to check that the inverse of the clover term is present
111  if (!clover->cloverInv) errorQuda("Clover inverse required for DiracCloverPC");
112  }
113 
115 
117 
119  {
120  if (&dirac != this) {
122  }
123  return *this;
124  }
125 
126  // Public method
128  const QudaParity parity) const
129  {
130  checkParitySpinor(in, out);
131 
132  ApplyClover(out, in, *clover, true, parity);
133  flops += 504ll*in.Volume();
134  }
135 
136  // apply hopping term, then clover: (A_ee^-1 D_eo) or (A_oo^-1 D_oe),
137  // and likewise for dagger: (A_ee^-1 D^dagger_eo) or (A_oo^-1 D^dagger_oe)
138  // NOTE - this isn't Dslash dagger since order should be reversed!
140  const QudaParity parity) const
141  {
142  checkParitySpinor(in, out);
143  checkSpinorAlias(in, out);
144 
146  flops += 1824ll*in.Volume();
147  }
148 
149  // xpay version of the above
151  const QudaParity parity, const ColorSpinorField &x,
152  const double &k) const
153  {
154  checkParitySpinor(in, out);
155  checkSpinorAlias(in, out);
156 
158  flops += 1872ll*in.Volume();
159  }
160 
161  // Apply the even-odd preconditioned clover-improved Dirac operator
163  {
164  double kappa2 = -kappa*kappa;
165  bool reset1 = newTmp(&tmp1, in);
166 
167  bool symmetric =(matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
168  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
169  QudaParity parity[2] = {static_cast<QudaParity>((1 + odd_bit) % 2), static_cast<QudaParity>((0 + odd_bit) % 2)};
170 
171  if (!symmetric) {
172 
173  // No need to change order of calls for dagger
174  // because the asymmetric operator is actually symmetric
175  // A_oo -D_oe A^{-1}_ee D_eo -> A_oo -D^\dag_oe A^{-1}_ee D^\dag_eo
176  // the pieces in Dslash and DslashXPay respect the dagger
177 
178  // DiracCloverPC::Dslash applies A^{-1}Dslash
179  Dslash(*tmp1, in, parity[0]);
180  // DiracClover::DslashXpay applies (A - kappa^2 D)
181  DiracClover::DslashXpay(out, *tmp1, parity[1], in, kappa2);
182  } else if (!dagger) { // symmetric preconditioning
183  // We need two cases because M = 1-ADAD and M^\dag = 1-D^\dag A D^dag A
184  // where A is actually a clover inverse.
185 
186  // This is the non-dag case: AD
187  Dslash(*tmp1, in, parity[0]);
188 
189  // Then x + AD (AD)
190  DslashXpay(out, *tmp1, parity[1], in, kappa2);
191  } else { // symmetric preconditioning, dagger
192 
193  // This is the dagger: 1 - DADA
194  // i) Apply A
195  CloverInv(out, in, parity[1]);
196  // ii) Apply A D => ADA
197  Dslash(*tmp1, out, parity[0]);
198  // iii) Apply x + D(ADA)
199  DiracWilson::DslashXpay(out, *tmp1, parity[1], in, kappa2);
200  }
201 
202  deleteTmp(&tmp1, reset1);
203  }
204 
206  {
207  // need extra temporary because of symmetric preconditioning dagger
208  // and for multi-gpu the input and output fields cannot alias
209  bool reset = newTmp(&tmp2, in);
210  M(*tmp2, in);
211  Mdag(out, *tmp2);
212  deleteTmp(&tmp2, reset);
213  }
214 
217  const QudaSolutionType solType) const
218  {
219  // we desire solution to preconditioned system
220  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
221  src = &b;
222  sol = &x;
223  return;
224  }
225 
226  bool reset = newTmp(&tmp1, b.Even());
227 
228  // we desire solution to full system
230  // src = A_ee^-1 (b_e + k D_eo A_oo^-1 b_o)
231  src = &(x.Odd());
232  CloverInv(*src, b.Odd(), QUDA_ODD_PARITY);
235  sol = &(x.Even());
236  } else if (matpcType == QUDA_MATPC_ODD_ODD) {
237  // src = A_oo^-1 (b_o + k D_oe A_ee^-1 b_e)
238  src = &(x.Even());
239  CloverInv(*src, b.Even(), QUDA_EVEN_PARITY);
241  CloverInv(*src, *tmp1, QUDA_ODD_PARITY);
242  sol = &(x.Odd());
244  // src = b_e + k D_eo A_oo^-1 b_o
245  src = &(x.Odd());
246  CloverInv(*tmp1, b.Odd(), QUDA_ODD_PARITY); // safe even when *tmp1 = b.odd
248  sol = &(x.Even());
250  // src = b_o + k D_oe A_ee^-1 b_e
251  src = &(x.Even());
252  CloverInv(*tmp1, b.Even(), QUDA_EVEN_PARITY); // safe even when *tmp1 = b.even
254  sol = &(x.Odd());
255  } else {
256  errorQuda("MatPCType %d not valid for DiracCloverPC", matpcType);
257  }
258 
259  // here we use final solution to store parity solution and parity source
260  // b is now up for grabs if we want
261 
262  deleteTmp(&tmp1, reset);
263 
264  }
265 
267  const QudaSolutionType solType) const
268  {
269  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
270  return;
271  }
272 
273  checkFullSpinor(x, b);
274 
275  bool reset = newTmp(&tmp1, b.Even());
276 
277  // create full solution
278 
281  // x_o = A_oo^-1 (b_o + k D_oe x_e)
284  } else if (matpcType == QUDA_MATPC_ODD_ODD ||
286  // x_e = A_ee^-1 (b_e + k D_eo x_o)
289  } else {
290  errorQuda("MatPCType %d not valid for DiracCloverPC", matpcType);
291  }
292 
293  deleteTmp(&tmp1, reset);
294 
295  }
296 
298  double kappa, double mass, double mu, double mu_factor) const {
300  errorQuda("Wilson-type operators only support aggregation coarsening");
301 
302  double a = - 2.0 * kappa * mu * T.Vectors().TwistFlavor();
304  }
305 
307  {
308  Dirac::prefetch(mem_space, stream);
309 
310  bool symmetric = (matpcType == QUDA_MATPC_EVEN_EVEN || matpcType == QUDA_MATPC_ODD_ODD) ? true : false;
311  int odd_bit = (matpcType == QUDA_MATPC_ODD_ODD || matpcType == QUDA_MATPC_ODD_ODD_ASYMMETRIC) ? 1 : 0;
312  QudaParity parity[2] = {static_cast<QudaParity>((1 + odd_bit) % 2), static_cast<QudaParity>((0 + odd_bit) % 2)};
313 
314  if (symmetric) {
316  } else {
319  }
320  }
321 
322 } // namespace quda
const ColorSpinorField & Odd() const
QudaTwistFlavorType TwistFlavor() const
const ColorSpinorField & Even() const
virtual void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const
Apply MdagM operator which may be optimized.
DiracClover & operator=(const DiracClover &dirac)
virtual void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
virtual ~DiracClover()
DiracClover(const DiracParam &param)
Definition: dirac_clover.cpp:7
virtual void reconstruct(ColorSpinorField &x, const ColorSpinorField &b, const QudaSolutionType) const
virtual void M(ColorSpinorField &out, const ColorSpinorField &in) const
Apply M for the dirac op. E.g. the Schur Complement operator.
void checkParitySpinor(const ColorSpinorField &, const ColorSpinorField &) const
Check parity spinors are usable (check geometry ?)
void Clover(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
void createCoarseOp(GaugeField &Y, GaugeField &X, const Transfer &T, double kappa, double mass=0., double mu=0., double mu_factor=0.) const
Create the coarse 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,...
cudaCloverField * clover
Definition: dirac_quda.h:470
virtual void prepare(ColorSpinorField *&src, ColorSpinorField *&sol, ColorSpinorField &x, ColorSpinorField &b, const QudaSolutionType) const
void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const
Apply MdagM operator which may be optimized.
void Dslash(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
apply 'dslash' operator for the DiracOp. This may be e.g. AD
void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
void M(ColorSpinorField &out, const ColorSpinorField &in) const
Apply M for the dirac op. E.g. the Schur Complement operator.
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=0., double mu=0., double mu_factor=0.) const
Create the coarse even-odd preconditioned clover operator. Unlike the Wilson operator,...
DiracCloverPC & operator=(const DiracCloverPC &dirac)
void CloverInv(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity 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,...
void reconstruct(ColorSpinorField &x, const ColorSpinorField &b, const QudaSolutionType) const
DiracCloverPC(const DiracParam &param)
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
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.
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 mu
quda::mgarray< double > mu_factor
GaugeCovDev * dirac
Definition: covdev_test.cpp:42
QudaParity parity
Definition: covdev_test.cpp:40
@ QUDA_CLOVER_DIRAC
Definition: enum_quda.h:293
@ QUDA_CLOVERPC_DIRAC
Definition: enum_quda.h:294
@ 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_MATPC_SOLUTION
Definition: enum_quda.h:159
@ QUDA_MATPCDAG_MATPC_SOLUTION
Definition: enum_quda.h:161
enum QudaParity_s QudaParity
void ApplyClover(ColorSpinorField &out, const ColorSpinorField &in, const CloverField &clover, bool inverse, int parity)
Apply clover-matrix field to a color-spinor field.
qudaStream_t * stream
void ApplyWilsonClover(ColorSpinorField &out, const ColorSpinorField &in, const GaugeField &U, const CloverField &A, double kappa, const ColorSpinorField &x, int parity, bool dagger, const int *comm_override, TimeProfile &profile)
Driver for applying the Wilson-clover stencil.
void ApplyWilsonCloverPreconditioned(ColorSpinorField &out, const ColorSpinorField &in, const GaugeField &U, const CloverField &A, double kappa, const ColorSpinorField &x, int parity, bool dagger, const int *comm_override, TimeProfile &profile)
Driver for applying the preconditioned Wilson-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