QUDA  1.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
dirac.cpp
Go to the documentation of this file.
1 #include <dirac_quda.h>
2 #include <dslash_quda.h>
3 #include <blas_quda.h>
4 
5 #include <iostream>
6 
7 namespace quda {
8 
9  // FIXME: At the moment, it's unsafe for more than one Dirac operator to be active unless
10  // they all have the same volume, etc. (used to initialize the various CUDA constants).
11 
13  gauge(param.gauge),
14  kappa(param.kappa),
15  mass(param.mass),
16  laplace3D(param.laplace3D),
17  matpcType(param.matpcType),
18  dagger(param.dagger),
19  flops(0),
20  tmp1(param.tmp1),
21  tmp2(param.tmp2),
22  type(param.type),
23  halo_precision(param.halo_precision),
24  profile("Dirac", false)
25  {
26  for (int i=0; i<4; i++) commDim[i] = param.commDim[i];
27  }
28 
30  gauge(dirac.gauge),
31  kappa(dirac.kappa),
32  laplace3D(dirac.laplace3D),
33  matpcType(dirac.matpcType),
34  dagger(dirac.dagger),
35  flops(0),
36  tmp1(dirac.tmp1),
37  tmp2(dirac.tmp2),
38  type(dirac.type),
39  halo_precision(dirac.halo_precision),
40  profile("Dirac", false)
41  {
42  for (int i=0; i<4; i++) commDim[i] = dirac.commDim[i];
43  }
44 
47  }
48 
50  {
51  if (&dirac != this) {
52  gauge = dirac.gauge;
53  kappa = dirac.kappa;
54  laplace3D = dirac.laplace3D;
55  matpcType = dirac.matpcType;
56  dagger = dirac.dagger;
57  flops = 0;
58  tmp1 = dirac.tmp1;
59  tmp2 = dirac.tmp2;
60 
61  for (int i=0; i<4; i++) commDim[i] = dirac.commDim[i];
62 
63  profile = dirac.profile;
64 
65  if (type != dirac.type) errorQuda("Trying to copy between incompatible types %d %d", type, dirac.type);
66  }
67  return *this;
68  }
69 
71  if (*tmp) return false;
73  param.create = QUDA_ZERO_FIELD_CREATE; // need to zero elements else padded region will be junk
74 
75  if (typeid(a) == typeid(cudaColorSpinorField)) *tmp = new cudaColorSpinorField(a, param);
76  else *tmp = new cpuColorSpinorField(param);
77 
78  return true;
79  }
80 
81  void Dirac::deleteTmp(ColorSpinorField **a, const bool &reset) const {
82  if (reset) {
83  delete *a;
84  *a = NULL;
85  }
86  }
87 
88 #define flip(x) (x) = ((x) == QUDA_DAG_YES ? QUDA_DAG_NO : QUDA_DAG_YES)
89 
91  {
92  flip(dagger);
93  M(out, in);
94  flip(dagger);
95  }
96 
98  {
99  flip(dagger);
100  MdagM(out, in);
101  flip(dagger);
102  }
103 
104 #undef flip
105 
107  {
109  in.Nspin() == 4) {
110  errorQuda("CUDA Dirac operator requires UKQCD basis, out = %d, in = %d",
111  out.GammaBasis(), in.GammaBasis());
112  }
113 
114  if (in.Precision() != out.Precision()) {
115  errorQuda("Input precision %d and output spinor precision %d don't match in dslash_quda",
116  in.Precision(), out.Precision());
117  }
118 
119  if (in.Stride() != out.Stride()) {
120  errorQuda("Input %d and output %d spinor strides don't match in dslash_quda",
121  in.Stride(), out.Stride());
122  }
123 
125  errorQuda("ColorSpinorFields are not single parity: in = %d, out = %d",
126  in.SiteSubset(), out.SiteSubset());
127  }
128 
129  if (!static_cast<const cudaColorSpinorField&>(in).isNative()) errorQuda("Input field is not in native order");
130  if (!static_cast<const cudaColorSpinorField&>(out).isNative()) errorQuda("Output field is not in native order");
131 
132  if (out.Ndim() != 5) {
133  if ((out.Volume() != gauge->Volume() && out.SiteSubset() == QUDA_FULL_SITE_SUBSET) ||
134  (out.Volume() != gauge->VolumeCB() && out.SiteSubset() == QUDA_PARITY_SITE_SUBSET) ) {
135  errorQuda("Spinor volume %d doesn't match gauge volume %d", out.Volume(), gauge->VolumeCB());
136  }
137  } else {
138  // Domain wall fermions, compare 4d volumes not 5d
139  if ((out.Volume()/out.X(4) != gauge->Volume() && out.SiteSubset() == QUDA_FULL_SITE_SUBSET) ||
140  (out.Volume()/out.X(4) != gauge->VolumeCB() && out.SiteSubset() == QUDA_PARITY_SITE_SUBSET) ) {
141  errorQuda("Spinor volume %d doesn't match gauge volume %d", out.Volume(), gauge->VolumeCB());
142  }
143  }
144  }
145 
147  {
149  errorQuda("ColorSpinorFields are not full fields: in = %d, out = %d",
150  in.SiteSubset(), out.SiteSubset());
151  }
152  }
153 
155  if (a.V() == b.V()) errorQuda("Aliasing pointers");
156  }
157 
158  // Dirac operator factory
160  {
161  if (param.type == QUDA_WILSON_DIRAC) {
162  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracWilson operator\n");
163  return new DiracWilson(param);
164  } else if (param.type == QUDA_WILSONPC_DIRAC) {
165  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracWilsonPC operator\n");
166  return new DiracWilsonPC(param);
167  } else if (param.type == QUDA_CLOVER_DIRAC) {
168  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracClover operator\n");
169  return new DiracClover(param);
170  } else if (param.type == QUDA_CLOVERPC_DIRAC) {
171  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracCloverPC operator\n");
172  return new DiracCloverPC(param);
173  } else if (param.type == QUDA_DOMAIN_WALL_DIRAC) {
174  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracDomainWall operator\n");
175  return new DiracDomainWall(param);
176  } else if (param.type == QUDA_DOMAIN_WALLPC_DIRAC) {
177  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracDomainWallPC operator\n");
178  return new DiracDomainWallPC(param);
179  } else if (param.type == QUDA_DOMAIN_WALL_4D_DIRAC) {
180  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracDomainWall4D operator\n");
181  return new DiracDomainWall4D(param);
182  } else if (param.type == QUDA_DOMAIN_WALL_4DPC_DIRAC) {
183  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracDomainWall4DPC operator\n");
184  return new DiracDomainWall4DPC(param);
185  } else if (param.type == QUDA_MOBIUS_DOMAIN_WALL_DIRAC) {
186  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracMobius operator\n");
187  return new DiracMobius(param);
188  } else if (param.type == QUDA_MOBIUS_DOMAIN_WALLPC_DIRAC) {
189  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracMobiusPC operator\n");
190  return new DiracMobiusPC(param);
191  } else if (param.type == QUDA_STAGGERED_DIRAC) {
192  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracStaggered operator\n");
193  return new DiracStaggered(param);
194  } else if (param.type == QUDA_STAGGEREDPC_DIRAC) {
195  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracStaggeredPC operator\n");
196  return new DiracStaggeredPC(param);
197  } else if (param.type == QUDA_ASQTAD_DIRAC) {
198  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracImprovedStaggered operator\n");
199  return new DiracImprovedStaggered(param);
200  } else if (param.type == QUDA_ASQTADPC_DIRAC) {
201  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracImprovedStaggeredPC operator\n");
202  return new DiracImprovedStaggeredPC(param);
203  } else if (param.type == QUDA_TWISTED_CLOVER_DIRAC) {
204  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracTwistedClover operator (%d flavor(s))\n", param.Ls);
205  if (param.Ls == 1) {
206  return new DiracTwistedClover(param, 4);
207  } else {
208  errorQuda("Cannot create DiracTwistedClover operator for %d flavors\n", param.Ls);
209  }
210  } else if (param.type == QUDA_TWISTED_CLOVERPC_DIRAC) {
211  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracTwistedCloverPC operator (%d flavor(s))\n", param.Ls);
212  if (param.Ls == 1) {
213  return new DiracTwistedCloverPC(param, 4);
214  } else {
215  errorQuda("Cannot create DiracTwistedCloverPC operator for %d flavors\n", param.Ls);
216  }
217  } else if (param.type == QUDA_TWISTED_MASS_DIRAC) {
218  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracTwistedMass operator (%d flavor(s))\n", param.Ls);
219  if (param.Ls == 1) return new DiracTwistedMass(param, 4);
220  else return new DiracTwistedMass(param, 5);
221  } else if (param.type == QUDA_TWISTED_MASSPC_DIRAC) {
223  printfQuda("Creating a DiracTwistedMassPC operator (%d flavor(s))\n", param.Ls);
224  if (param.Ls == 1)
225  return new DiracTwistedMassPC(param, 4);
226  else
227  return new DiracTwistedMassPC(param, 5);
228  } else if (param.type == QUDA_COARSE_DIRAC) {
229  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracCoarse operator\n");
230  return new DiracCoarse(param);
231  } else if (param.type == QUDA_COARSEPC_DIRAC) {
232  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a DiracCoarsePC operator\n");
233  return new DiracCoarsePC(param);
234  } else if (param.type == QUDA_GAUGE_COVDEV_DIRAC) {
235  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a GaugeCovDev operator\n");
236  return new GaugeCovDev(param);
237  } else if (param.type == QUDA_GAUGE_LAPLACE_DIRAC) {
238  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a GaugeLaplace operator\n");
239  return new GaugeLaplace(param);
240  } else if (param.type == QUDA_GAUGE_LAPLACEPC_DIRAC) {
241  if (getVerbosity() >= QUDA_DEBUG_VERBOSE) printfQuda("Creating a GaugeLaplacePC operator\n");
242  return new GaugeLaplacePC(param);
243  } else {
244  errorQuda("Unsupported Dirac type %d", param.type);
245  }
246 
247  return nullptr;
248  }
249 
250  // Count the number of stencil applications per dslash application.
252  {
253  int steps = 0;
254  switch (type)
255  {
256  case QUDA_COARSE_DIRAC: // single fused operator
259  steps = 1;
260  break;
261  case QUDA_WILSON_DIRAC:
262  case QUDA_CLOVER_DIRAC:
266  case QUDA_ASQTAD_DIRAC:
269  steps = 2; // For D_{eo} and D_{oe} piece.
270  break;
271  case QUDA_WILSONPC_DIRAC:
272  case QUDA_CLOVERPC_DIRAC:
277  case QUDA_ASQTADPC_DIRAC:
280  case QUDA_COARSEPC_DIRAC:
282  steps = 2;
283  break;
284  default:
285  errorQuda("Unsupported Dslash type %d.\n", type);
286  steps = 0;
287  break;
288  }
289 
290  return steps;
291  }
292 
293 } // namespace quda
cudaColorSpinorField * tmp2
Dirac(const DiracParam &param)
Definition: dirac.cpp:12
Even-odd preconditioned Gauge Laplace operator.
Definition: dirac_quda.h:1047
QudaDiracType type
Definition: dirac_quda.h:124
unsigned long long flops
Definition: dirac_quda.h:121
cudaColorSpinorField * tmp1
QudaVerbosity getVerbosity()
Definition: util_quda.cpp:21
double kappa
Definition: test_util.cpp:1647
#define errorQuda(...)
Definition: util_quda.h:121
#define flip(x)
Definition: dirac.cpp:88
cudaGaugeField * gauge
Definition: dirac_quda.h:115
virtual void checkFullSpinor(const ColorSpinorField &, const ColorSpinorField &) const
Definition: dirac.cpp:146
cudaColorSpinorField * tmp
Definition: covdev_test.cpp:44
void deleteTmp(ColorSpinorField **, const bool &reset) const
Definition: dirac.cpp:81
QudaGammaBasis GammaBasis() const
TimeProfile profile
Definition: dirac_quda.h:132
int getStencilSteps() const
Definition: dirac.cpp:251
virtual ~Dirac()
Definition: dirac.cpp:45
virtual void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const =0
QudaGaugeParam param
Definition: pack_test.cpp:17
bool newTmp(ColorSpinorField **, const ColorSpinorField &) const
Definition: dirac.cpp:70
int commDim[QUDA_MAX_DIM]
Definition: dirac_quda.h:130
void MMdag(ColorSpinorField &out, const ColorSpinorField &in) const
Definition: dirac.cpp:97
void checkSpinorAlias(const ColorSpinorField &, const ColorSpinorField &) const
Definition: dirac.cpp:154
int commDim[QUDA_MAX_DIM]
Definition: dirac_quda.h:44
cpuColorSpinorField * in
void Print()
Definition: timer.cpp:7
QudaSiteSubset SiteSubset() const
QudaDiracType type
Definition: dirac_quda.h:22
QudaDagType dagger
Definition: dirac_quda.h:120
int Volume() const
double kappa
Definition: dirac_quda.h:116
QudaMatPCType matpcType
Definition: dirac_quda.h:119
void Mdag(ColorSpinorField &out, const ColorSpinorField &in) const
Definition: dirac.cpp:90
virtual void M(ColorSpinorField &out, const ColorSpinorField &in) const =0
Dirac & operator=(const Dirac &dirac)
Definition: dirac.cpp:49
GaugeCovDev * dirac
Definition: covdev_test.cpp:73
cpuColorSpinorField * out
int laplace3D
Definition: test_util.cpp:1622
Full Gauge Laplace operator. Although not a Dirac operator per se, it&#39;s a linear operator so it&#39;s con...
Definition: dirac_quda.h:1022
#define printfQuda(...)
Definition: util_quda.h:115
int VolumeCB() const
unsigned long long flops
Definition: blas_quda.cu:22
const int * X() const
static Dirac * create(const DiracParam &param)
Definition: dirac.cpp:159
virtual void checkParitySpinor(const ColorSpinorField &, const ColorSpinorField &) const
Definition: dirac.cpp:106
QudaPrecision Precision() const
QudaDagType dagger
Definition: test_util.cpp:1620
Full Covariant Derivative operator. Although not a Dirac operator per se, it&#39;s a linear operator so i...
Definition: dirac_quda.h:1069
ColorSpinorField * tmp2
Definition: dirac_quda.h:123
ColorSpinorField * tmp1
Definition: dirac_quda.h:122