QUDA  v1.1.0
A library for QCD on GPUs
covdev_reference.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <string.h>
5 
6 #include <host_utils.h>
7 #include <misc.h>
8 #include <covdev_reference.h>
9 #include <dslash_reference.h>
10 
11 #include <quda_internal.h>
12 #include <quda.h>
13 #include <util_quda.h>
14 #include <blas_quda.h>
15 
16 extern void *memset(void *s, int c, size_t n);
17 
18 // covdevReference()
19 //
20 // if oddBit is zero: calculate even parity spinor elements (using odd parity spinor)
21 // if oddBit is one: calculate odd parity spinor elements
22 //
23 // if daggerBit is zero: perform ordinary covariant derivative operator
24 // if daggerBit is one: perform hermitian covariant derivative operator
25 //
26 template<typename Float>
28 {
29  int i, j;
30 
31  for (i = 0;i < 3; i++){
32  for(j=0;j < 3; j++){
33  printf("(%10f,%10f) \t", link[i*3*2 + j*2], link[i*3*2 + j*2 + 1]);
34  }
35  printf("\n");
36  }
37  printf("\n");
38  return;
39 }
40 
41 
42 template <typename sFloat, typename gFloat>
43 void covdevReference(sFloat *res, gFloat **link, sFloat *spinorField,
44  int oddBit, int daggerBit, int mu)
45 {
46  for (int i = 0; i < Vh * spinor_site_size; i++) res[i] = 0.0;
47 
48  gFloat *linkEven[4], *linkOdd[4];
49 
50  for (int dir = 0; dir < 4; dir++) {
51  linkEven[dir] = link[dir];
52  linkOdd[dir] = link[dir] + Vh * gauge_site_size;
53  }
54 
55  for (int sid = 0; sid < Vh; sid++) {
56  int offset = spinor_site_size * sid;
57 
58  sFloat gaugedSpinor[spinor_site_size];
59 
60  gFloat *lnk = gaugeLink(sid, mu, oddBit, linkEven, linkOdd, 1);
61  sFloat *spinor = spinorNeighbor(sid, mu, oddBit, spinorField, 1);
62 
63  if (daggerBit) {
64  for (int s = 0; s < 4; s++)
65  su3Tmul(&gaugedSpinor[s*6], lnk, &spinor[s*6]);
66  } else {
67  for (int s = 0; s < 4; s++)
68  su3Mul (&gaugedSpinor[s*6], lnk, &spinor[s*6]);
69  }
70 
71  sum(&res[offset], &res[offset], gaugedSpinor, spinor_site_size);
72  } // 4-d volume
73 }
74 
75 void covdev_dslash(void *res, void **link, void *spinorField, int oddBit, int daggerBit, int mu,
76  QudaPrecision sPrecision, QudaPrecision gPrecision)
77 {
78 
79  if (sPrecision == QUDA_DOUBLE_PRECISION) {
80  if (gPrecision == QUDA_DOUBLE_PRECISION){
81  covdevReference((double*)res, (double**)link, (double*)spinorField, oddBit, daggerBit, mu);
82  } else {
83  covdevReference((double*)res, (float**) link, (double*)spinorField, oddBit, daggerBit, mu);
84  }
85  }
86  else {
87  if (gPrecision == QUDA_DOUBLE_PRECISION){
88  covdevReference((float*)res, (double**)link, (float*)spinorField, oddBit, daggerBit, mu);
89  } else {
90  covdevReference((float*)res, (float**) link, (float*)spinorField, oddBit, daggerBit, mu);
91  }
92  }
93 }
94 
95 template <typename sFloat, typename gFloat>
96 void Mat(sFloat *out, gFloat **link, sFloat *in, int daggerBit, int mu)
97 {
98  sFloat *inEven = in;
99  sFloat *inOdd = in + Vh * spinor_site_size;
100  sFloat *outEven = out;
101  sFloat *outOdd = out + Vh * spinor_site_size;
102 
103  // full dslash operator
104  covdevReference(outOdd, link, inEven, 1, daggerBit, mu);
105  covdevReference(outEven, link, inOdd, 0, daggerBit, mu);
106 }
107 
108 
109 void mat(void *out, void **link, void *in, int dagger_bit, int mu,
110  QudaPrecision sPrecision, QudaPrecision gPrecision)
111 {
112 
113  if (sPrecision == QUDA_DOUBLE_PRECISION){
114  if (gPrecision == QUDA_DOUBLE_PRECISION) {
115  Mat((double*)out, (double**)link, (double*)in, dagger_bit, mu);
116  } else {
117  Mat((double*)out, (float**) link, (double*)in, dagger_bit, mu);
118  }
119  } else {
120  if (gPrecision == QUDA_DOUBLE_PRECISION){
121  Mat((float*)out, (double**)link, (float*)in, dagger_bit, mu);
122  } else {
123  Mat((float*)out, (float**) link, (float*)in, dagger_bit, mu);
124  }
125  }
126 }
127 
128 
129 
130 template <typename sFloat, typename gFloat>
131 void Matdagmat(sFloat *out, gFloat **link, sFloat *in, int daggerBit, int mu, sFloat* tmp, QudaParity parity)
132 {
133  switch(parity){
134  case QUDA_EVEN_PARITY:
135  {
136  sFloat *inEven = in;
137  sFloat *outEven = out;
138  covdevReference(tmp, link, inEven, 1, daggerBit, mu);
139  covdevReference(outEven, link, tmp, 0, daggerBit, mu);
140  break;
141  }
142  case QUDA_ODD_PARITY:
143  {
144  sFloat *inOdd = in;
145  sFloat *outOdd = out;
146  covdevReference(tmp, link, inOdd, 0, daggerBit, mu);
147  covdevReference(outOdd, link, tmp, 1, daggerBit, mu);
148  break;
149  }
150 
151  default:
152  fprintf(stderr, "ERROR: invalid parity in %s,line %d\n", __FUNCTION__, __LINE__);
153  break;
154  }
155 
156 }
157 
158 
159 
160 void matdagmat(void *out, void **link, void *in, int dagger_bit, int mu,
161  QudaPrecision sPrecision, QudaPrecision gPrecision, void *tmp, QudaParity parity)
162 {
163  if (sPrecision == QUDA_DOUBLE_PRECISION) {
164  if (gPrecision == QUDA_DOUBLE_PRECISION) {
165  Matdagmat((double*)out, (double**)link, (double*)in, dagger_bit, mu, (double*)tmp, parity);
166  } else {
167  Matdagmat((double*)out, (float**) link, (double*)in, dagger_bit, mu, (double*)tmp, parity);
168  }
169  } else {
170  if (gPrecision == QUDA_DOUBLE_PRECISION){
171  Matdagmat((float*)out, (double**)link, (float*)in, dagger_bit, mu, (float*)tmp, parity);
172  } else {
173  Matdagmat((float*)out, (float**) link, (float*)in, dagger_bit, mu, (float*)tmp, parity);
174  }
175  }
176 }
177 
178 #ifdef MULTI_GPU
179 
180 template <typename sFloat, typename gFloat>
181 void covdevReference_mg4dir(sFloat *res, gFloat **link, gFloat **ghostLink, sFloat *spinorField,
182  sFloat **fwd_nbr_spinor, sFloat **back_nbr_spinor, int oddBit, int daggerBit, int mu)
183 {
184  for (int i = 0; i < Vh * spinor_site_size; i++) res[i] = 0.0;
185 
186  gFloat *linkEven[4], *linkOdd[4];
187  gFloat *ghostLinkEven[4], *ghostLinkOdd[4];
188 
189  for (int dir = 0; dir < 4; dir++) {
190  linkEven[dir] = link[dir];
191  linkOdd[dir] = link[dir] + Vh * gauge_site_size;
192 
193  ghostLinkEven[dir] = ghostLink[dir];
194  ghostLinkOdd[dir] = ghostLink[dir] + (faceVolume[dir] / 2) * gauge_site_size;
195  }
196 
197  for (int sid = 0; sid < Vh; sid++) {
198  int offset = spinor_site_size * sid;
199 
200  gFloat *lnk = gaugeLink_mg4dir(sid, mu, oddBit, linkEven, linkOdd, ghostLinkEven, ghostLinkOdd, 1, 1);
201  sFloat *spinor = spinorNeighbor_mg4dir(sid, mu, oddBit, spinorField, fwd_nbr_spinor, back_nbr_spinor, 1, 1);
202 
203  sFloat gaugedSpinor[spinor_site_size];
204 
205  if (daggerBit) {
206  for (int s = 0; s < 4; s++)
207  su3Tmul(&gaugedSpinor[s*6], lnk, &spinor[s*6]);
208  } else {
209  for (int s = 0; s < 4; s++)
210  su3Mul (&gaugedSpinor[s*6], lnk, &spinor[s*6]);
211  }
212  sum(&res[offset], &res[offset], gaugedSpinor, spinor_site_size);
213  } // 4-d volume
214 }
215 
216 void covdev_dslash_mg4dir(cpuColorSpinorField* out, void **link, void** ghostLink,
217  cpuColorSpinorField* in, int oddBit, int daggerBit, int mu,
218  QudaPrecision sPrecision, QudaPrecision gPrecision)
219 {
220  QudaParity otherparity = QUDA_INVALID_PARITY;
221  if (oddBit == QUDA_EVEN_PARITY) {
222  otherparity = QUDA_ODD_PARITY;
223  } else if (oddBit == QUDA_ODD_PARITY) {
224  otherparity = QUDA_EVEN_PARITY;
225  } else {
226  errorQuda("ERROR: full parity not supported in function %s", __FUNCTION__);
227  }
228  const int nFace = 1;
229 
230  in->exchangeGhost(otherparity, nFace, daggerBit);
231 
232  void** fwd_nbr_spinor = in->fwdGhostFaceBuffer;
233  void** back_nbr_spinor = in->backGhostFaceBuffer;
234 
235  if (sPrecision == QUDA_DOUBLE_PRECISION) {
236  if (gPrecision == QUDA_DOUBLE_PRECISION) {
237  covdevReference_mg4dir((double*)out->V(), (double**)link, (double**)ghostLink, (double*)in->V(),
238  (double**)fwd_nbr_spinor, (double**)back_nbr_spinor, oddBit, daggerBit, mu);
239  } else {
240  covdevReference_mg4dir((double*)out->V(), (float**) link, (float**) ghostLink, (double*)in->V(),
241  (double**)fwd_nbr_spinor, (double**)back_nbr_spinor, oddBit, daggerBit, mu);
242  }
243  } else {
244  if (gPrecision == QUDA_DOUBLE_PRECISION) {
245  covdevReference_mg4dir((float*)out->V(), (double**)link, (double**)ghostLink, (float*)in->V(),
246  (float**)fwd_nbr_spinor, (float**)back_nbr_spinor, oddBit, daggerBit, mu);
247  } else {
248  covdevReference_mg4dir((float*)out->V(), (float**)link, (float**)ghostLink, (float*)in->V(),
249  (float**)fwd_nbr_spinor, (float**)back_nbr_spinor, oddBit, daggerBit, mu);
250  }
251  }
252 
253 }
254 
255 template <typename sFloat, typename gFloat>
256 void Mat_mg4dir(cpuColorSpinorField *out, gFloat **link, gFloat **ghostLink, cpuColorSpinorField *in, int daggerBit, int mu)
257 {
258  const int nFace = 1;
259  {
260  cpuColorSpinorField &inEven = static_cast<cpuColorSpinorField&>(in->Even());
261  cpuColorSpinorField &outOdd = static_cast<cpuColorSpinorField&>(out->Odd());
262 
263  inEven.exchangeGhost(QUDA_EVEN_PARITY, nFace, daggerBit);
264 
265  covdevReference_mg4dir(reinterpret_cast<sFloat*>(outOdd.V()), link, ghostLink,
266  reinterpret_cast<sFloat*>(inEven.V()),
267  reinterpret_cast<sFloat**>(inEven.fwdGhostFaceBuffer),
268  reinterpret_cast<sFloat**>(inEven.backGhostFaceBuffer),
269  1, daggerBit, mu);
270  }
271 
272  {
273  cpuColorSpinorField &inOdd = static_cast<cpuColorSpinorField&>(in->Odd());
274  cpuColorSpinorField &outEven = static_cast<cpuColorSpinorField&>(out->Even());
275 
276  inOdd.exchangeGhost(QUDA_ODD_PARITY, nFace, daggerBit);
277 
278  covdevReference_mg4dir(reinterpret_cast<sFloat*>(outEven.V()), link, ghostLink,
279  reinterpret_cast<sFloat*>(inOdd.V()),
280  reinterpret_cast<sFloat**>(inOdd.fwdGhostFaceBuffer),
281  reinterpret_cast<sFloat**>(inOdd.backGhostFaceBuffer),
282  0, daggerBit, mu);
283  }
284 }
285 
286 
287 void mat_mg4dir(cpuColorSpinorField *out, void **link, void **ghostLink, cpuColorSpinorField *in, int dagger_bit, int mu,
288  QudaPrecision sPrecision, QudaPrecision gPrecision)
289 {
290 
291  if (sPrecision == QUDA_DOUBLE_PRECISION){
292  if (gPrecision == QUDA_DOUBLE_PRECISION) {
293  Mat_mg4dir<double, double>(out, (double**)link, (double**) ghostLink, in, dagger_bit, mu);
294  } else {
295  Mat_mg4dir<double, float> (out, (float**) link, (float**) ghostLink, in, dagger_bit, mu);
296  }
297  } else {
298  if (gPrecision == QUDA_DOUBLE_PRECISION){
299  Mat_mg4dir<float, double> (out, (double**)link, (double**) ghostLink, in, dagger_bit, mu);
300  } else {
301  Mat_mg4dir<float, float> (out, (float**) link, (float**) ghostLink, in, dagger_bit, mu);
302  }
303  }
304 }
305 
306 
307 
308 void matdagmat_mg4dir(cpuColorSpinorField* out, void **link, void** ghostLink, cpuColorSpinorField* in,
309  int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision,
311 {
312  //assert sPrecision and gPrecision must be the same
313  if (sPrecision != gPrecision){
314  errorQuda("Spinor precision and gPrecison is not the same");
315  }
316 
317  QudaParity otherparity = QUDA_INVALID_PARITY;
318  if (parity == QUDA_EVEN_PARITY){
319  otherparity = QUDA_ODD_PARITY;
320  } else if (parity == QUDA_ODD_PARITY) {
321  otherparity = QUDA_EVEN_PARITY;
322  } else {
323  errorQuda("ERROR: full parity not supported in function %s\n", __FUNCTION__);
324  }
325 
326  covdev_dslash_mg4dir(tmp, link, ghostLink, in, otherparity, dagger_bit, mu, sPrecision, gPrecision);
327 
328  covdev_dslash_mg4dir(out, link, ghostLink, tmp, parity, dagger_bit, mu, sPrecision, gPrecision);
329 }
330 
331 #endif
332 
const ColorSpinorField & Odd() const
const ColorSpinorField & Even() const
void exchangeGhost(QudaParity parity, int nFace, int dagger, const MemoryLocation *pack_destination=nullptr, const MemoryLocation *halo_location=nullptr, bool gdr_send=false, bool gdr_recv=false, QudaPrecision ghost_precision=QUDA_INVALID_PRECISION) const
This is a unified ghost exchange function for doing a complete halo exchange regardless of the type o...
static void * fwdGhostFaceBuffer[QUDA_MAX_DIM]
static void * backGhostFaceBuffer[QUDA_MAX_DIM]
double mu
int Vh
Definition: host_utils.cpp:38
void covdevReference(sFloat *res, gFloat **link, sFloat *spinorField, int oddBit, int daggerBit, int mu)
void Mat(sFloat *out, gFloat **link, sFloat *in, int daggerBit, int mu)
void display_link_internal(Float *link)
void Matdagmat(sFloat *out, gFloat **link, sFloat *in, int daggerBit, int mu, sFloat *tmp, QudaParity parity)
void matdagmat(void *out, void **link, void *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision, void *tmp, QudaParity parity)
void covdev_dslash(void *res, void **link, void *spinorField, int oddBit, int daggerBit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
void mat(void *out, void **link, void *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
void * memset(void *s, int c, size_t n)
void matdagmat_mg4dir(cpuColorSpinorField *out, void **link, void **ghostLink, cpuColorSpinorField *in, int dagger_bit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision, cpuColorSpinorField *tmp, QudaParity parity)
void covdev_dslash_mg4dir(cpuColorSpinorField *out, void **link, void **ghostLink, cpuColorSpinorField *in, int oddBit, int daggerBit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
void mat_mg4dir(cpuColorSpinorField *out, void **link, void **ghostLink, cpuColorSpinorField *in, int daggerBit, int mu, QudaPrecision sPrecision, QudaPrecision gPrecision)
void ** ghostLink
Definition: covdev_test.cpp:38
QudaParity parity
Definition: covdev_test.cpp:40
cudaColorSpinorField * tmp
Definition: covdev_test.cpp:34
cpuColorSpinorField * spinor
Definition: covdev_test.cpp:31
enum QudaPrecision_s QudaPrecision
@ 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_DOUBLE_PRECISION
Definition: enum_quda.h:65
enum QudaParity_s QudaParity
#define gauge_site_size
Definition: face_gauge.cpp:34
int faceVolume[4]
Definition: host_utils.cpp:41
#define spinor_site_size
Definition: host_utils.h:9
FloatingPoint< float > Float
__host__ __device__ T sum(const array< T, s > &a)
Definition: utility.h:76
Main header file for the QUDA library.
#define errorQuda(...)
Definition: util_quda.h:120