QUDA  v1.1.0
A library for QCD on GPUs
dirac_domain_wall.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <dirac_quda.h>
3 #include <blas_quda.h>
4 
5 namespace quda {
6 
8  DiracWilson(param, 5),
9  m5(param.m5),
10  kappa5(0.5 / (5.0 + m5)),
11  Ls(param.Ls)
12  {
13  }
14 
17  m5(dirac.m5),
18  kappa5(0.5 / (5.0 + m5)),
19  Ls(dirac.Ls)
20  {
21  }
22 
24 
26  {
27  if (&dirac != this) {
29  m5 = dirac.m5;
30  kappa5 = dirac.kappa5;
31  }
32  return *this;
33  }
34 
36  {
37  if (in.Ndim() != 5 || out.Ndim() != 5) errorQuda("Wrong number of dimensions\n");
38  if (in.X(4) != Ls) errorQuda("Expected Ls = %d, not %d\n", Ls, in.X(4));
39  if (out.X(4) != Ls) errorQuda("Expected Ls = %d, not %d\n", Ls, out.X(4));
40  }
41 
43  const QudaParity parity) const
44  {
45  checkDWF(out, in);
46  checkParitySpinor(in, out);
47  checkSpinorAlias(in, out);
48 
49  ApplyDomainWall5D(out, in, *gauge, 0.0, mass, in, parity, dagger, commDim, profile);
50 
51  long long Ls = in.X(4);
52  long long bulk = (Ls-2)*(in.Volume()/Ls);
53  long long wall = 2*in.Volume()/Ls;
54  flops += 1320LL*(long long)in.Volume() + 96LL*bulk + 120LL*wall;
55  }
56 
58  const QudaParity parity, const ColorSpinorField &x,
59  const double &k) const
60  {
61  checkDWF(out, in);
62  checkParitySpinor(in, out);
63  checkSpinorAlias(in, out);
64 
65  ApplyDomainWall5D(out, in, *gauge, k, mass, x, parity, dagger, commDim, profile);
66 
67  long long Ls = in.X(4);
68  long long bulk = (Ls-2)*(in.Volume()/Ls);
69  long long wall = 2*in.Volume()/Ls;
70  flops += (1320LL+48LL)*(long long)in.Volume() + 96LL*bulk + 120LL*wall;
71  }
72 
74  {
75  checkFullSpinor(out, in);
76 
78 
79  long long Ls = in.X(4);
80  long long bulk = (Ls - 2) * (in.Volume() / Ls);
81  long long wall = 2 * in.Volume() / Ls;
82  flops += (1320LL + 48LL) * (long long)in.Volume() + 96LL * bulk + 120LL * wall;
83  }
84 
86  {
87  checkFullSpinor(out, in);
88 
89  bool reset = newTmp(&tmp1, in);
90 
91  M(*tmp1, in);
92  Mdag(out, *tmp1);
93 
94  deleteTmp(&tmp1, reset);
95  }
96 
99  const QudaSolutionType solType) const
100  {
101  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
102  errorQuda("Preconditioned solution requires a preconditioned solve_type");
103  }
104 
105  src = &b;
106  sol = &x;
107  }
108 
110  const QudaSolutionType solType) const
111  {
112  // do nothing
113  }
114 
117  {
118 
119  }
120 
123  {
124 
125  }
126 
128  {
129 
130  }
131 
133  {
134  if (&dirac != this) {
136  }
137 
138  return *this;
139  }
140 
141  // Apply the even-odd preconditioned clover-improved Dirac operator
143  {
144  checkDWF(out, in);
145  double kappa2 = -kappa5*kappa5;
146 
147  bool reset = newTmp(&tmp1, in);
148 
150  Dslash(*tmp1, in, QUDA_ODD_PARITY);
151  DslashXpay(out, *tmp1, QUDA_EVEN_PARITY, in, kappa2);
152  } else if (matpcType == QUDA_MATPC_ODD_ODD) {
154  DslashXpay(out, *tmp1, QUDA_ODD_PARITY, in, kappa2);
155  } else {
156  errorQuda("MatPCType %d not valid for DiracDomainWallPC", matpcType);
157  }
158 
159  deleteTmp(&tmp1, reset);
160  }
161 
163  {
164  //M(out, in);
165  //Mdag(out, out);
166  bool reset = newTmp(&tmp2, in);
167  M(*tmp2, in);
168  Mdag(out, *tmp2);
169  deleteTmp(&tmp2, reset);
170  }
171 
174  const QudaSolutionType solType) const
175  {
176  // we desire solution to preconditioned system
177  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
178  src = &b;
179  sol = &x;
180  } else {
181  // we desire solution to full system
183  // src = b_e + k D_eo b_o
184  DslashXpay(x.Odd(), b.Odd(), QUDA_EVEN_PARITY, b.Even(), kappa5);
185  src = &(x.Odd());
186  sol = &(x.Even());
187  } else if (matpcType == QUDA_MATPC_ODD_ODD) {
188  // src = b_o + k D_oe b_e
189  DslashXpay(x.Even(), b.Even(), QUDA_ODD_PARITY, b.Odd(), kappa5);
190  src = &(x.Even());
191  sol = &(x.Odd());
192  } else {
193  errorQuda("MatPCType %d not valid for DiracDomainWallPC", matpcType);
194  }
195  // here we use final solution to store parity solution and parity source
196  // b is now up for grabs if we want
197  }
198 
199  }
200 
202  const QudaSolutionType solType) const
203  {
204  if (solType == QUDA_MATPC_SOLUTION || solType == QUDA_MATPCDAG_MATPC_SOLUTION) {
205  return;
206  }
207 
208  // create full solution
209 
210  checkFullSpinor(x, b);
212  // x_o = b_o + k D_oe x_e
213  DslashXpay(x.Odd(), x.Even(), QUDA_ODD_PARITY, b.Odd(), kappa5);
214  } else if (matpcType == QUDA_MATPC_ODD_ODD) {
215  // x_e = b_e + k D_eo x_o
216  DslashXpay(x.Even(), x.Odd(), QUDA_EVEN_PARITY, b.Even(), kappa5);
217  } else {
218  errorQuda("MatPCType %d not valid for DiracDomainWallPC", matpcType);
219  }
220  }
221 
222 
223 } // namespace quda
const ColorSpinorField & Odd() const
const ColorSpinorField & Even() const
const int * X() const
void Dslash(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity) const
apply 'dslash' operator for the DiracOp. This may be e.g. AD
DiracDomainWall & operator=(const DiracDomainWall &dirac)
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.
DiracDomainWall(const DiracParam &param)
void DslashXpay(ColorSpinorField &out, const ColorSpinorField &in, const QudaParity parity, const ColorSpinorField &x, const double &k) const
Xpay version of Dslash.
virtual void prepare(ColorSpinorField *&src, ColorSpinorField *&sol, ColorSpinorField &x, 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 checkDWF(const ColorSpinorField &out, const ColorSpinorField &in) const
void MdagM(ColorSpinorField &out, const ColorSpinorField &in) const
Apply MdagM operator which may be optimized.
void reconstruct(ColorSpinorField &x, const ColorSpinorField &b, const QudaSolutionType) const
DiracDomainWallPC(const DiracParam &param)
DiracDomainWallPC & operator=(const DiracDomainWallPC &dirac)
void prepare(ColorSpinorField *&src, ColorSpinorField *&sol, ColorSpinorField &x, ColorSpinorField &b, const QudaSolutionType) const
void M(ColorSpinorField &out, const ColorSpinorField &in) const
Apply M for the dirac op. E.g. the Schur Complement operator.
unsigned long long flops
Definition: dirac_quda.h:150
bool newTmp(ColorSpinorField **, const ColorSpinorField &) const
Definition: dirac.cpp:72
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)
double m5
GaugeCovDev * dirac
Definition: covdev_test.cpp:42
QudaParity parity
Definition: covdev_test.cpp:40
@ QUDA_EVEN_PARITY
Definition: enum_quda.h:284
@ QUDA_ODD_PARITY
Definition: enum_quda.h:284
@ QUDA_INVALID_PARITY
Definition: enum_quda.h:284
enum QudaSolutionType_s QudaSolutionType
@ QUDA_MATPC_ODD_ODD
Definition: enum_quda.h:217
@ QUDA_MATPC_EVEN_EVEN
Definition: enum_quda.h:216
@ QUDA_MATPC_SOLUTION
Definition: enum_quda.h:159
@ QUDA_MATPCDAG_MATPC_SOLUTION
Definition: enum_quda.h:161
enum QudaParity_s QudaParity
int Ls
Definition: host_utils.cpp:48
double kappa5
Definition: host_utils.cpp:51
void ApplyDomainWall5D(ColorSpinorField &out, const ColorSpinorField &in, const GaugeField &U, double a, double m_f, const ColorSpinorField &x, int parity, bool dagger, const int *comm_override, TimeProfile &profile)
Driver for applying the Domain-wall 5-d stencil to a 5-d vector with 5-d preconditioned data order.
QudaGaugeParam param
Definition: pack_test.cpp:18
#define errorQuda(...)
Definition: util_quda.h:120