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