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