QUDA  1.0.0
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
contract_reference.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <blas_reference.h>
4 #include <quda_internal.h>
5 #include "color_spinor_field.h"
6 #include <test_util.h>
7 
8 extern int Z[4];
9 extern int Vh;
10 extern int V;
11 
12 using namespace quda;
13 using namespace std;
14 
15 template <typename Float> void contractDegrandRossi(Float *h_result_)
16 {
17 
18  // Put data in complex form
19  complex<Float> temp[16];
20  complex<Float> *h_result = (complex<Float> *)(Float *)(h_result_);
21  complex<Float> I(0.0, 1.0);
22 
23  for (int site = 0; site < V; site++) {
24 
25  // Spin contract: <\phi(x)_{\mu} \Gamma_{mu,nu}^{rho,tau} \phi(y)_{\nu}>
26  // The rho index runs slowest.
27  // Layout is defined in enum_quda.h: G_idx = 4*rho + tau
28 
29  int idx = 16 * site;
30 
31  int G_idx = 0;
32  // G_idx = 0: I
33  temp[G_idx] = 0.0;
34  temp[G_idx] += h_result[idx + 4 * 0 + 0];
35  temp[G_idx] += h_result[idx + 4 * 1 + 1];
36  temp[G_idx] += h_result[idx + 4 * 2 + 2];
37  temp[G_idx] += h_result[idx + 4 * 3 + 3];
38  G_idx++;
39 
40  // G_idx = 1: \gamma_1
41  temp[G_idx] = 0.0;
42  temp[G_idx] += I * h_result[idx + 4 * 0 + 3];
43  temp[G_idx] += I * h_result[idx + 4 * 1 + 2];
44  temp[G_idx] -= I * h_result[idx + 4 * 2 + 1];
45  temp[G_idx] -= I * h_result[idx + 4 * 3 + 0];
46  G_idx++;
47 
48  // G_idx = 2: \gamma_2
49  temp[G_idx] = 0.0;
50  temp[G_idx] -= h_result[idx + 4 * 0 + 3];
51  temp[G_idx] += h_result[idx + 4 * 1 + 2];
52  temp[G_idx] += h_result[idx + 4 * 2 + 1];
53  temp[G_idx] -= h_result[idx + 4 * 3 + 0];
54  G_idx++;
55 
56  // G_idx = 3: \gamma_3
57  temp[G_idx] = 0.0;
58  temp[G_idx] += I * h_result[idx + 4 * 0 + 2];
59  temp[G_idx] -= I * h_result[idx + 4 * 1 + 3];
60  temp[G_idx] -= I * h_result[idx + 4 * 2 + 0];
61  temp[G_idx] += I * h_result[idx + 4 * 3 + 1];
62  G_idx++;
63 
64  // G_idx = 4: \gamma_4
65  temp[G_idx] = 0.0;
66  temp[G_idx] += h_result[idx + 4 * 0 + 2];
67  temp[G_idx] += h_result[idx + 4 * 1 + 3];
68  temp[G_idx] += h_result[idx + 4 * 2 + 0];
69  temp[G_idx] += h_result[idx + 4 * 3 + 1];
70  G_idx++;
71 
72  // G_idx = 5: \gamma_5
73  temp[G_idx] = 0.0;
74  temp[G_idx] += h_result[idx + 4 * 0 + 0];
75  temp[G_idx] += h_result[idx + 4 * 1 + 1];
76  temp[G_idx] -= h_result[idx + 4 * 2 + 2];
77  temp[G_idx] -= h_result[idx + 4 * 3 + 3];
78  G_idx++;
79 
80  // G_idx = 6: \gamma_5\gamma_1
81  temp[G_idx] = 0.0;
82  temp[G_idx] += I * h_result[idx + 4 * 0 + 3];
83  temp[G_idx] += I * h_result[idx + 4 * 1 + 2];
84  temp[G_idx] += I * h_result[idx + 4 * 2 + 1];
85  temp[G_idx] += I * h_result[idx + 4 * 3 + 0];
86  G_idx++;
87 
88  // G_idx = 7: \gamma_5\gamma_2
89  temp[G_idx] = 0.0;
90  temp[G_idx] -= h_result[idx + 4 * 0 + 3];
91  temp[G_idx] += h_result[idx + 4 * 1 + 2];
92  temp[G_idx] -= h_result[idx + 4 * 2 + 1];
93  temp[G_idx] += h_result[idx + 4 * 3 + 0];
94  G_idx++;
95 
96  // G_idx = 8: \gamma_5\gamma_3
97  temp[G_idx] = 0.0;
98  temp[G_idx] += I * h_result[idx + 4 * 0 + 2];
99  temp[G_idx] -= I * h_result[idx + 4 * 1 + 3];
100  temp[G_idx] += I * h_result[idx + 4 * 2 + 0];
101  temp[G_idx] -= I * h_result[idx + 4 * 3 + 1];
102  G_idx++;
103 
104  // G_idx = 9: \gamma_5\gamma_4
105  temp[G_idx] = 0.0;
106  temp[G_idx] += h_result[idx + 4 * 0 + 2];
107  temp[G_idx] += h_result[idx + 4 * 1 + 3];
108  temp[G_idx] -= h_result[idx + 4 * 2 + 0];
109  temp[G_idx] -= h_result[idx + 4 * 3 + 1];
110  G_idx++;
111 
112  // G_idx = 10: (i/2) * [\gamma_1, \gamma_2]
113  temp[G_idx] = 0.0;
114  temp[G_idx] += h_result[idx + 4 * 0 + 0];
115  temp[G_idx] -= h_result[idx + 4 * 1 + 1];
116  temp[G_idx] += h_result[idx + 4 * 2 + 2];
117  temp[G_idx] -= h_result[idx + 4 * 3 + 3];
118  G_idx++;
119 
120  // G_idx = 11: (i/2) * [\gamma_1, \gamma_3]
121  temp[G_idx] = 0.0;
122  temp[G_idx] -= I * h_result[idx + 4 * 0 + 2];
123  temp[G_idx] -= I * h_result[idx + 4 * 1 + 3];
124  temp[G_idx] += I * h_result[idx + 4 * 2 + 0];
125  temp[G_idx] += I * h_result[idx + 4 * 3 + 1];
126  G_idx++;
127 
128  // G_idx = 12: (i/2) * [\gamma_1, \gamma_4]
129  temp[G_idx] = 0.0;
130  temp[G_idx] -= h_result[idx + 4 * 0 + 1];
131  temp[G_idx] -= h_result[idx + 4 * 1 + 0];
132  temp[G_idx] += h_result[idx + 4 * 2 + 3];
133  temp[G_idx] += h_result[idx + 4 * 3 + 2];
134  G_idx++;
135 
136  // G_idx = 13: (i/2) * [\gamma_2, \gamma_3]
137  temp[G_idx] = 0.0;
138  temp[G_idx] += h_result[idx + 4 * 0 + 1];
139  temp[G_idx] += h_result[idx + 4 * 1 + 0];
140  temp[G_idx] += h_result[idx + 4 * 2 + 3];
141  temp[G_idx] += h_result[idx + 4 * 3 + 2];
142  G_idx++;
143 
144  // G_idx = 14: (i/2) * [\gamma_2, \gamma_3]
145  temp[G_idx] = 0.0;
146  temp[G_idx] -= I * h_result[idx + 4 * 0 + 1];
147  temp[G_idx] += I * h_result[idx + 4 * 1 + 0];
148  temp[G_idx] += I * h_result[idx + 4 * 2 + 3];
149  temp[G_idx] -= I * h_result[idx + 4 * 3 + 2];
150  G_idx++;
151 
152  // G_idx = 14: (i/2) * [\gamma_2, \gamma_3]
153  temp[G_idx] = 0.0;
154  temp[G_idx] -= h_result[idx + 4 * 0 + 0];
155  temp[G_idx] -= h_result[idx + 4 * 1 + 1];
156  temp[G_idx] += h_result[idx + 4 * 2 + 2];
157  temp[G_idx] += h_result[idx + 4 * 3 + 3];
158  G_idx++;
159 
160  // Replace data in h_result with spin contracted data
161  for (int i = 0; i < 16; i++) h_result[idx + i] = temp[i];
162  }
163 }
164 
165 template <typename Float> void contractColor(Float *spinorX, Float *spinorY, Float *h_result)
166 {
167 
168  Float re = 0.0, im = 0.0;
169 
170  // Conjugate spinorX
171  for (int i = 0; i < V; i++) {
172  for (int s1 = 0; s1 < 4; s1++) {
173  for (int s2 = 0; s2 < 4; s2++) {
174 
175  re = im = 0.0;
176  for (int c = 0; c < 3; c++) {
177  re += (((Float *)spinorX)[24 * i + 6 * s1 + 2 * c + 0] * ((Float *)spinorY)[24 * i + 6 * s2 + 2 * c + 0]
178  + ((Float *)spinorX)[24 * i + 6 * s1 + 2 * c + 1] * ((Float *)spinorY)[24 * i + 6 * s2 + 2 * c + 1]);
179 
180  im += (((Float *)spinorX)[24 * i + 6 * s1 + 2 * c + 0] * ((Float *)spinorY)[24 * i + 6 * s2 + 2 * c + 1]
181  - ((Float *)spinorX)[24 * i + 6 * s1 + 2 * c + 1] * ((Float *)spinorY)[24 * i + 6 * s2 + 2 * c + 0]);
182  }
183 
184  ((Float *)h_result)[2 * (i * 16 + 4 * s1 + s2) + 0] = re;
185  ((Float *)h_result)[2 * (i * 16 + 4 * s1 + s2) + 1] = im;
186  }
187  }
188  }
189 }
190 
191 template <typename Float>
192 int contraction_reference(Float *spinorX, Float *spinorY, Float *d_result, QudaContractType cType, int X[])
193 {
194 
195  int faults = 0;
196  Float tol = (sizeof(Float) == sizeof(double) ? 1e-9 : 2e-5);
197  void *h_result = malloc(V * 2 * 16 * sizeof(Float));
198 
199  // compute spin elementals
200  contractColor(spinorX, spinorY, (Float *)h_result);
201 
202  // Apply gamma insertion on host spin elementals
203  if (cType == QUDA_CONTRACT_TYPE_DR) contractDegrandRossi((Float *)h_result);
204 
205  // compare each contraction
206  for (int j = 0; j < 16; j++) {
207  bool pass = true;
208  for (int i = 0; i < V; i++) {
209  if (abs(((Float *)h_result)[32 * i + 2 * j] - ((Float *)d_result)[32 * i + 2 * j]) > tol) {
210  faults++;
211  pass = false;
212  // printfQuda("Contraction %d %d failed\n", i, j);
213  } else {
214  // printfQuda("Contraction %d %d passed\n", i, j);
215  }
216  // printfQuda("%.16f %.16f\n", ((Float*)h_result)[32*i + 2*j],((Float*)d_result)[32*i + 2*j]);
217  if (abs(((Float *)h_result)[32 * i + 2 * j + 1] - ((Float *)d_result)[32 * i + 2 * j + 1]) > tol) {
218  faults++;
219  pass = false;
220  // printfQuda("Contraction %d %d failed\n", i, j);
221  } else {
222  // printfQuda("Contraction %d %d passed\n", i, j);
223  }
224  // printfQuda("%.16f %.16f\n", ((Float*)h_result)[32*i+2*j+1],((Float*)d_result)[32*i+2*j+1]);
225  }
226  if (pass)
227  printfQuda("Contraction %d passed\n", j);
228  else
229  printfQuda("Contraction %d failed\n", j);
230  }
231 
232  free(h_result);
233  return faults;
234 };
Z
int Z[4]
Definition: test_util.cpp:26
color_spinor_field.h
test_util.h
std
STL namespace.
tol
double tol
Definition: test_util.cpp:1656
contraction_reference
int contraction_reference(Float *spinorX, Float *spinorY, Float *d_result, QudaContractType cType, int X[])
Definition: contract_reference.h:192
blas_reference.h
X
int X[4]
Definition: covdev_test.cpp:70
V
int V
Definition: test_util.cpp:27
contractDegrandRossi
void contractDegrandRossi(Float *h_result_)
Definition: contract_reference.h:15
printfQuda
#define printfQuda(...)
Definition: util_quda.h:115
QudaContractType
enum QudaContractType_s QudaContractType
quda::abs
__host__ __device__ ValueType abs(ValueType x)
Definition: complex_quda.h:125
contractColor
void contractColor(Float *spinorX, Float *spinorY, Float *h_result)
Definition: contract_reference.h:165
QUDA_CONTRACT_TYPE_DR
Definition: enum_quda.h:498
quda_internal.h
Vh
int Vh
Definition: test_util.cpp:28
Generated by   doxygen 1.8.13