QUDA  v1.1.0
A library for QCD on GPUs
dslash_reference.cpp
Go to the documentation of this file.
1 // QUDA header (for HeavyQuarkResidualNorm)
2 #include <blas_quda.h>
3 
4 // External headers
5 #include <misc.h>
6 #include <host_utils.h>
7 #include <dslash_reference.h>
8 #include <wilson_dslash_reference.h>
9 #include <domain_wall_dslash_reference.h>
10 #include <staggered_dslash_reference.h>
11 #include <command_line_params.h>
12 
13 // Overload for workflows without multishift
14 void verifyInversion(void *spinorOut, void *spinorIn, void *spinorCheck, QudaGaugeParam &gauge_param,
15  QudaInvertParam &inv_param, void **gauge, void *clover, void *clover_inv)
16 {
17  void **spinorOutMulti = nullptr;
18  verifyInversion(spinorOut, spinorOutMulti, spinorIn, spinorCheck, gauge_param, inv_param, gauge, clover, clover_inv);
19 }
20 
21 void verifyInversion(void *spinorOut, void **spinorOutMulti, void *spinorIn, void *spinorCheck,
22  QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover,
23  void *clover_inv)
24 {
25 
26  if (dslash_type == QUDA_DOMAIN_WALL_DSLASH || dslash_type == QUDA_DOMAIN_WALL_4D_DSLASH
27  || dslash_type == QUDA_MOBIUS_DWF_DSLASH || dslash_type == QUDA_MOBIUS_DWF_EOFA_DSLASH) {
28  verifyDomainWallTypeInversion(spinorOut, spinorOutMulti, spinorIn, spinorCheck, gauge_param, inv_param, gauge,
29  clover, clover_inv);
30  } else if (dslash_type == QUDA_WILSON_DSLASH || dslash_type == QUDA_CLOVER_WILSON_DSLASH
31  || dslash_type == QUDA_TWISTED_MASS_DSLASH || dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
32  verifyWilsonTypeInversion(spinorOut, spinorOutMulti, spinorIn, spinorCheck, gauge_param, inv_param, gauge, clover,
33  clover_inv);
34  } else {
35  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
36  }
37 }
38 
39 void verifyDomainWallTypeInversion(void *spinorOut, void **spinorOutMulti, void *spinorIn, void *spinorCheck,
40  QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover,
41  void *clover_inv)
42 {
43  if (inv_param.solution_type == QUDA_MAT_SOLUTION) {
44  if (dslash_type == QUDA_DOMAIN_WALL_DSLASH) {
45  dw_mat(spinorCheck, gauge, spinorOut, kappa5, inv_param.dagger, inv_param.cpu_prec, gauge_param, inv_param.mass);
46  } else if (dslash_type == QUDA_DOMAIN_WALL_4D_DSLASH) {
47  dw_4d_mat(spinorCheck, gauge, spinorOut, kappa5, inv_param.dagger, inv_param.cpu_prec, gauge_param, inv_param.mass);
48  } else if (dslash_type == QUDA_MOBIUS_DWF_DSLASH) {
49  double _Complex *kappa_b = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
50  double _Complex *kappa_c = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
51  for (int xs = 0; xs < Lsdim; xs++) {
52  kappa_b[xs] = 1.0 / (2 * (inv_param.b_5[xs] * (4.0 + inv_param.m5) + 1.0));
53  kappa_c[xs] = 1.0 / (2 * (inv_param.c_5[xs] * (4.0 + inv_param.m5) - 1.0));
54  }
55  mdw_mat(spinorCheck, gauge, spinorOut, kappa_b, kappa_c, inv_param.dagger, inv_param.cpu_prec, gauge_param,
56  inv_param.mass, inv_param.b_5, inv_param.c_5);
57  free(kappa_b);
58  free(kappa_c);
59  } else if (dslash_type == QUDA_MOBIUS_DWF_EOFA_DSLASH) {
60  mdw_eofa_mat(spinorCheck, gauge, spinorOut, inv_param.dagger, inv_param.cpu_prec, gauge_param, inv_param.mass,
61  inv_param.m5, (__real__ inv_param.b_5[0]), (__real__ inv_param.c_5[0]), inv_param.mq1, inv_param.mq2,
62  inv_param.mq3, inv_param.eofa_pm, inv_param.eofa_shift);
63  } else {
64  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
65  }
66 
67  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
68  ax(0.5 / kappa5, spinorCheck, V * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
69  }
70 
71  } else if (inv_param.solution_type == QUDA_MATPC_SOLUTION) {
72 
73  // DOMAIN_WALL START
74  if (dslash_type == QUDA_DOMAIN_WALL_DSLASH) {
75  dw_matpc(spinorCheck, gauge, spinorOut, kappa5, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
76  inv_param.mass);
77  } else if (dslash_type == QUDA_DOMAIN_WALL_4D_DSLASH) {
78  dw_4d_matpc(spinorCheck, gauge, spinorOut, kappa5, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
79  inv_param.mass);
80  } else if (dslash_type == QUDA_MOBIUS_DWF_DSLASH) {
81  double _Complex *kappa_b = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
82  double _Complex *kappa_c = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
83  for (int xs = 0; xs < Lsdim; xs++) {
84  kappa_b[xs] = 1.0 / (2 * (inv_param.b_5[xs] * (4.0 + inv_param.m5) + 1.0));
85  kappa_c[xs] = 1.0 / (2 * (inv_param.c_5[xs] * (4.0 + inv_param.m5) - 1.0));
86  }
87  mdw_matpc(spinorCheck, gauge, spinorOut, kappa_b, kappa_c, inv_param.matpc_type, 0, inv_param.cpu_prec,
88  gauge_param, inv_param.mass, inv_param.b_5, inv_param.c_5);
89  free(kappa_b);
90  free(kappa_c);
91  // DOMAIN_WALL END
92  } else if (dslash_type == QUDA_MOBIUS_DWF_EOFA_DSLASH) {
93  mdw_eofa_matpc(spinorCheck, gauge, spinorOut, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
94  inv_param.mass, inv_param.m5, (__real__ inv_param.b_5[0]), (__real__ inv_param.c_5[0]),
95  inv_param.mq1, inv_param.mq2, inv_param.mq3, inv_param.eofa_pm, inv_param.eofa_shift);
96  } else {
97  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
98  }
99 
100  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
101  ax(0.25 / (kappa5 * kappa5), spinorCheck, V * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
102  }
103 
104  } else if (inv_param.solution_type == QUDA_MATPCDAG_MATPC_SOLUTION) {
105 
106  void *spinorTmp = malloc(V * spinor_site_size * host_spinor_data_type_size * inv_param.Ls);
107  ax(0, spinorCheck, V * spinor_site_size, inv_param.cpu_prec);
108 
109  // DOMAIN_WALL START
110  if (dslash_type == QUDA_DOMAIN_WALL_DSLASH) {
111  dw_matpc(spinorTmp, gauge, spinorOut, kappa5, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
112  inv_param.mass);
113  dw_matpc(spinorCheck, gauge, spinorTmp, kappa5, inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param,
114  inv_param.mass);
115  } else if (dslash_type == QUDA_DOMAIN_WALL_4D_DSLASH) {
116  dw_4d_matpc(spinorTmp, gauge, spinorOut, kappa5, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
117  inv_param.mass);
118  dw_4d_matpc(spinorCheck, gauge, spinorTmp, kappa5, inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param,
119  inv_param.mass);
120  } else if (dslash_type == QUDA_MOBIUS_DWF_DSLASH) {
121  double _Complex *kappa_b = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
122  double _Complex *kappa_c = (double _Complex *)malloc(Lsdim * sizeof(double _Complex));
123  for (int xs = 0; xs < Lsdim; xs++) {
124  kappa_b[xs] = 1.0 / (2 * (inv_param.b_5[xs] * (4.0 + inv_param.m5) + 1.0));
125  kappa_c[xs] = 1.0 / (2 * (inv_param.c_5[xs] * (4.0 + inv_param.m5) - 1.0));
126  }
127  mdw_matpc(spinorTmp, gauge, spinorOut, kappa_b, kappa_c, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
128  inv_param.mass, inv_param.b_5, inv_param.c_5);
129  mdw_matpc(spinorCheck, gauge, spinorTmp, kappa_b, kappa_c, inv_param.matpc_type, 1, inv_param.cpu_prec,
130  gauge_param, inv_param.mass, inv_param.b_5, inv_param.c_5);
131  free(kappa_b);
132  free(kappa_c);
133  // DOMAIN_WALL END
134  } else if (dslash_type == QUDA_MOBIUS_DWF_EOFA_DSLASH) {
135  mdw_eofa_matpc(spinorTmp, gauge, spinorOut, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param,
136  inv_param.mass, inv_param.m5, (__real__ inv_param.b_5[0]), (__real__ inv_param.c_5[0]),
137  inv_param.mq1, inv_param.mq2, inv_param.mq3, inv_param.eofa_pm, inv_param.eofa_shift);
138  mdw_eofa_matpc(spinorCheck, gauge, spinorTmp, inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param,
139  inv_param.mass, inv_param.m5, (__real__ inv_param.b_5[0]), (__real__ inv_param.c_5[0]),
140  inv_param.mq1, inv_param.mq2, inv_param.mq3, inv_param.eofa_pm, inv_param.eofa_shift);
141 
142  } else {
143  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
144  }
145 
146  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
147  errorQuda("Mass normalization %s not implemented", get_mass_normalization_str(inv_param.mass_normalization));
148  }
149 
150  free(spinorTmp);
151  } else {
152  errorQuda("Solution type %s not implemented", get_solution_str(inv_param.solution_type));
153  }
154 
155  int vol = inv_param.solution_type == QUDA_MAT_SOLUTION ? V : Vh;
156  mxpy(spinorIn, spinorCheck, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
157  double nrm2 = norm_2(spinorCheck, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
158  double src2 = norm_2(spinorIn, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
159  double l2r = sqrt(nrm2 / src2);
160 
161  printfQuda("Residuals: (L2 relative) tol %9.6e, QUDA = %9.6e, host = %9.6e; (heavy-quark) tol %9.6e, QUDA = %9.6e\n",
162  inv_param.tol, inv_param.true_res, l2r, inv_param.tol_hq, inv_param.true_res_hq);
163 }
164 
165 void verifyWilsonTypeInversion(void *spinorOut, void **spinorOutMulti, void *spinorIn, void *spinorCheck,
166  QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover,
167  void *clover_inv)
168 {
169  if (multishift > 1) {
170  // ONLY WILSON/CLOVER/TWISTED TYPES
171  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
172  errorQuda("Mass normalization %s not implemented", get_mass_normalization_str(inv_param.mass_normalization));
173  }
174 
175  void *spinorTmp = malloc(Vh * spinor_site_size * host_spinor_data_type_size * inv_param.Ls);
176  printfQuda("Host residuum checks: \n");
177  for (int i = 0; i < inv_param.num_offset; i++) {
178  ax(0, spinorCheck, Vh * spinor_site_size, inv_param.cpu_prec);
179 
180  if (dslash_type == QUDA_TWISTED_MASS_DSLASH) {
181  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET) {
182  int tm_offset = Vh * spinor_site_size;
183  void *out0 = spinorCheck;
184  void *out1 = (char *)out0 + tm_offset * cpu_prec;
185 
186  void *tmp0 = spinorTmp;
187  void *tmp1 = (char *)tmp0 + tm_offset * cpu_prec;
188 
189  void *in0 = spinorOutMulti[i];
190  void *in1 = (char *)in0 + tm_offset * cpu_prec;
191 
192  tm_ndeg_matpc(tmp0, tmp1, gauge, in0, in1, inv_param.kappa, inv_param.mu, inv_param.epsilon,
193  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
194  tm_ndeg_matpc(out0, out1, gauge, tmp0, tmp1, inv_param.kappa, inv_param.mu, inv_param.epsilon,
195  inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
196  } else {
197  tm_matpc(spinorTmp, gauge, spinorOutMulti[i], inv_param.kappa, inv_param.mu, inv_param.twist_flavor,
198  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
199  tm_matpc(spinorCheck, gauge, spinorTmp, inv_param.kappa, inv_param.mu, inv_param.twist_flavor,
200  inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
201  }
202  } else if (dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
203  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET)
204  errorQuda("Twisted mass solution type %s not supported", get_flavor_str(inv_param.twist_flavor));
205  tmc_matpc(spinorTmp, gauge, spinorOutMulti[i], clover, clover_inv, inv_param.kappa, inv_param.mu,
206  inv_param.twist_flavor, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
207  tmc_matpc(spinorCheck, gauge, spinorTmp, clover, clover_inv, inv_param.kappa, inv_param.mu,
208  inv_param.twist_flavor, inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
209  } else if (dslash_type == QUDA_WILSON_DSLASH) {
210  wil_matpc(spinorTmp, gauge, spinorOutMulti[i], inv_param.kappa, inv_param.matpc_type, 0, inv_param.cpu_prec,
211  gauge_param);
212  wil_matpc(spinorCheck, gauge, spinorTmp, inv_param.kappa, inv_param.matpc_type, 1, inv_param.cpu_prec,
213  gauge_param);
214  } else if (dslash_type == QUDA_CLOVER_WILSON_DSLASH) {
215  clover_matpc(spinorTmp, gauge, clover, clover_inv, spinorOutMulti[i], inv_param.kappa, inv_param.matpc_type, 0,
216  inv_param.cpu_prec, gauge_param);
217  clover_matpc(spinorCheck, gauge, clover, clover_inv, spinorTmp, inv_param.kappa, inv_param.matpc_type, 1,
218  inv_param.cpu_prec, gauge_param);
219  } else {
220  printfQuda("Domain wall not supported for multi-shift\n");
221  exit(-1);
222  }
223 
224  axpy(inv_param.offset[i], spinorOutMulti[i], spinorCheck, Vh * spinor_site_size, inv_param.cpu_prec);
225  mxpy(spinorIn, spinorCheck, Vh * spinor_site_size, inv_param.cpu_prec);
226  double nrm2 = norm_2(spinorCheck, Vh * spinor_site_size, inv_param.cpu_prec);
227  double src2 = norm_2(spinorIn, Vh * spinor_site_size, inv_param.cpu_prec);
228  double l2r = sqrt(nrm2 / src2);
229 
230  printfQuda("Shift %2d residuals: (L2 relative) tol %9.6e, QUDA = %9.6e, host = %9.6e; (heavy-quark) tol %9.6e, "
231  "QUDA = %9.6e\n",
232  i, inv_param.tol_offset[i], inv_param.true_res_offset[i], l2r, inv_param.tol_hq_offset[i],
233  inv_param.true_res_hq_offset[i]);
234  }
235  free(spinorTmp);
236 
237  } else {
238  // Non-multishift workflow
239  if (inv_param.solution_type == QUDA_MAT_SOLUTION) {
240  if (dslash_type == QUDA_TWISTED_MASS_DSLASH) {
241  if (inv_param.twist_flavor == QUDA_TWIST_SINGLET) {
242  tm_mat(spinorCheck, gauge, spinorOut, inv_param.kappa, inv_param.mu, inv_param.twist_flavor, 0,
243  inv_param.cpu_prec, gauge_param);
244  } else {
245  int tm_offset = V * spinor_site_size;
246  void *evenOut = spinorCheck;
247  void *oddOut = (char *)evenOut + tm_offset * cpu_prec;
248 
249  void *evenIn = spinorOut;
250  void *oddIn = (char *)evenIn + tm_offset * cpu_prec;
251 
252  tm_ndeg_mat(evenOut, oddOut, gauge, evenIn, oddIn, inv_param.kappa, inv_param.mu, inv_param.epsilon, 0,
253  inv_param.cpu_prec, gauge_param);
254  }
255  } else if (dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
256  tmc_mat(spinorCheck, gauge, clover, spinorOut, inv_param.kappa, inv_param.mu, inv_param.twist_flavor, 0,
257  inv_param.cpu_prec, gauge_param);
258  } else if (dslash_type == QUDA_WILSON_DSLASH) {
259  wil_mat(spinorCheck, gauge, spinorOut, inv_param.kappa, 0, inv_param.cpu_prec, gauge_param);
260  } else if (dslash_type == QUDA_CLOVER_WILSON_DSLASH) {
261  clover_mat(spinorCheck, gauge, clover, spinorOut, inv_param.kappa, 0, inv_param.cpu_prec, gauge_param);
262  } else {
263  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
264  }
265  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
266  if (dslash_type == QUDA_TWISTED_MASS_DSLASH && twist_flavor == QUDA_TWIST_NONDEG_DOUBLET) {
267  ax(0.5 / inv_param.kappa, spinorCheck, 2 * V * spinor_site_size, inv_param.cpu_prec);
268  // CAREFULL
269  } else {
270  ax(0.5 / inv_param.kappa, spinorCheck, V * spinor_site_size, inv_param.cpu_prec);
271  }
272  }
273 
274  } else if (inv_param.solution_type == QUDA_MATPC_SOLUTION) {
275 
276  if (dslash_type == QUDA_TWISTED_MASS_DSLASH) {
277  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET) {
278  int tm_offset = Vh * spinor_site_size;
279  void *out0 = spinorCheck;
280  void *out1 = (char *)out0 + tm_offset * cpu_prec;
281 
282  void *in0 = spinorOut;
283  void *in1 = (char *)in0 + tm_offset * cpu_prec;
284 
285  tm_ndeg_matpc(out0, out1, gauge, in0, in1, inv_param.kappa, inv_param.mu, inv_param.epsilon,
286  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
287  } else {
288  tm_matpc(spinorCheck, gauge, spinorOut, inv_param.kappa, inv_param.mu, inv_param.twist_flavor,
289  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
290  }
291  } else if (dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
292  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET)
293  errorQuda("Twisted mass solution type %s not supported", get_flavor_str(inv_param.twist_flavor));
294  tmc_matpc(spinorCheck, gauge, spinorOut, clover, clover_inv, inv_param.kappa, inv_param.mu,
295  inv_param.twist_flavor, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
296  } else if (dslash_type == QUDA_WILSON_DSLASH) {
297  wil_matpc(spinorCheck, gauge, spinorOut, inv_param.kappa, inv_param.matpc_type, 0, inv_param.cpu_prec,
298  gauge_param);
299  } else if (dslash_type == QUDA_CLOVER_WILSON_DSLASH) {
300  clover_matpc(spinorCheck, gauge, clover, clover_inv, spinorOut, inv_param.kappa, inv_param.matpc_type, 0,
301  inv_param.cpu_prec, gauge_param);
302  } else {
303  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
304  }
305 
306  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
307  ax(0.25 / (inv_param.kappa * inv_param.kappa), spinorCheck, Vh * spinor_site_size, inv_param.cpu_prec);
308  }
309 
310  } else if (inv_param.solution_type == QUDA_MATPCDAG_MATPC_SOLUTION) {
311 
312  void *spinorTmp = malloc(V * spinor_site_size * host_spinor_data_type_size * inv_param.Ls);
313  ax(0, spinorCheck, V * spinor_site_size, inv_param.cpu_prec);
314 
315  if (dslash_type == QUDA_TWISTED_MASS_DSLASH) {
316  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET) {
317  int tm_offset = Vh * spinor_site_size;
318  void *out0 = spinorCheck;
319  void *out1 = (char *)out0 + tm_offset * cpu_prec;
320 
321  void *tmp0 = spinorTmp;
322  void *tmp1 = (char *)tmp0 + tm_offset * cpu_prec;
323 
324  void *in0 = spinorOut;
325  void *in1 = (char *)in0 + tm_offset * cpu_prec;
326 
327  tm_ndeg_matpc(tmp0, tmp1, gauge, in0, in1, inv_param.kappa, inv_param.mu, inv_param.epsilon,
328  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
329  tm_ndeg_matpc(out0, out1, gauge, tmp0, tmp1, inv_param.kappa, inv_param.mu, inv_param.epsilon,
330  inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
331  } else {
332  tm_matpc(spinorTmp, gauge, spinorOut, inv_param.kappa, inv_param.mu, inv_param.twist_flavor,
333  inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
334  tm_matpc(spinorCheck, gauge, spinorTmp, inv_param.kappa, inv_param.mu, inv_param.twist_flavor,
335  inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
336  }
337  } else if (dslash_type == QUDA_TWISTED_CLOVER_DSLASH) {
338  if (inv_param.twist_flavor != QUDA_TWIST_SINGLET)
339  errorQuda("Twisted mass solution type %s not supported", get_flavor_str(inv_param.twist_flavor));
340  tmc_matpc(spinorTmp, gauge, spinorOut, clover, clover_inv, inv_param.kappa, inv_param.mu,
341  inv_param.twist_flavor, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
342  tmc_matpc(spinorCheck, gauge, spinorTmp, clover, clover_inv, inv_param.kappa, inv_param.mu,
343  inv_param.twist_flavor, inv_param.matpc_type, 1, inv_param.cpu_prec, gauge_param);
344  } else if (dslash_type == QUDA_WILSON_DSLASH) {
345  wil_matpc(spinorTmp, gauge, spinorOut, inv_param.kappa, inv_param.matpc_type, 0, inv_param.cpu_prec, gauge_param);
346  wil_matpc(spinorCheck, gauge, spinorTmp, inv_param.kappa, inv_param.matpc_type, 1, inv_param.cpu_prec,
347  gauge_param);
348  } else if (dslash_type == QUDA_CLOVER_WILSON_DSLASH) {
349  clover_matpc(spinorTmp, gauge, clover, clover_inv, spinorOut, inv_param.kappa, inv_param.matpc_type, 0,
350  inv_param.cpu_prec, gauge_param);
351  clover_matpc(spinorCheck, gauge, clover, clover_inv, spinorTmp, inv_param.kappa, inv_param.matpc_type, 1,
352  inv_param.cpu_prec, gauge_param);
353  } else {
354  errorQuda("Unsupported dslash_type=%s", get_dslash_str(dslash_type));
355  }
356 
357  if (inv_param.mass_normalization == QUDA_MASS_NORMALIZATION) {
358  errorQuda("Mass normalization %s not implemented", get_mass_normalization_str(inv_param.mass_normalization));
359  }
360 
361  free(spinorTmp);
362  } else {
363  errorQuda("Solution type %s not implemented", get_solution_str(inv_param.solution_type));
364  }
365 
366  int vol = inv_param.solution_type == QUDA_MAT_SOLUTION ? V : Vh;
367  mxpy(spinorIn, spinorCheck, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
368  double nrm2 = norm_2(spinorCheck, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
369  double src2 = norm_2(spinorIn, vol * spinor_site_size * inv_param.Ls, inv_param.cpu_prec);
370  double l2r = sqrt(nrm2 / src2);
371 
372  printfQuda(
373  "Residuals: (L2 relative) tol %9.6e, QUDA = %9.6e, host = %9.6e; (heavy-quark) tol %9.6e, QUDA = %9.6e\n",
374  inv_param.tol, inv_param.true_res, l2r, inv_param.tol_hq, inv_param.true_res_hq);
375  }
376 }
377 
378 void verifyStaggeredInversion(quda::ColorSpinorField *tmp, quda::ColorSpinorField *ref, quda::ColorSpinorField *in,
379  quda::ColorSpinorField *out, double mass, void *qdp_fatlink[], void *qdp_longlink[],
380  void **ghost_fatlink, void **ghost_longlink, QudaGaugeParam &gauge_param,
381  QudaInvertParam &inv_param, int shift)
382 {
383 
384  switch (test_type) {
385  case 0: // full parity solution, full parity system
386  case 1: // full parity solution, solving EVEN EVEN prec system
387  case 2: // full parity solution, solving ODD ODD prec system
388 
389  // In QUDA, the full staggered operator has the sign convention
390  // {{m, -D_eo},{-D_oe,m}}, while the CPU verify function does not
391  // have the minus sign. Passing in QUDA_DAG_YES solves this
392  // discrepancy.
393  staggeredDslash(reinterpret_cast<quda::cpuColorSpinorField *>(&ref->Even()), qdp_fatlink, qdp_longlink,
394  ghost_fatlink, ghost_longlink, reinterpret_cast<quda::cpuColorSpinorField *>(&out->Odd()),
395  QUDA_EVEN_PARITY, QUDA_DAG_YES, inv_param.cpu_prec, gauge_param.cpu_prec, dslash_type);
396  staggeredDslash(reinterpret_cast<quda::cpuColorSpinorField *>(&ref->Odd()), qdp_fatlink, qdp_longlink,
397  ghost_fatlink, ghost_longlink, reinterpret_cast<quda::cpuColorSpinorField *>(&out->Even()),
398  QUDA_ODD_PARITY, QUDA_DAG_YES, inv_param.cpu_prec, gauge_param.cpu_prec, dslash_type);
399 
400  if (dslash_type == QUDA_LAPLACE_DSLASH) {
401  xpay(out->V(), kappa, ref->V(), ref->Length(), gauge_param.cpu_prec);
402  ax(0.5 / kappa, ref->V(), ref->Length(), gauge_param.cpu_prec);
403  } else {
404  axpy(2 * mass, out->V(), ref->V(), ref->Length(), gauge_param.cpu_prec);
405  }
406  break;
407 
408  case 3: // even parity solution, solving EVEN system
409  case 4: // odd parity solution, solving ODD system
410  case 5: // multi mass CG, even parity solution, solving EVEN system
411  case 6: // multi mass CG, odd parity solution, solving ODD system
412 
413  staggeredMatDagMat(ref, qdp_fatlink, qdp_longlink, ghost_fatlink, ghost_longlink, out, mass, 0, inv_param.cpu_prec,
414  gauge_param.cpu_prec, tmp,
415  (test_type == 3 || test_type == 5) ? QUDA_EVEN_PARITY : QUDA_ODD_PARITY, dslash_type);
416  break;
417  }
418 
419  int len = 0;
420  if (solution_type == QUDA_MAT_SOLUTION || solution_type == QUDA_MATDAG_MAT_SOLUTION) {
421  len = V;
422  } else {
423  len = Vh;
424  }
425 
426  mxpy(in->V(), ref->V(), len * stag_spinor_site_size, inv_param.cpu_prec);
427  double nrm2 = norm_2(ref->V(), len * stag_spinor_site_size, inv_param.cpu_prec);
428  double src2 = norm_2(in->V(), len * stag_spinor_site_size, inv_param.cpu_prec);
429  double hqr = sqrt(quda::blas::HeavyQuarkResidualNorm(*out, *ref).z);
430  double l2r = sqrt(nrm2 / src2);
431 
432  if (multishift == 1) {
433  printfQuda("Residuals: (L2 relative) tol %9.6e, QUDA = %9.6e, host = %9.6e; (heavy-quark) tol %9.6e, QUDA = %9.6e, "
434  "host = %9.6e\n",
435  inv_param.tol, inv_param.true_res, l2r, inv_param.tol_hq, inv_param.true_res_hq, hqr);
436  } else {
437  printfQuda("Shift %2d residuals: (L2 relative) tol %9.6e, QUDA = %9.6e, host = %9.6e; (heavy-quark) tol %9.6e, "
438  "QUDA = %9.6e, host = %9.6e\n",
439  shift, inv_param.tol_offset[shift], inv_param.true_res_offset[shift], l2r,
440  inv_param.tol_hq_offset[shift], inv_param.true_res_hq_offset[shift], hqr);
441  // Empirical: if the cpu residue is more than 1 order the target accuracy, then it fails to converge
442  if (sqrt(nrm2 / src2) > 10 * inv_param.tol_offset[shift]) {
443  printfQuda("Shift %2d has empirically failed to converge\n", shift);
444  }
445  }
446 }
blas_quda.h
quda::ColorSpinorField::Odd
const ColorSpinorField & Odd() const
Definition: color_spinor_field.cpp:578
quda::ColorSpinorField::Even
const ColorSpinorField & Even() const
Definition: color_spinor_field.cpp:570
clover_mat
void clover_mat(void *out, void **gauge, void *clover, void *in, double kappa, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: clover_reference.cpp:147
tmc_mat
void tmc_mat(void *out, void **gauge, void *clover, void *in, double kappa, double mu, QudaTwistFlavorType flavor, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: clover_reference.cpp:267
clover_matpc
void clover_matpc(void *out, void **gauge, void *clover, void *clover_inv, void *in, double kappa, QudaMatPCType matpc_type, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: clover_reference.cpp:90
tmc_matpc
void tmc_matpc(void *out, void **gauge, void *in, void *clover, void *cInv, double kappa, double mu, QudaTwistFlavorType flavor, QudaMatPCType matpc_type, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: clover_reference.cpp:294
kappa
double kappa
Definition: command_line_params.cpp:72
mass
double mass
Definition: command_line_params.cpp:71
test_type
int test_type
Definition: command_line_params.cpp:57
twist_flavor
QudaTwistFlavorType twist_flavor
Definition: command_line_params.cpp:91
solution_type
QudaSolutionType solution_type
Definition: command_line_params.cpp:95
dslash_type
QudaDslashType dslash_type
Definition: command_line_params.cpp:41
multishift
int multishift
Definition: command_line_params.cpp:67
Lsdim
int Lsdim
Definition: command_line_params.cpp:39
command_line_params.h
Vh
int Vh
Definition: host_utils.cpp:38
V
int V
Definition: host_utils.cpp:37
tmp
cudaColorSpinorField * tmp
Definition: covdev_test.cpp:34
spinorOut
cpuColorSpinorField * spinorOut
Definition: covdev_test.cpp:31
gauge_param
QudaGaugeParam gauge_param
Definition: covdev_test.cpp:26
inv_param
QudaInvertParam inv_param
Definition: covdev_test.cpp:27
dw_mat
void dw_mat(void *out, void **gauge, void *in, double kappa, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm)
Definition: domain_wall_dslash_reference.cpp:991
dw_4d_matpc
void dw_4d_matpc(void *out, void **gauge, void *in, double kappa, QudaMatPCType matpc_type, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm)
Definition: domain_wall_dslash_reference.cpp:1154
mdw_matpc
void mdw_matpc(void *out, void **gauge, void *in, double _Complex *kappa_b, double _Complex *kappa_c, QudaMatPCType matpc_type, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm, double _Complex *b5, double _Complex *c5)
Definition: domain_wall_dslash_reference.cpp:1194
dw_4d_mat
void dw_4d_mat(void *out, void **gauge, void *in, double kappa, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm)
Definition: domain_wall_dslash_reference.cpp:1005
mdw_mat
void mdw_mat(void *out, void **gauge, void *in, double _Complex *kappa_b, double _Complex *kappa_c, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm, double _Complex *b5, double _Complex *c5)
Definition: domain_wall_dslash_reference.cpp:1022
mdw_eofa_matpc
void mdw_eofa_matpc(void *out, void **gauge, void *in, QudaMatPCType matpc_type, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm, double m5, double b, double c, double mq1, double mq2, double mq3, int eofa_pm, double eofa_shift)
Definition: domain_wall_dslash_reference.cpp:1267
dw_matpc
void dw_matpc(void *out, void **gauge, void *in, double kappa, QudaMatPCType matpc_type, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm)
Definition: domain_wall_dslash_reference.cpp:1134
mdw_eofa_mat
void mdw_eofa_mat(void *out, void **gauge, void *in, int dagger, QudaPrecision precision, QudaGaugeParam &gauge_param, double mferm, double m5, double b, double c, double mq1, double mq2, double mq3, int eofa_pm, double eofa_shift)
Definition: domain_wall_dslash_reference.cpp:1069
domain_wall_dslash_reference.h
verifyInversion
void verifyInversion(void *spinorOut, void *spinorIn, void *spinorCheck, QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover, void *clover_inv)
Definition: dslash_reference.cpp:14
verifyStaggeredInversion
void verifyStaggeredInversion(quda::ColorSpinorField *tmp, quda::ColorSpinorField *ref, quda::ColorSpinorField *in, quda::ColorSpinorField *out, double mass, void *qdp_fatlink[], void *qdp_longlink[], void **ghost_fatlink, void **ghost_longlink, QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, int shift)
Definition: dslash_reference.cpp:378
verifyDomainWallTypeInversion
void verifyDomainWallTypeInversion(void *spinorOut, void **spinorOutMulti, void *spinorIn, void *spinorCheck, QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover, void *clover_inv)
Definition: dslash_reference.cpp:39
verifyWilsonTypeInversion
void verifyWilsonTypeInversion(void *spinorOut, void **spinorOutMulti, void *spinorIn, void *spinorCheck, QudaGaugeParam &gauge_param, QudaInvertParam &inv_param, void **gauge, void *clover, void *clover_inv)
Definition: dslash_reference.cpp:165
dslash_reference.h
QUDA_WILSON_DSLASH
@ QUDA_WILSON_DSLASH
Definition: enum_quda.h:90
QUDA_TWISTED_CLOVER_DSLASH
@ QUDA_TWISTED_CLOVER_DSLASH
Definition: enum_quda.h:100
QUDA_MOBIUS_DWF_DSLASH
@ QUDA_MOBIUS_DWF_DSLASH
Definition: enum_quda.h:95
QUDA_CLOVER_WILSON_DSLASH
@ QUDA_CLOVER_WILSON_DSLASH
Definition: enum_quda.h:91
QUDA_TWISTED_MASS_DSLASH
@ QUDA_TWISTED_MASS_DSLASH
Definition: enum_quda.h:99
QUDA_DOMAIN_WALL_DSLASH
@ QUDA_DOMAIN_WALL_DSLASH
Definition: enum_quda.h:93
QUDA_MOBIUS_DWF_EOFA_DSLASH
@ QUDA_MOBIUS_DWF_EOFA_DSLASH
Definition: enum_quda.h:96
QUDA_LAPLACE_DSLASH
@ QUDA_LAPLACE_DSLASH
Definition: enum_quda.h:101
QUDA_DOMAIN_WALL_4D_DSLASH
@ QUDA_DOMAIN_WALL_4D_DSLASH
Definition: enum_quda.h:94
QUDA_MASS_NORMALIZATION
@ QUDA_MASS_NORMALIZATION
Definition: enum_quda.h:227
QUDA_DAG_YES
@ QUDA_DAG_YES
Definition: enum_quda.h:223
QUDA_EVEN_PARITY
@ QUDA_EVEN_PARITY
Definition: enum_quda.h:284
QUDA_ODD_PARITY
@ QUDA_ODD_PARITY
Definition: enum_quda.h:284
QUDA_MATPC_SOLUTION
@ QUDA_MATPC_SOLUTION
Definition: enum_quda.h:159
QUDA_MATDAG_MAT_SOLUTION
@ QUDA_MATDAG_MAT_SOLUTION
Definition: enum_quda.h:158
QUDA_MATPCDAG_MATPC_SOLUTION
@ QUDA_MATPCDAG_MATPC_SOLUTION
Definition: enum_quda.h:161
QUDA_MAT_SOLUTION
@ QUDA_MAT_SOLUTION
Definition: enum_quda.h:157
QUDA_TWIST_SINGLET
@ QUDA_TWIST_SINGLET
Definition: enum_quda.h:400
QUDA_TWIST_NONDEG_DOUBLET
@ QUDA_TWIST_NONDEG_DOUBLET
Definition: enum_quda.h:401
norm_2
double norm_2(void *v, int len, QudaPrecision precision)
Definition: host_blas.cpp:48
host_spinor_data_type_size
size_t host_spinor_data_type_size
Definition: host_utils.cpp:66
cpu_prec
QudaPrecision & cpu_prec
Definition: host_utils.cpp:57
kappa5
double kappa5
Definition: host_utils.cpp:51
host_utils.h
stag_spinor_site_size
#define stag_spinor_site_size
Definition: host_utils.h:10
spinor_site_size
#define spinor_site_size
Definition: host_utils.h:9
get_solution_str
const char * get_solution_str(QudaSolutionType type)
Definition: misc.cpp:215
get_flavor_str
const char * get_flavor_str(QudaTwistFlavorType type)
Definition: misc.cpp:248
get_dslash_str
const char * get_dslash_str(QudaDslashType type)
Definition: misc.cpp:118
get_mass_normalization_str
const char * get_mass_normalization_str(QudaMassNormalization type)
Definition: misc.cpp:186
quda::blas::HeavyQuarkResidualNorm
double3 HeavyQuarkResidualNorm(ColorSpinorField &x, ColorSpinorField &r)
quda::blas::xpay
void xpay(ColorSpinorField &x, double a, ColorSpinorField &y)
Definition: blas_quda.h:45
quda::blas::ax
void ax(double a, ColorSpinorField &x)
quda::blas::axpy
void axpy(double a, ColorSpinorField &x, ColorSpinorField &y)
Definition: blas_quda.h:43
quda::blas::mxpy
void mxpy(ColorSpinorField &x, ColorSpinorField &y)
Definition: blas_quda.h:42
quda::sqrt
__host__ __device__ ValueType sqrt(ValueType x)
Definition: complex_quda.h:120
staggeredMatDagMat
void staggeredMatDagMat(ColorSpinorField *out, void **fatlink, void **longlink, void **ghost_fatlink, void **ghost_longlink, ColorSpinorField *in, double mass, int dagger_bit, QudaPrecision sPrecision, QudaPrecision gPrecision, ColorSpinorField *tmp, QudaParity parity, QudaDslashType dslash_type)
Definition: staggered_dslash_reference.cpp:170
staggeredDslash
void staggeredDslash(ColorSpinorField *out, void **fatlink, void **longlink, void **ghost_fatlink, void **ghost_longlink, ColorSpinorField *in, int oddBit, int daggerBit, QudaPrecision sPrecision, QudaPrecision gPrecision, QudaDslashType dslash_type)
Definition: staggered_dslash_reference.cpp:126
staggered_dslash_reference.h
QudaGaugeParam_s
Definition: quda.h:31
QudaGaugeParam_s::cpu_prec
QudaPrecision cpu_prec
Definition: quda.h:46
QudaInvertParam_s
Definition: quda.h:98
QudaInvertParam_s::solution_type
QudaSolutionType solution_type
Definition: quda.h:228
QudaInvertParam_s::tol_hq
double tol_hq
Definition: quda.h:140
QudaInvertParam_s::mass_normalization
QudaMassNormalization mass_normalization
Definition: quda.h:232
QudaInvertParam_s::tol_hq_offset
double tol_hq_offset[QUDA_MAX_MULTI_SHIFT]
Definition: quda.h:206
QudaInvertParam_s::mq2
double mq2
Definition: quda.h:128
QudaInvertParam_s::true_res
double true_res
Definition: quda.h:143
QudaInvertParam_s::num_offset
int num_offset
Definition: quda.h:186
QudaInvertParam_s::mass
double mass
Definition: quda.h:109
QudaInvertParam_s::m5
double m5
Definition: quda.h:112
QudaInvertParam_s::twist_flavor
QudaTwistFlavorType twist_flavor
Definition: quda.h:134
QudaInvertParam_s::eofa_pm
int eofa_pm
Definition: quda.h:126
QudaInvertParam_s::mq1
double mq1
Definition: quda.h:127
QudaInvertParam_s::mu
double mu
Definition: quda.h:131
QudaInvertParam_s::b_5
double_complex b_5[QUDA_MAX_DWF_LS]
Definition: quda.h:115
QudaInvertParam_s::true_res_hq
double true_res_hq
Definition: quda.h:144
QudaInvertParam_s::eofa_shift
double eofa_shift
Definition: quda.h:125
QudaInvertParam_s::offset
double offset[QUDA_MAX_MULTI_SHIFT]
Definition: quda.h:200
QudaInvertParam_s::dagger
QudaDagType dagger
Definition: quda.h:231
QudaInvertParam_s::Ls
int Ls
Definition: quda.h:113
QudaInvertParam_s::epsilon
double epsilon
Definition: quda.h:132
QudaInvertParam_s::cpu_prec
QudaPrecision cpu_prec
Definition: quda.h:237
QudaInvertParam_s::true_res_offset
double true_res_offset[QUDA_MAX_MULTI_SHIFT]
Definition: quda.h:209
QudaInvertParam_s::matpc_type
QudaMatPCType matpc_type
Definition: quda.h:230
QudaInvertParam_s::c_5
double_complex c_5[QUDA_MAX_DWF_LS]
Definition: quda.h:116
QudaInvertParam_s::tol_offset
double tol_offset[QUDA_MAX_MULTI_SHIFT]
Definition: quda.h:203
QudaInvertParam_s::mq3
double mq3
Definition: quda.h:129
QudaInvertParam_s::tol
double tol
Definition: quda.h:138
QudaInvertParam_s::true_res_hq_offset
double true_res_hq_offset[QUDA_MAX_MULTI_SHIFT]
Definition: quda.h:215
QudaInvertParam_s::kappa
double kappa
Definition: quda.h:110
printfQuda
#define printfQuda(...)
Definition: util_quda.h:114
errorQuda
#define errorQuda(...)
Definition: util_quda.h:120
wil_mat
void wil_mat(void *out, void **gauge, void *in, double kappa, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:294
tm_mat
void tm_mat(void *out, void **gauge, void *in, double kappa, double mu, QudaTwistFlavorType flavor, int dagger_bit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:309
tm_ndeg_mat
void tm_ndeg_mat(void *evenOut, void *oddOut, void **gauge, void *evenIn, void *oddIn, double kappa, double mu, double epsilon, int daggerBit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:554
wil_matpc
void wil_matpc(void *outEven, void **gauge, void *inEven, double kappa, QudaMatPCType matpc_type, int daggerBit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:332
tm_ndeg_matpc
void tm_ndeg_matpc(void *outEven1, void *outEven2, void **gauge, void *inEven1, void *inEven2, double kappa, double mu, double epsilon, QudaMatPCType matpc_type, int daggerBit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:484
tm_matpc
void tm_matpc(void *outEven, void **gauge, void *inEven, double kappa, double mu, QudaTwistFlavorType flavor, QudaMatPCType matpc_type, int daggerBit, QudaPrecision precision, QudaGaugeParam &gauge_param)
Definition: wilson_dslash_reference.cpp:356
wilson_dslash_reference.h
Generated on Thu Oct 28 2021 16:10:27 for QUDA by doxygen 1.9.1