su2hmc
Loading...
Searching...
No Matches
diagnostics.c
1#ifdef DIAGNOSTIC
2#include <assert.h>
3#include <complex.h>
4#include <float.h>
5#include <clover.h>
6#include <matrices.h>
7#include <su2hmc.h>
8#include <string.h>
9
10int Diagnostics(int istart, Complex *u[2], Complex *ut[2],Complex_f *ut_f[2],\
11 unsigned int *iu, unsigned int *id, int *hu, int *hd, double *dk[2], float *dk_f[2],\
12 const unsigned short gamin[16], const Complex gamval[20], const Complex_f gamval_f[20],\
13 const Complex *sigval, const Complex_f *sigval_f, const unsigned short *sigin,
14 Complex_f jqq,float akappa,float beta, float c_sw, double ancg){
15 /*
16 * Routine to check if the multiplication routines are working or not
17 * How I hope this will work is that
18 * 1) Initialise the system
19 * 2) Just after the initialisation of the system but before anything
20 * else call this routine using the C Preprocessor.
21 * 3) Give dummy values for the fields and then do work with them
22 * Caveats? Well this could get messy if we call something we didn't
23 * realise was being called and hadn't initialised it properly (Congradq
24 * springs to mind straight away)
25 */
26 const char *funcname = "Diagnostics";
27
28 //Initialise the arrays being used. Just going to assume MKL is being
29 //used here will also assert the number of flavours for now to avoid issues
30 //later
31 assert(nf==1);
32 printf("FLT_EVAL_METHOD is %i. Check online for what this means\n", FLT_EVAL_METHOD);
33
34 const unsigned short nclov=6;
35 unsigned int itercg=0;
36 Complex_f *clover_f[nc], *hLeaves[ndim][nc]; Complex *clover[nc];
37 Bilinear_a Xmn[nclov];
38 Complex *ut_save[nc];
39#ifdef USE_GPU
40 int device=-1;
41 cudaGetDevice(&device);
42 Complex *xi,*R1,*Phi,*X0,*X1, *smallPhi;
43 Complex_f *X0_f, *X1_f, *xi_f, *R1_f, *Phi_f;
44 double *dSdpi,*pp;
45 //Some of these strictly do not have a halo. To make things easier I'm giving them one anyway and adjusting the
46 //output compared to what might be expected in the main code ((void **)kvol vs kvolHalo mainly)
47 cudaMallocManaged((void **)clover+0,6*kvol*sizeof((void **)Complex),cudaMemAttachGlobal);
48 cudaMallocManaged((void **)clover+1,6*kvol*sizeof((void **)Complex),cudaMemAttachGlobal);
49 cudaMallocManaged((void **)clover_f+0,6*kvol*sizeof((void **)Complex_f),cudaMemAttachGlobal);
50 cudaMallocManaged((void **)clover_f+1,6*kvol*sizeof((void **)Complex_f),cudaMemAttachGlobal);
51 cudaMallocManaged((void **)&R1,kfermHalo*sizeof((void **)Complex),cudaMemAttachGlobal);
52 cudaMallocManaged((void **)&xi,kfermHalo*sizeof((void **)Complex),cudaMemAttachGlobal);
53 cudaMallocManaged((void **)&R1_f,kfermHalo*sizeof((void **)Complex_f),cudaMemAttachGlobal);
54 cudaMallocManaged((void **)&xi_f,kfermHalo*sizeof((void **)Complex_f),cudaMemAttachGlobal);
55 cudaMallocManaged((void **)&Phi,nf*kferm*sizeof((void **)Complex),cudaMemAttachGlobal);
56 cudaMallocManaged((void **)&smallPhi,kferm2*sizeof((void **)Complex),cudaMemAttachGlobal);
57 cudaMallocManaged((void **)&Phi_f,nf*kferm*sizeof((void **)Complex_f),cudaMemAttachGlobal);
58 cudaMallocManaged((void **)&X0,kferm2Halo*sizeof((void **)Complex),cudaMemAttachGlobal);
59 cudaMallocManaged((void **)&X1,kferm2Halo*sizeof((void **)Complex),cudaMemAttachGlobal);
60 cudaMallocManaged((void **)&X0_f,kferm2Halo*sizeof((void **)Complex_f),cudaMemAttachGlobal);
61 cudaMallocManaged((void **)&X1_f,kferm2Halo*sizeof((void **)Complex_f),cudaMemAttachGlobal);
62 cudaMallocManaged((void **)&X2_f,kferm2Halo*sizeof((void **)Complex_f),cudaMemAttachGlobal);
63 cudaMallocManaged((void **)&pp,kmom*sizeof((void **)double),cudaMemAttachGlobal);
64 cudaMallocManaged((void **)&dSdpi,kmom*sizeof((void **)double),cudaMemAttachGlobal);
65 cudaMallocManaged((void **)&ut_save[0],ndim*kvolHalo*sizeof(Complex),cudaMemAttachGlobal);
66 cudaMallocManaged((void **)&ut_save[1],ndim*kvolHalo*sizeof(Complex),cudaMemAttachGlobal);
67 for(unsigned short i=0;i<ndim;i++){
68 cudaMallocManaged((void **)hLeaves[i]+0,kvol*ndim*sizeof(Complex_f),cudaMemAttachGlobal);
69 cudaMallocManaged((void **)hLeaves[i]+1,kvol*ndim*sizeof(Complex_f),cudaMemAttachGlobal);
70 }
71 for(unsigned short clov=0;clov<nclov;clov++){
72 cudaMallocManaged((void**)(Xmn[clov.diag]),2*kvolHalo*sizeof(float),cudaMemAttachGlobal);
73 cudaMallocManaged((void**)(Xmn[clov.offd]),kvolHalo*sizeof(Complex_f),cudaMemAttachGlobal);
74 }
75#else
76 clover[0]=aligned_alloc(AVX,6*kvol*sizeof(Complex));
77 clover[1]=aligned_alloc(AVX,6*kvol*sizeof(Complex));
78 for(unsigned short mu=0;mu<ndim;mu++){
79 hLeaves[mu][0]=(Complex_f *)aligned_alloc(AVX,ndim*kvol*sizeof(Complex_f));
80 hLeaves[mu][1]=(Complex_f *)aligned_alloc(AVX,ndim*kvol*sizeof(Complex_f));
81 }
82 Complex *R1= aligned_alloc(AVX,kfermHalo*sizeof(Complex));
83 Complex *xi= aligned_alloc(AVX,kfermHalo*sizeof(Complex));
84 Complex_f *R1_f= aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
85 Complex_f *xi_f= aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
86 Complex *smallPhi= aligned_alloc(AVX,kferm2*sizeof(Complex));
87 Complex *Phi= aligned_alloc(AVX,nf*kferm*sizeof(Complex));
88 Complex_f *Phi_f= aligned_alloc(AVX,nf*kferm*sizeof(Complex_f));
89 Complex *X0= aligned_alloc(AVX,nf*kferm2Halo*sizeof(Complex));
90 Complex *X1= aligned_alloc(AVX,kferm2Halo*sizeof(Complex));
91 double *pp = aligned_alloc(AVX,kmom*sizeof(double));
92 Complex_f *X0_f= aligned_alloc(AVX,nf*kferm2Halo*sizeof(Complex_f));
93 Complex_f *X1_f= aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
94 Complex_f *X2_f= (Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
95 double *dSdpi = aligned_alloc(AVX,kmom*sizeof(double));
96 ut_save[0] = aligned_alloc(AVX, ndim*kvolHalo*sizeof(Complex));
97 ut_save[1] = aligned_alloc(AVX, ndim*kvolHalo*sizeof(Complex));
98 for(unsigned short clov=0;clov<nclov;clov++){
99 Xmn[clov].diag=(float *)aligned_alloc(AVX,2*kvolHalo*sizeof(float));
100 Xmn[clov].offd=(Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f));
101 }
102#endif
103 //Trial fields shouldn't get modified (except for gauge_update
104 switch(istart){
105 //Got gauge fields from file or random so print them
106 case(1):
107#pragma omp parallel sections
108 {
109#pragma omp section
110 {
111 FILE *trial_out = fopen("gauge_t", "w");
112 for(unsigned int i=0;i<(kvol+halo);i++){
113 if(i<kvol)
114 fprintf(trial_out,"Site %d:\n",i);
115 else
116 fprintf(trial_out,"Halo site %d:\n",i);
117 for(unsigned short mu=0;mu<ndim;mu++)
118 fprintf(trial_out,"Dir %d:\t%.3f+%.3fI\t%.3f+%.3fI\n", mu,\
119 creal(ut[0][i+mu*kvolHalo]),cimag(ut[0][i+mu*kvolHalo]),\
120 creal(ut[1][i+mu*kvolHalo]),cimag(ut[1][i+mu*kvolHalo]));
121 fprintf(trial_out,"\n");
122 }
123 fclose(trial_out);
124 }
125#pragma omp section
126 {
127 FILE *trial_out_f = fopen("gauge_t_f", "w");
128 for(unsigned int i=0;i<(kvol+halo);i++){
129 if(i<kvol)
130 fprintf(trial_out_f,"Site %d:\n",i);
131 else
132 fprintf(trial_out_f,"Halo site %d:\n",i);
133 for(unsigned short mu=0;mu<ndim;mu++)
134 fprintf(trial_out_f,"Dir %d:\t%.3f+%.3fI\t%.3f+%.3fI\n", mu,\
135 creal(ut_f[0][i+mu*kvolHalo]),cimag(ut_f[0][i+mu*kvolHalo]),\
136 creal(ut_f[1][i+mu*kvolHalo]),cimag(ut_f[1][i+mu*kvolHalo]));
137 fprintf(trial_out_f,"\n");
138 }
139 fclose(trial_out_f);
140 }
141 }
142 break;
143 default:
144 //Cold start as a default. Don't need to print
145 //NOTE: Single link set non unity
146 if(!rank)
147 printf("Cold Start\n");
148 u[0][0]=1+0*I; u[1][0]=0+0*I;
149 u[0][1+kvolHalo]=1+0*I; u[1][1+kvolHalo]=0+0*I;
150#pragma omp parallel for
151 for(unsigned short mu=0;mu<ndim;mu++){
152 memcpy(ut[0]+mu*kvolHalo,u[0]+mu*kvol,kvol*sizeof(Complex));
153 memcpy(ut[1]+mu*kvolHalo,u[1]+mu*kvol,kvol*sizeof(Complex));
154 }
155 break;
156 }
157 //Ensure reunitarisation is working
158 Reunitarise(ut);
159 for(unsigned short mu=0;mu<ndim;mu++)
160 for(unsigned int i=0;i<kvol;i++){
161 double diff = 1-fabs(creal(ut[0][i+mu*kvolHalo]*conj(ut[0][i+mu*kvolHalo])+ut[1][i+mu*kvolHalo]*conj(ut[1][i+mu*kvolHalo])));
162 if(diff >1e-6){
163 fprintf(stderr,"Error %i in %s: Gauge links not correctly reuniterised for site %i and direction %d. Diff %e"\
164 "\nExiting...\n\n",REUNIERR,funcname,i,mu,diff);
165 exit(REUNIERR);
166 }
167 }
168 //Check precision change works
169 ComplexConvert(ut_f[0],ut[0],kvol,true,ndim);
170 for(unsigned short mu=0;mu<ndim;mu++)
171 for(unsigned int i=0;i<kvol;i++){
172 Complex diff =ut_f[0][i+kvolHalo*mu]-ut[0][i+kvolHalo*mu];
173 if(fabs(creal(diff))>1e-6||fabs(cimag(diff))>1e-6){
174 fprintf(stderr,"Error %i in %s: Gauge links not correctly converted to float for site %i and direction %d. Diff %e+I%e"\
175 "\nExiting...\n\n",CONVERR,funcname,i,mu,creal(diff),cimag(diff));
176 exit(CONVERR);
177 }
178 }
179 //Repeat in the opposite direction.
180 ComplexConvert(ut_f[0],ut[0],kvol,false,ndim);
181 for(unsigned short mu=0;mu<ndim;mu++)
182 for(unsigned int i=0;i<kvol;i++){
183 Complex diff =ut_f[0][i+kvolHalo*mu]-ut[0][i+kvolHalo*mu];
184 if(fabs(creal(diff))>1e-6||fabs(cimag(diff))>1e-6){
185 fprintf(stderr,"Error %i in %s: Gauge links not correctly converted to double for site %i and direction %d. Diff %e+I%e"\
186 "\nExiting...\n\n",CONVERR,funcname,i,mu,creal(diff),cimag(diff));
187 exit(CONVERR);
188 }
189 }
190 //Gauge halo exchange.
191 Trial_Exchange(ut,ut_f);
192 //TODO: Figure out a test for this. May require a second lattice to be copied over in full...
193
194#pragma omp parallel sections
195 {
196#pragma omp section
197 {
198 FILE *dk4m_File = fopen("dk0","w");
199 for(int i=0;i<kvol;i+=4)
200 fprintf(dk4m_File,"%f\t%f\t%f\t%f\n",dk[0][i],dk[0][i+1],dk[0][i+2],dk[0][i+3]);
201 }
202#pragma omp section
203 {
204 FILE *dk4p_File = fopen("dk1","w");
205 for(int i=0;i<kvol;i+=4)
206 fprintf(dk4p_File,"%f\t%f\t%f\t%f\n",dk[1][i],dk[1][i+1],dk[1][i+2],dk[1][i+3]);
207 }
208 }
209
210 const int na=0;
211#ifdef LOAD_HIREP_PF
212 /* Load HiRep pseudofermion into Phi upper Gor'kov block (flavor 0, j=0..nc*ndirac-1).
213 * Format: flat binary, no header, T(slow)->X->Y->Z(fast) order,
214 * nc*ndirac=8 Complex per site, Dirac-outer/colour-inner.
215 * Lower Gor'kov block (j=nc*ndirac..nc*ngorkov-1) stays zero.
216 * Gamma note: spatial-spatial clover force matches HiRep; temporal-spatial differs by -1. */
217 {
218 FILE *pf_in = fopen("hirep_pf.bin", "rb");
219 if (!pf_in)
220 printf("hirep_pf.bin not found; running with default Phi\n");
221 else {
222 memset(Phi, 0, nf*kferm*sizeof(Complex));
223 Complex spinor[nc*ndirac];
224 for (unsigned int i = 0; i < kvol; i++) {
225 if (fread(spinor, sizeof(Complex), nc*ndirac, pf_in) != (size_t)(nc*ndirac)) {
226 fprintf(stderr, "%s: short read from hirep_pf.bin\n", funcname);
227 break;
228 }
229 for (int j = 0; j < nc*ndirac; j++){
230 Phi[i + kvol*j] = spinor[j];
231 R1[i + kvolHalo*j] = spinor[j];
232 }
233 }
234 fclose(pf_in);
235 double norm2 = 0;
236 for (int k = 0; k < nc*ndirac*kvol; k++)
237 norm2 += creal(Phi[k])*creal(Phi[k]) + cimag(Phi[k])*cimag(Phi[k]);
238 printf("HiRep PF loaded into Phi (sq.norm=%.6e)\n", norm2);
239 }
240 }
241
242#else
243 /*
244 Gauss_d(pp,kmom,0,1);
245 Gauss_c(R1_f, kferm, 0, 1/sqrt(2)); Gauss_c(Phi_f, kferm, 0, 1/sqrt(2));
246 Gauss_c(xi_f, kferm, 0, 1/sqrt(2));
247 */
248#pragma omp parallel for simd aligned(Phi,xi,R1:AVX)
249 for(unsigned int i=0;i<kvol;i++)
250 for(unsigned short j=0;j<nc*ngorkov;j++){
251 Phi[i+j*kvol]=0.0f+0.0*I; xi[i+j*kvolHalo]=0.0f+0.0*I; R1[i+j*kvolHalo]=0.0f+0.0*I;
252 }
253
254 Phi[0+(0*nc+0)*kvol]=1.0f+0.0*I; xi[0+(0*nc+0)*kvolHalo]=1.0f+0.0*I; R1[0+(0*nc+0)*kvolHalo]=1.0f+0.0*I;
255 //Phi[0+(1*nc+0)*kvol]=1.0f+0.0*I; xi[0+(1*nc+0)*kvolHalo]=1.0f+0.0*I; R1[0+(1*nc+0)*kvolHalo]=1.0f+0.0*I;
256#endif
257 ComplexConvert(Phi_f,Phi,kferm,true,1);
258 ComplexConvert(xi_f,xi,kvol,true,ngorkov);
259 ComplexConvert(R1_f,R1,kvol,true,ngorkov);
260
261 //Gauss_c(X0_f, kferm2, 0, 1/sqrt(2)); Gauss_c(X1_f, kferm2, 0, 1/sqrt(2));
262#pragma omp parallel for simd aligned(X0,X1:AVX)
263 for(unsigned int i=0;i<kvol;i++)
264 for(unsigned short j=0;j<ndirac;j++)
265 {
266 X0_f[i+j*kvolHalo]=1; xi_f[i+j*kvolHalo]=1;
267 }
268
269 ComplexConvert(X0_f,X0,kvol,false,ndirac);
270 ComplexConvert(X1_f,X1,kvol,false,ndirac);
271#pragma omp parallel for simd aligned(pp:AVX)
272 for(unsigned int i=0;i<kmom;i++)
273 pp[i]=0;
274
275 //Random nomalised momentum field
276 Gauss_d(dSdpi,kmom,0,1/sqrt(2));
277#pragma omp for simd aligned(dSdpi:AVX) nowait
278 for(int i=0; i<kmom; i+=4){
279 double norm = sqrt(dSdpi[i]*dSdpi[i]+dSdpi[i+1]*dSdpi[i+1]+dSdpi[i+2]*dSdpi[i+2]+dSdpi[i+3]*dSdpi[i+3]);
280 dSdpi[i]/=norm; dSdpi[i+1]/=norm; dSdpi[i+2]/=norm;dSdpi[i+3]/=norm;
281 }
282 FILE *input, *output;
283 FILE *input_f, *output_f;
284 FILE *input_diff, *output_diff;
285 for(int test = 0; test<=17; test++){
286 switch(test){
287 case(0): //UpDownPart
288 input = fopen("PreUpDownPart","w");
289 for(int i=0; i<kvol; i++){
290 fprintf(input,"Site %d:\t",i);
291 for(unsigned short j=0;j<nc*ndirac;j++){
292 fprintf(input,"%.5e+%.5ei\t", creal(R1[i+j*kvol]),cimag(R1[i+j*kvol]));
293 }
294 fprintf(input,"\n");
295 }
296 fclose(input);
297 UpDownPart(na,X0,R1);
298 output = fopen("UpDownPart","w");
299 for(unsigned int i=0; i<kvol; i++){
300 fprintf(output,"Site %d:\t",i);
301 for(unsigned short j=0;j<nc*ndirac;j++){
302 fprintf(output,"%.5e+%.5ei\t", creal(X0[i+j*kvol]),cimag(X0[i+j*kvol]));
303 }
304 fprintf(output,"\n");
305 }
306 fclose(output);
307 for(unsigned short idirac=0;idirac<ndirac;idirac++)
308 for(unsigned short ic=0;ic<nc;ic++)
309 for(unsigned int i=0;i<kvol;i++){
310 if(X0[i+kvol*(ic+nc*(idirac+ndirac*na))]!=R1[i+kvol*(ic+nc*idirac)]){
311 fprintf(stderr,"Error %i in %s: Up/down partitioning failed for site %d colour %d and dirac spinor %d."
312 "\nExiting...\n\n",UDPERR,funcname,i,ic,idirac);
313 exit(UDPERR);
314 }
315 }
316 break;
317 case(1): //Dslash
318 ComplexConvert(R1_f,R1,kvol,false,nc*ngorkov);
319 memset(xi,0,kfermHalo*sizeof(Complex)); memset(xi_f,0,kfermHalo*sizeof(Complex_f));
320 //NOTE: Each line corresponds to one lattice direction, in the form of colour 0, colour 1.
321 //Each block to one lattice site
322 input = fopen("dslash_in", "w"); input_f = fopen("dslash_f_in", "w"); input_diff = fopen("dslash_diff_in", "w");
323#ifdef USE_GPU
324 cudaDeviceSynchronise();
325#endif
326 for(unsigned int i = 0; i< kvol; i++){
327 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
328 for(unsigned short j=0;j<nc*ngorkov;j++){
329 fprintf(input, "%.3f+%.3fI\t",creal(R1[i+j*kvolHalo]),cimag(R1[i+j*kvolHalo]));
330 fprintf(input_f, "%.3f+%.3fI\t", creal(R1_f[i+j*kvolHalo]),cimag(R1_f[i+j*kvolHalo]));
331 fprintf(input_diff,"%.3f+%.3fI\t", creal(R1[i+j*kvolHalo]-R1_f[i+j*kvolHalo]),cimag(R1[i+j*kvolHalo]-R1_f[i+j*kvolHalo]));
332 }
333 fprintf(input, "\n\n"); fprintf(input_f,"\n\n"); fprintf(input_diff,"\n\n");
334 }
335 fclose(input); fclose(input_f); fclose(input_diff);
336 Dslash(xi,R1,ut,iu,id,gamval,gamin,dk,jqq,akappa);
337 Dslash_f(xi_f,R1_f,ut_f,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
338#ifdef USE_GPU
339 cudaDeviceSynchronise();
340#endif
341 output = fopen("dslash", "w"); output_f = fopen("dslash_f", "w"); output_diff = fopen("dslash_diff", "w");
342 for(unsigned int i = 0; i< kvol; i++){
343 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
344 for(unsigned short j=0;j<nc*ngorkov;j++){
345 fprintf(output, "%.3f+%.3fI\t",creal(xi[i+j*kvolHalo]),cimag(xi[i+j*kvolHalo]));
346 fprintf(output_f, "%.3f+%.3fI\t", creal(xi_f[i+j*kvolHalo]),cimag(xi_f[i+j*kvolHalo]));
347 Complex diff = xi[i+j*kvolHalo]-xi_f[i+j*kvolHalo];
348 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
349 fprintf(stderr,"Error %i in %s: Single and double disagree for Dslash site %i and spinor/color %d. Difference %e+%ei"\
350 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
351 fclose(output);fclose(output_f);fclose(output_diff);
352 exit(CONVERR);
353 }
354 else
355 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
356 }
357 fprintf(output, "\n\n"); fprintf(output_f,"\n\n"); fprintf(output_diff,"\n\n");
358 }
359 fclose(output); fclose(output_f); fclose(output_diff);
360 break;
361 case(2): //Dslashd
362 ComplexConvert(R1_f,R1,kvol,false,nc*ngorkov);
363 memset(xi,0,kfermHalo*sizeof(Complex)); memset(xi_f,0,kfermHalo*sizeof(Complex_f));
364 //NOTE: Each line corresponds to one lattice direction, in the form of colour 0, colour 1.
365 //Each block to one lattice site
366 input = fopen("dslashd_in", "w"); input_f = fopen("dslashd_f_in", "w"); input_diff = fopen("dslashd_diff_in", "w");
367#ifdef USE_GPU
368 cudaDeviceSynchronise();
369#endif
370 for(unsigned int i = 0; i< kvol; i++){
371 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
372 for(unsigned short j=0;j<nc*ngorkov;j++){
373 fprintf(input, "%.3f+%.3fI\t",creal(R1[i+j*kvolHalo]),cimag(R1[i+j*kvolHalo]));
374 fprintf(input_f, "%.3f+%.3fI\t", creal(R1_f[i+j*kvolHalo]),cimag(R1_f[i+j*kvolHalo]));
375 fprintf(input_diff,"%.3f+%.3fI\t", creal(R1[i+j*kvolHalo]-R1_f[i+j*kvolHalo]),cimag(R1[i+j*kvolHalo]-R1_f[i+j*kvolHalo]));
376 }
377 fprintf(input, "\n\n"); fprintf(input_f,"\n\n"); fprintf(input_diff,"\n\n");
378 }
379 fclose(input); fclose(input_f);fclose(input_diff);
380 Dslashd(xi,R1,ut,iu,id,gamval,gamin,dk,jqq,akappa);
381 Dslashd_f(xi_f,R1_f,ut_f,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
382#ifdef USE_GPU
383 cudaDeviceSynchronise();
384#endif
385 output = fopen("dslashd", "w"); output_f = fopen("dslashd_f", "w"); output_diff = fopen("dslashd_diff", "w");
386 for(unsigned int i = 0; i< kvol; i++){
387 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
388 //Note. The output of Dslashd should not have a halo. Whilst xi is defined with one we do not use it here
389 //so stride is kvol, not kvolHalo
390 for(unsigned short j=0;j<nc*ngorkov;j++){
391 fprintf(output, "%.3f+%.3fI\t",creal(xi[i+j*kvol]),cimag(xi[i+j*kvol]));
392 fprintf(output_f, "%.3f+%.3fI\t", creal(xi_f[i+j*kvol]),cimag(xi_f[i+j*kvol]));
393 Complex diff = xi[i+j*kvol]-xi_f[i+j*kvol];
394 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
395 fprintf(stderr,"Error %i in %s: Single and double disagree for Dslashd site %i and spinor/color %d. Difference %e+%ei"\
396 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
397 fclose(output);fclose(output_f);fclose(output_diff);
398 exit(CONVERR);
399 }
400 else
401 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
402 }
403 fprintf(output, "\n\n"); fprintf(output_f,"\n\n"); fprintf(output_diff,"\n\n");
404 }
405 input = fopen("dslashd_in", "w"); input_f = fopen("dslashd_f_in", "w"); input_diff = fopen("dslashd_diff_in", "w");
406 break;
407 case(3): //Hdslash
408 //NOTE: Each line corresponds to one lattice direction, in the form of colour 0, colour 1.
409 //Each block to one lattice site
410 ComplexConvert(X0_f,X0,kvol,false,nc*ndirac);
411 memset(X1,0,kferm2Halo*sizeof(Complex)); memset(X1_f,0,kferm2Halo*sizeof(Complex_f));
412 input = fopen("hdslash_in", "w"); input_f = fopen("hdslash_f_in", "w"); input_diff = fopen("hdslash_diff_in", "w");
413 for(unsigned int i = 0; i< kvol; i++){
414 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
415 for(unsigned short j=0;j<nc*ndirac;j++){
416 fprintf(input, "%.3f+%.3fI\t",creal(X0[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]));
417 fprintf(input_f, "%.3f+%.3fI\t", creal(X0_f[i+j*kvolHalo]),cimag(X0_f[i+j*kvolHalo]));
418 fprintf(input_diff,"%.3f+%.3fI\t", creal(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]));
419 }
420 fprintf(input, "\n\n"); fprintf(input_f,"\n\n"); fprintf(input_diff,"\n\n");
421 }
422 fclose(input);fclose(input_f);fclose(input_diff);
423 Hdslash(X1,X0,ut,iu,id,gamval,gamin,dk,akappa);
424 Hdslash_f(X1_f,X0_f,ut_f,iu,id,gamval_f,gamin,dk_f,akappa);
425#ifdef USE_GPU
426 cudaDeviceSynchronise();
427#endif
428 output = fopen("hdslash", "w"); output_f = fopen("hdslash_f", "w"); output_diff = fopen("hdslash_diff", "w");
429 for(unsigned int i = 0; i< kvol; i++){
430 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
431 //Note. The output of Dslashd should not have a halo. Whilst xi is defined with one we do not use it here
432 //so stride is kvol, not kvolHalo
433 for(unsigned short j=0;j<nc*ndirac;j++){
434 fprintf(output, "%.3f+%.3fI\t",creal(X1[i+j*kvolHalo]),cimag(X1[i+j*kvolHalo]));
435 fprintf(output_f, "%.3f+%.3fI\t", creal(X1_f[i+j*kvolHalo]),cimag(X1_f[i+j*kvolHalo]));
436 Complex diff = X1[i+j*kvolHalo]-X1_f[i+j*kvolHalo];
437 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
438 fprintf(stderr,"Error %i in %s: Single and double disagree for Hdslash site %i and spinor/color %d. Difference %e+%ei"\
439 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
440 fclose(output);fclose(output_f);fclose(output_diff);
441 exit(CONVERR);
442 }
443 else
444 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
445 }
446 fprintf(output, "\n\n"); fprintf(output_f,"\n\n"); fprintf(output_diff,"\n\n");
447 }
448 fclose(output);fclose(output_f);fclose(output_diff);
449 break;
450 case(4): //Hdslashd
451 ComplexConvert(X0_f,X0,kvol,false,nc*ndirac);
452 memset(X1,0,kferm2Halo*sizeof(Complex)); memset(X1_f,0,kferm2Halo*sizeof(Complex_f));
453 input = fopen("hdslashd_in", "w"); input_f = fopen("hdslashd_f_in", "w"); input_diff = fopen("hdslashd_diff_in", "w");
454#ifdef USE_GPU
455 cudaDeviceSynchronise();
456#endif
457 for(unsigned int i = 0; i< kvol; i++){
458 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
459 for(unsigned short j=0;j<nc*ndirac;j++){
460 fprintf(input, "%.3f+%.3fI\t",creal(X0[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]));
461 fprintf(input_f, "%.3f+%.3fI\t", creal(X0_f[i+j*kvolHalo]),cimag(X0_f[i+j*kvolHalo]));
462 fprintf(input_diff,"%.3f+%.3fI\t", creal(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]));
463 }
464 fprintf(input, "\n\n"); fprintf(input_f,"\n\n"); fprintf(input_diff,"\n\n");
465 }
466 fclose(input);fclose(input_f);fclose(input_diff);
467 Hdslashd(X1,X0,ut,iu,id,gamval,gamin,dk,akappa);
468 Hdslashd_f(X1_f,X0_f,ut_f,iu,id,gamval_f,gamin,dk_f,akappa);
469#ifdef USE_GPU
470 cudaDeviceSynchronise();
471#endif
472 output = fopen("hdslashd", "w"); output_f = fopen("hdslashd_f", "w"); output_diff = fopen("hdslashd_diff", "w");
473 for(unsigned int i = 0; i< kvol; i++){
474 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
475 //Note. The output of Dslashd should not have a halo. Whilst xi is defined with one we do not use it here
476 //so stride is kvol, not kvolHalo
477 for(unsigned short j=0;j<nc*ndirac;j++){
478 fprintf(output, "%.3f+%.3fI\t",creal(X1[i+j*kvol]),cimag(X1[i+j*kvol]));
479 fprintf(output_f, "%.3f+%.3fI\t", creal(X1_f[i+j*kvol]),cimag(X1_f[i+j*kvol]));
480 Complex diff = X1[i+j*kvol]-X1_f[i+j*kvol];
481 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
482 fprintf(stderr,"Error %i in %s: Single and double disagree for Hdslashd site %i and spinor/color %d. Difference %e+%ei"\
483 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
484 fclose(output);fclose(output_f);fclose(output_diff);
485 exit(CONVERR);
486 }
487 else
488 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
489 }
490 fprintf(output, "\n\n"); fprintf(output_f,"\n\n"); fprintf(output_diff,"\n\n");
491 }
492 fclose(output);fclose(output_f);fclose(output_diff);
493 break;
494 case(5): //Clover
495 if(c_sw==0)
496 break;
497 //Should really make Leaves a seperate case. But too much effort for now
498 output = fopen("Leaves","w");
499 for(unsigned int i=0;i<kvol;i++){
500 fprintf(output,"Site %d\n",i);
501 for(unsigned short mu=0;mu<ndim-1;mu++)
502 for(unsigned short nu=mu+1;nu<ndim;nu++)
503 if(mu!=nu){
504 unsigned short clov = (mu==0) ? nu-1 :mu+nu;
505 fprintf(output,"Clover %d\n",clov);
506 Complex_f Leaves[nc];
507 for(unsigned short leaf =0;leaf<ndim;leaf++){
508 Leaf(Leaves,ut_f,iu,id,i,mu,nu,leaf);
509 fprintf(output,"Leaf %d: Leaf0 = %e+I%e Leaf1=%e+I%e\n",leaf,\
510 crealf(Leaves[0]),cimagf(Leaves[0]),crealf(Leaves[1]),cimagf(Leaves[1]));
511 }
512 }
513 fprintf(output,"\n");
514 }
515 fclose(output);
516 Clover(clover_f,ut_f,iu,id);
517 output=fopen("Clover","w");
518 for(unsigned int i=0;i<kvol;i++){
519 fprintf(output,"Site %d\n",i);
520 for(unsigned short mu=0;mu<ndim-1;mu++)
521 for(unsigned short nu=mu+1;nu<ndim;nu++)
522 if(mu!=nu){
523 unsigned short clov = (mu==0) ? nu-1 :mu+nu;
524 fprintf(output,"mu %d nu %d Clover1 %e+i%e Clover2 %e+i%e\n",mu,nu,\
525 crealf(clover_f[0][i+kvol*clov]), cimagf(clover_f[0][i+kvol*clov]), crealf(clover_f[1][i+kvol*clov]),\
526 cimagf(clover_f[1][i+kvol*clov]));
527 }
528 fprintf(output,"\n");
529 }
530 fclose(output);
531 //Clover correct, Convert works so get it in double here for everywhere else
532 ComplexConvert(clover_f[0],clover[0],6*kvol,false,1);
533 ComplexConvert(clover_f[1],clover[1],6*kvol,false,1);
534 break;
535 case(6): //ByClover
536 if(c_sw==0)
537 break;
538 ComplexConvert(R1_f,R1,kvol,false,nc*ngorkov);
539 memset(xi,0,kfermHalo*sizeof(Complex)); memset(xi_f,0,kfermHalo*sizeof(Complex_f));
540 //NOTE: Each line corresponds to one lattice direction, in the form of colour 0, colour 1.
541 //Each block to one lattice site
542 input = fopen("byclover_in", "w"); input_f = fopen("byclover_f_in", "w"); input_diff = fopen("byclover_diff_in", "w");
543#ifdef USE_GPU
544 cudaDeviceSynchronise();
545#endif
546 for(unsigned int i = 0; i< kvol; i++){
547 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
548 for(unsigned short j=0;j<ngorkov;j++){
549 for(unsigned short c=0;c<nc;c++){
550 fprintf(input, "%.3f+%.3fI\t",creal(R1[i+(j*nc+c)*kvolHalo]),cimag(R1[i+(j*nc+c)*kvolHalo]));
551 fprintf(input_f, "%.3f+%.3fI\t", creal(R1_f[i+(j*nc+c)*kvolHalo]),cimag(R1_f[i+(j*nc+c)*kvolHalo]));
552 fprintf(input_diff,"%.3f+%.3fI\t", creal(R1[i+(j*nc+c)*kvolHalo]-R1_f[i+(j*nc+c)*kvolHalo]),cimag(R1[i+(j*nc+c)*kvolHalo]-R1_f[i+(j*nc+c)*kvolHalo]));
553 }
554 fprintf(input, "\n"); fprintf(input_f,"\n"); fprintf(input_diff,"\n");
555 }
556 fprintf(input, "\n"); fprintf(input_f,"\n"); fprintf(input_diff,"\n");
557 }
558 fclose(input); fclose(input_f); fclose(input_diff);
559 ByClover(xi,R1,clover,sigval,akappa,sigin,false);
560 ByClover_f(xi_f,R1_f,clover_f,sigval_f,akappa,sigin,false);
561#ifdef USE_GPU
562 cudaDeviceSynchronise();
563#endif
564 output = fopen("byclover", "w"); output_f = fopen("byclover_f", "w"); output_diff = fopen("byclover_diff", "w");
565 for(unsigned int i = 0; i< kvol; i++){
566 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
567 for(unsigned short j=0;j<ngorkov;j++){
568 for(unsigned short c=0;c<nc;c++){
569 fprintf(output, "%.3f+%.3fI\t",creal(xi[i+(j*nc+c)*kvolHalo]),cimag(xi[i+(j*nc+c)*kvolHalo]));
570 fprintf(output_f, "%.3f+%.3fI\t", creal(xi_f[i+(j*nc+c)*kvolHalo]),cimag(xi_f[i+(j*nc+c)*kvolHalo]));
571 Complex diff = xi[i+(j*nc+c)*kvolHalo]-xi_f[i+(j*nc+c)*kvolHalo];
572 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
573 fprintf(stderr,"Error %i in %s: Single and double disagree for ByClover site %i and spinor/color %d. Difference %e+%ei"\
574 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
575 fclose(output);fclose(output_f);fclose(output_diff);
576 exit(CONVERR);
577 }
578 else
579 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
580 }
581 fprintf(output, "\n"); fprintf(output_f,"\n"); fprintf(output_diff,"\n");
582 }
583 fprintf(output, "\n"); fprintf(output_f,"\n"); fprintf(output_diff,"\n");
584 }
585 fclose(output); fclose(output_f); fclose(output_diff);
586 break;
587 case(7): //HbyClover
588 if(c_sw==0)
589 break;
590 ComplexConvert(X0_f,X0,kvol,false,nc*ndirac);
591 memset(X1,0,kferm2Halo*sizeof(Complex)); memset(X1_f,0,kferm2Halo*sizeof(Complex_f));
592 input = fopen("hbyclover_in", "w"); input_f = fopen("hbyclover_f_in", "w"); input_diff = fopen("hbyclover_diff_in", "w");
593#ifdef USE_GPU
594 cudaDeviceSynchronise();
595#endif
596 for(unsigned int i = 0; i< kvol; i++){
597 fprintf(input, "Site %d:\n",i); fprintf(input_f, "Site %d:\n",i); fprintf(input_diff, "Site %d:\n",i);
598 for(unsigned short j=0;j<nc*ndirac;j++){
599 fprintf(input, "%.3f+%.3fI\t",creal(X0[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]));
600 fprintf(input_f, "%.3f+%.3fI\t", creal(X0_f[i+j*kvolHalo]),cimag(X0_f[i+j*kvolHalo]));
601 fprintf(input_diff,"%.3f+%.3fI\t", creal(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]),cimag(X0[i+j*kvolHalo]-X0_f[i+j*kvolHalo]));
602 }
603 fprintf(input, "\n\n"); fprintf(input_f,"\n\n"); fprintf(input_diff,"\n\n");
604 }
605 fclose(input);fclose(input_f);fclose(input_diff);
606 HbyClover(X1,X0,clover,sigval,akappa,sigin,false);
607 HbyClover_f(X1_f,X0_f,clover_f,sigval_f,akappa,sigin,false);
608#ifdef USE_GPU
609 cudaDeviceSynchronise();
610#endif
611 output = fopen("hbyclover", "w"); output_f = fopen("hbyclover_f", "w"); output_diff = fopen("hbyclover_diff", "w");
612 for(unsigned int i = 0; i< kvol; i++){
613 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i); fprintf(output_diff, "Site %d:\n",i);
614 //Note. The output of Dslashd should not have a halo. Whilst xi is defined with one we do not use it here
615 //so stride is kvol, not kvolHalo
616 for(unsigned short j=0;j<nc*ndirac;j++){
617 fprintf(output, "%.3f+%.3fI\t",creal(X1[i+j*kvol]),cimag(X1[i+j*kvol]));
618 fprintf(output_f, "%.3f+%.3fI\t", creal(X1_f[i+j*kvol]),cimag(X1_f[i+j*kvol]));
619 Complex diff = X1[i+j*kvol]-X1_f[i+j*kvol];
620 if(fabs(creal(diff))>5e-6 || fabs(cimag(diff))>5e-6){
621 fprintf(stderr,"Error %i in %s: Single and double disagree for HbyClover site %i and spinor/color %d. Difference %e+%ei"\
622 "\nExiting...\n\n",CONVERR,funcname,i,j,creal(diff),cimag(diff));
623 fclose(output);fclose(output_f);fclose(output_diff);
624 exit(CONVERR);
625 }
626 else
627 fprintf(output_diff,"%.3f+%.3fI\t", creal(diff),cimag(diff));
628 }
629 fprintf(output, "\n\n"); fprintf(output_f,"\n\n"); fprintf(output_diff,"\n\n");
630 }
631 fclose(output);fclose(output_f);fclose(output_diff);
632 break;
633 case(8): //Filling smallPhi
634 memset(smallPhi,0,kferm2*sizeof(Complex));
635 Fill_Small_Phi(na,smallPhi,Phi);
636 for(unsigned int i = 0; i<kvol;i++)
637 for(unsigned short idirac = 0; idirac<ndirac; idirac++)
638 for(unsigned short ic= 0; ic<nc; ic++)
639 // PHI_index=i*16+j*2+k;
640 if(cabs(smallPhi[i+kvol*(ic+nc*idirac)]-Phi[i+kvol*(ic+nc*(idirac+ngorkov*na))])>1e-6){
641 fprintf(stderr,"Error %i in %s: Failed to fill small phi correctly.\nExiting\n\n.",SPHIERR,funcname);
642 exit(SPHIERR);
643 }
644 break;
645 case(9): //Congradq
646 memset(X1,0,kferm2Halo*sizeof(Complex));
647 itercg=0;
648 if(Congradq(0,rescga,X1,smallPhi,ut,ut_f,clover_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,c_sw,&itercg)){
649 fprintf(stderr,"Error %i in %s: Congradq failed to converge.\nExiting\n\n",ITERLIM,funcname);
650 exit(ITERLIM);
651 }
652 //Not part of the Congrad test. But we need to know X1_f and X2_f later.
653 ComplexConvert(X1_f,X1,kvol,true,nc*ndirac);
654 Hdslash_f(X2_f,X1_f,ut_f,iu,id,gamval_f,gamin,dk_f,akappa);
655 if(c_sw)
656 HbyClover_f(X2_f,X1_f,clover_f,sigval_f,akappa,sigin,false);
657 output=fopen("X1_f","w"); output_f=fopen("X2_f","w");
658 for(unsigned int i = 0; i< kvol; i++){
659 fprintf(output, "Site %d:\n",i); fprintf(output_f, "Site %d:\n",i);
660 for(unsigned short c=0;c<nc;c++){
661 fprintf(output,"c %d",c);
662 fprintf(output_f,"c %d",c);
663 for(unsigned short j=0;j<ndirac;j++){
664 fprintf(output, "\t%.3f+%.3fI",creal(X1_f[i+kvolHalo*(c+nc*j)]),cimag(X1_f[i+kvolHalo*(c+nc*j)]));
665 fprintf(output_f, "\t%.3f+%.3fI",creal(X2_f[i+kvolHalo*(c+nc*j)]),cimag(X2_f[i+kvolHalo*(c+nc*j)]));
666 }
667 fprintf(output,"\n"); fprintf(output_f,"\n");
668 }
669 fprintf(output, "\n\n"); fprintf(output_f, "\n\n");
670 }
671 fclose(output); fclose(output_f);
672 break;
673 case(10): //Hamilton
674 memset(X1,0,kferm2Halo*sizeof(Complex));
675 double h,s,ancgh; h=s=ancgh=0;
676 Hamilton(&h,&s,rescgg,pp,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
677 jqq,akappa,beta,c_sw,&ancgh,0);
678 output = fopen("Hamiltonian", "w");
679 fprintf(output,"h=%e\ts=%e Congrad Iterations %.4e\n\n",h,s,ancgh);
680 for(unsigned int i = 0; i< kvol; i++){
681 fprintf(output, "Site %d:\n",i);
682 for(unsigned short j=0;j<nc*ndirac;j++){
683 fprintf(output, "%.3f+%.3fI\t",creal(X1[i+j*kvolHalo]),cimag(X1[i+j*kvolHalo]));
684 }
685 fprintf(output, "\n\n");
686 }
687 fclose(output);
688 break;
689 case(11): //Gauge Force
690 memset(dSdpi,0,kmom*sizeof(double));
691#ifdef USE_GPU
692 //cudaMemPrefetchAsync(dSdpi,kmom*sizeof(double),device,NULL);
693#endif
694 //Isolate Gauge force contribution
695 memset(dSdpi,0,kmom*sizeof(double));
696 Gauge_force(dSdpi,ut_f,iu,id,beta);
697#ifdef USE_GPU
698 cudaDeviceSynchronise();
699#endif
700 output = fopen("Gauge_Force","w");
701 for(unsigned int i = 0; i< kvol; i++){
702 fprintf(output,"Site %d:\n",i);
703 for(unsigned short gen=0;gen<nadj;gen++){
704 fprintf(output,"Gen %d:\n",gen);
705 for(unsigned int mu=0;mu<ndim;mu++){
706 fprintf(output, "%.3e\t", dSdpi[i+kvol*(gen*ndim+mu)]);
707 }
708 fprintf(output,"\n");
709 }
710 fprintf(output,"\n");
711 }
712 fclose(output);
713 break;
714 //Two force cases because of the flag. This also tests the conjugate gradient works okay
715 case(12): //Wilson Force
716 if(nproc>1){
717 fprintf(stderr,"Error %i in %s: MPI force diagnostic not implemented yet.\n\n"\
718 "Breaking and moving to next test",NOIMPL,funcname);
719 break;
720 }
721 //Isolate wilson force contribution
722 memset(dSdpi,0,kmom*sizeof(double));
723 for(unsigned short mu=0;mu<ndim-1;mu++)
724 Force_s(dSdpi,ut_f,X1_f,X2_f,gamval_f,iu,gamin,akappa,mu);
725 Force_t(dSdpi,ut_f,X1_f,X2_f,gamval_f,dk_f,iu,gamin,akappa);
726 output = fopen("Wilson_Force","w");
727 for(unsigned int i = 0; i< kvol; i++){
728 fprintf(output,"Site %d:\n",i);
729 for(unsigned short gen=0;gen<nadj;gen++){
730 fprintf(output,"Gen %d:\n",gen);
731 for(unsigned int short mu=0;mu<ndim;mu++){
732 fprintf(output, "%.3e\t", dSdpi[i+kvol*(gen*ndim+mu)]);
733 }
734 fprintf(output,"\n");
735 }
736 fprintf(output,"\n");
737 }
738 fclose(output);
739 break;
740 case(13): //Clover Half Leaves
741 if(c_sw==0)
742 break;
743 output=fopen("Half_leaves","w");
744 for(unsigned short mu=0;mu<ndim-1;mu++)
745 for(unsigned short nu=mu+1;nu<ndim;nu++){
746 Half_Leaves(hLeaves[mu],ut_f,iu,id,mu,nu);
747 Half_Leaves(hLeaves[nu],ut_f,iu,id,nu,mu);
748 fprintf(output,"mu %d nu %d\n",mu,nu);
749 for(unsigned int i=0;i<kvol;i++){
750 fprintf(output,"Site %d\n",i);
751 for(unsigned short leaf=0;leaf<ndim;leaf++)
752 fprintf(output,"leaf %d: mu-nu: hLeaf1 %e+i%e hLeaf2 %e+i%e\tnu-mu: hLeaf1 %e+i%e hLeaf2 %e+i%e\n",leaf,\
753 crealf(hLeaves[mu][0][i+kvol*leaf]), cimagf(hLeaves[mu][0][i+kvol*leaf]),\
754 crealf(hLeaves[mu][1][i+kvol*leaf]), cimagf(hLeaves[mu][1][i+kvol*leaf]),\
755 crealf(hLeaves[nu][0][i+kvol*leaf]), cimagf(hLeaves[nu][0][i+kvol*leaf]),\
756 crealf(hLeaves[nu][1][i+kvol*leaf]), cimagf(hLeaves[nu][1][i+kvol*leaf]));
757 }
758 fprintf(output,"\n");
759 }
760 fclose(output);
761 break;
762 case(14): //Xmunu
763 if(nproc>1){
764 fprintf(stderr,"Error %i in %s: MPI clover force not implemented yet.\n\n"\
765 "Breaking and moving to next test",NOIMPL,funcname);
766 break;
767 }
768 //Don't test if no clover.
769 if(c_sw==0)
770 break;
771 for(unsigned short mu=0;mu<ndim;mu++)
772 for(unsigned short nu=0;nu<ndim;nu++)
773 if(mu!=nu){
774 unsigned short clov = (mu==0) ? nu-1 : mu+nu;
775 CalcXmunu(Xmn[clov],X1_f,X2_f,sigval_f,sigin,mu,nu);
776 }
777 output = fopen("Xmunu","w");
778 for(unsigned int i=0;i<kvol;i++) {
779 fprintf(output,"Site %d\n",i);
780 for(unsigned short mu=0;mu<ndim;mu++)
781 for(unsigned short nu=0;nu<ndim;nu++){
782 unsigned short clov = (mu==0) ? nu-1 : mu+nu;
783 if(mu!=nu){
784 fprintf(output,"mu %d nu %d:\n",mu,nu);
785 fprintf(output,"%.3e\t%.3e+i%.3e\n",Xmn[clov].diag[i],creal(Xmn[clov].offd[i]),cimag(Xmn[clov].offd[i]));
786 fprintf(output,"%.3e+i%.3e\t%.3e\n",creal(conjf(Xmn[clov].offd[i])),cimag(conjf(Xmn[clov].offd[i])),Xmn[clov].diag[i+kvolHalo]);
787 }
788 }
789 fprintf(output,"\n");
790 }
791 fclose(output);
792 break;
793 case(15): //Clover Force
794 if(nproc>1){
795 fprintf(stderr,"Error %i in %s: MPI clover force not implemented yet.\n\n"\
796 "Breaking and moving to next test",NOIMPL,funcname);
797 break;
798 }
799 //Don't test if no clover.
800 if(c_sw==0)
801 break;
802 memset(dSdpi,0,kmom*sizeof(double));
803 Clov_Force(dSdpi,ut_f,X1_f,X2_f,sigval_f,sigin,iu,id,akappa);
804 output = fopen("Clover_Force","w");
805 for(unsigned int i = 0; i< kvol; i++){
806 fprintf(output,"Site %d:\n",i);
807 for(unsigned short gen=0;gen<nadj;gen++){
808 fprintf(output,"Gen %d:\n",gen);
809 for(unsigned short mu=0;mu<ndim;mu++){
810 fprintf(output, "%.3e\t", dSdpi[i+kvol*(gen*ndim+mu)]);
811 }
812 fprintf(output,"\n");
813 }
814 fprintf(output,"\n");
815 }
816 fclose(output);
817 break;
818 case(16): //Congradp
819 itercg=0;
820 if(Congradp(0, respbp, Phi, R1,ut,ut_f,clover_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,c_sw,&itercg)){
821 fprintf(stderr,"Error %i in %s: Congradp failed to converge.\nExiting\n\n",ITERLIM,funcname);
822 exit(ITERLIM);
823 }
824 break;
825 case(17): //Finite difference check. Produced by Claude Code Opus 4.7
826 //Build clover
827 if(c_sw){
828 free(clover_f[0]); free(clover_f[1]);
829 Clover(clover_f, ut_f, iu, id);
830 }
831 //Gaussian @f$\xi@f$ → R (ngorkov, with halo stride)
832 for(unsigned short j=0;j<nc*ngorkov;j++)
833 Gauss_c(xi_f+j*kvolHalo, kvol, 0, 1/sqrt(2));
834 //@f$\Phi=M^\dagger\xi@f$
835 Dslashd_f(R1_f, xi_f, ut_f, iu, id, gamval_f, gamin, dk_f, jqq, akappa);
836 if(c_sw)
837 ByClover_f(R1_f, xi_f, clover_f, sigval_f, akappa, sigin, true);
838 //Convert and store @f$\Phi$, populate X0 with upper half
839 for(int i=0;i<kferm;i++) R1[i] = (Complex)R1_f[i];
840 memcpy(Phi, R1, kferm*sizeof(Complex));
841 UpDownPart(0, X0, R1);
842
843 //Save original gauge fields
844 memcpy(ut_save[0], ut[0], ndim*kvolHalo*sizeof(Complex));
845 memcpy(ut_save[1], ut[1], ndim*kvolHalo*sizeof(Complex));
846
847 double h0, s0, h1, s1, ancgt=0;
848
849 //(1) Baseline: S(U_0). Use res=rescgg for tight CG.
850 memset(pp, 0, kmom*sizeof(double)); //hp=0 so s = S
851 memset(X1, 0, kferm2Halo*sizeof(Complex));
852 Hamilton(&h0, &s0, rescgg, pp, X0, X1, Phi, ut, ut_f, iu, id,
853 gamval, gamval_f, gamin, sigval, sigval_f, sigin,
854 dk, dk_f, jqq, akappa, beta, c_sw, &ancgt, 0);
855 //Hamilton wrote the CG solution to X0 → first Force call can use iflag=1
856
857 //(2) Force at U_0 (same X0, same Phi → same action functional)
858 memset(dSdpi, 0, kmom*sizeof(double));
859 Force(dSdpi, 1, rescgg, X0, X1, Phi, ut, ut_f, iu, id,
860 gamval, gamval_f, gamin, sigval, sigval_f, sigin,
861 dk, dk_f, jqq, akappa, beta, c_sw, &ancgt);
862
863 //(3) |dSdpi|^2
864 double fnorm2 = 0;
865 for(int i=0; i<kmom; i++) fnorm2 += dSdpi[i]*dSdpi[i];
866 if(nproc>1) Par_dsum(&fnorm2);
867
868 //(4) Sweep \varepsilon: take the force as the momentum direction
869 char output_name[64];
870 sprintf(output_name,"Force_Action_Check_%1.2f",c_sw);
871 output = fopen(output_name,"w");
872 fprintf(output,"|dSdpi|^2 = %.10e\n", fnorm2);
873 fprintf(output,"eps\tdS_num\tdS_ana\tratio\t(num-ana)\n");
874
875 for(int k=0; k<126; k++){
876 double eps = 1e-2 - k*(1.0/12800.0); // 1e-2, 5e-3, 2.5e-3, ...
877 memcpy(pp, dSdpi, kmom*sizeof(double)); // pp = force direction
878
879 //Restore U and move U by @f$\varepsilon@f$ along pp
880 memcpy(ut[0], ut_save[0], ndim*kvolHalo*sizeof(Complex));
881 memcpy(ut[1], ut_save[1], ndim*kvolHalo*sizeof(Complex));
882 Gauge_Update(eps, pp, ut, ut_f); // U ← exp(i \varepsilon pp T)U
883
884 //@f$S(U_\varepsilon)@f$. Fresh X0 copy so the CG initial guess is deterministic.
885 //Use Phi unchanged — same pseudofermion functional.
886 memset(X1, 0, kferm2Halo*sizeof(Complex));
887 Hamilton(&h1, &s1, rescgg, pp, X0, X1, Phi, ut, ut_f, iu, id,
888 gamval, gamval_f, gamin, sigval, sigval_f, sigin,
889 dk, dk_f, jqq, akappa, beta, c_sw, &ancgt, 0);
890 // d pp was overwritten by Hamilton? No — Hamilton doesn't modify pp.
891 // But hp = |pp|^2/2 is nonzero now. Use s1, not h1.
892
893 double dS_num = s1 - s0;
894 double dS_ana = eps * fnorm2;
895 fprintf(output,"%.5e\t%.10e\t%.10e\t%.6f\t%.3e\n",
896 eps, dS_num, dS_ana, dS_num/dS_ana,
897 (dS_num - dS_ana)/(eps*eps));
898 }
899 fclose(output);
900
901 break;
902
903 }
904 }
905 //George Michael's favourite bit of the code
906#ifdef USE_GPU
907 //Make a routine that does this for us
908 cudaFree(dk[0]); cudaFree(dk[1]); cudaFree(R1); cudaFree(dSdpi); cudaFree(pp);
909 cudaFree(Phi); cudaFree(ut[0]); cudaFree(ut[1]);
910 cudaFree(Phi_f); cudaFree(xi_f); cudaFree(R1_f);
911 cudaFree(clover[0]); cudaFree(clover[1]);
912 cudaFree(clover_f[0]); cudaFree(clover_f[1]);
913 cudaFree(X0); cudaFree(X1); cudaFree(u[0]); cudaFree(u[1]);
914 cudaFree(X0_f); cudaFree(X1_f); cudaFree(ut_f[0]); cudaFree(ut_f[1]);
915 cudaFree(X2_f);
916 for(unsigned short i=0;i<ndim;i++){
917 cudaFree(hLeaves[i][0]); cudaFree(hLeaves[i][1]);
918 }
919 for(unsigned short clov=0;clov<nclov;clov++){
920 cudaFree(Xmn[clov].diag); cudaFree(Xmn[clov].offd);
921 }
922 cudaFree(id); cudaFree(iu); cudaFree(hd); cudaFree(hu);
923 cudaFree(ut_save[0]); cudaFree(ut_save[1]);
924#else
925 free(dk[0]); free(dk[1]); free(R1); free(dSdpi); free(pp);
926 free(Phi); free(ut[0]); free(ut[1]); free(xi);
927 free(Phi_f); free(xi_f); free(R1_f);
928 free(clover[0]); free(clover[1]);
929 free(clover_f[0]); free(clover_f[1]);
930 for(unsigned short i=0;i<ndim;i++){
931 free(hLeaves[i][0]); free(hLeaves[i][1]);
932 }
933 for(unsigned short clov=0;clov<nclov;clov++){
934 free(Xmn[clov].diag); free(Xmn[clov].offd);
935 }
936 free(X0); free(X1); free(u[0]); free(u[1]);
937 free(X2_f);
938 free(id); free(iu); free(hd); free(hu);
939 free(ut_save[0]); free(ut_save[1]);
940 free(pcoord);
941#endif
942
943#if(nproc>1)
944 MPI_Finalise();
945#endif
946 exit(0);
947}
948#endif
clover.h
Routines needed for Clover improved wilson fermions.
hd
unsigned int * hd
Down halo indices.
Definition coord.c:11
hu
unsigned int * hu
Up halo indices.
Definition coord.c:11
SPHIERR
#define SPHIERR
Up/down partitioning failed.
Definition errorcodes.h:132
REUNIERR
#define REUNIERR
Gauge link reunitarisation failed.
Definition errorcodes.h:126
ITERLIM
#define ITERLIM
Exceeded max number of iterations.
Definition errorcodes.h:137
NOIMPL
#define NOIMPL
Not implemented.
Definition errorcodes.h:165
CONVERR
#define CONVERR
Failed to convert precision correctly.
Definition errorcodes.h:128
UDPERR
#define UDPERR
Up/down partitioning failed.
Definition errorcodes.h:130
CalcXmunu
void CalcXmunu(Bilinear_a Xmunu, Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, const unsigned short mu, const unsigned short nu)
Gets for the clover force.
Definition clover.c:437
Clov_Force
void Clov_Force(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, const Complex_f *sigval, const unsigned short *sigin, unsigned int *iu, unsigned int *id, const float akappa)
Gets the clover contribution to the force.
Definition clover.c:534
HbyClover
void HbyClover(Complex *phi, Complex *r, Complex *clover[2], Complex *sigval, const float akappa, unsigned short *sigin, bool dag)
Clover analogue of the Dslash operation. This version acts on all flavours similar to Dslash and Dsla...
Definition clover.c:287
ByClover
void ByClover(Complex *phi, Complex *r, Complex *clover[2], Complex *sigval, const float akappa, unsigned short *sigin, bool dag)
Clover analogue of the Dslash operation. This version acts on all flavours similar to Dslash and Dsla...
Definition clover.c:243
ByClover_f
void ByClover_f(Complex_f *phi, Complex_f *r, Complex_f *clover[2], Complex_f *sigval, const float akappa, unsigned short *sigin, bool dag)
Clover analogue of the Dslash operation. This version acts on all flavours similar to Dslash and Dsla...
Definition clover.c:332
HbyClover_f
void HbyClover_f(Complex_f *phi, Complex_f *r, Complex_f *clover[2], Complex_f *sigval, const float akappa, unsigned short *sigin, bool dag)
Clover analogue of the Dslash operation. This version acts on all flavours similar to Dslash and Dsla...
Definition clover.c:376
Leaf
void Leaf(Complex_f Leaves[nc], Complex_f *ut[nc], unsigned int *iu, unsigned int *id, unsigned int i, const unsigned short mu, const unsigned short nu, const unsigned short leaf)
Calculates a leaf for a clover term.
Definition clover.c:135
Half_Leaves
void Half_Leaves(Complex_f *hLeaves[2], Complex_f *ut[2], unsigned int *iu, unsigned int *id, const unsigned short mu, const unsigned short nu)
Calculates the products of the first two links in a plaquette.
Definition clover.c:122
Clover
void Clover(Complex_f *clover[2], Complex_f *ut[2], unsigned int *iu, unsigned int *id)
Calculates the clovers in all directions at all sites.
Definition clover.c:203
Dslash_f
int Dslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa)
Evaluates in single precision.
Definition matrices.c:425
Dslashd_f
int Dslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], Complex_f jqq, float akappa)
Evaluates in single precision.
Definition matrices.c:544
Hdslashd_f
int Hdslashd_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa)
Evaluates in single precision.
Definition matrices.c:746
Hdslashd
int Hdslashd(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa)
Evaluates in double precision.
Definition matrices.c:338
Hdslash
int Hdslash(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], float akappa)
Evaluates in double precision.
Definition matrices.c:254
Dslashd
int Dslashd(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa)
Evaluates in double precision.
Definition matrices.c:135
Hdslash_f
int Hdslash_f(Complex_f *phi, Complex_f *r, Complex_f *ut[nc], unsigned int *iu, unsigned int *id, Complex_f gamval[20], const unsigned short gamin[16], float *dk[nc], float akappa)
Evaluates in single precision.
Definition matrices.c:663
Dslash
int Dslash(Complex *phi, Complex *r, Complex *ut[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], const unsigned short gamin[16], double *dk[nc], Complex_f jqq, float akappa)
Evaluates in double precision.
Definition matrices.c:16
ComplexConvert
int ComplexConvert(Complex_f *a, Complex *b, const unsigned int len, const bool dtof, const unsigned short stride)
takes an array of complex float and double precision numbers and converts the precision
Definition coord.c:420
conj
__device__ __forceinline__ T conj(const T &z)
Complex Conjugation.
Definition cusu2hmc.cu:33
UpDownPart
int UpDownPart(const unsigned int na, Complex *X0, Complex *R1)
Up/Down partitioning of the pseudofermion field.
Definition su2hmc.c:327
Reunitarise
int Reunitarise(Complex *ut[2])
Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.
Definition su2hmc.c:342
Fill_Small_Phi
int Fill_Small_Phi(int na, Complex *smallPhi, Complex *Phi)
Copies necessary (2*4*kvol) elements of Phi into a vector variable.
Definition su2hmc.c:311
Congradp
int Congradp(int na, double res, Complex *Phi, Complex *xi, Complex *ud[2], Complex_f *ut[2], Complex_f *clover_f[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float c_sw, int *itercg)
Matrix Inversion via Conjugate Gradient (no up/down flavour partitioning). Solves The matrix multipl...
Definition congrad.c:735
Force_t
void Force_t(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], float *dk[2], unsigned int *iu, const unsigned short gamin[16], float akappa)
Calculates the force at each intermediate time.
Definition force.c:156
Force_s
void Force_s(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20], unsigned int *iu, const unsigned short gamin[16], const float akappa, const unsigned short mu)
Calculates the force at each intermediate time.
Definition force.c:81
Congradq
int Congradq(int na, double res, Complex *X1, Complex *r, Complex *ud[2], Complex_f *ut[2], Complex_f *clover_f[nc], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float c_sw, int *itercg)
Matrix Inversion via Conjugate Gradient (up/down flavour partitioning). Solves Implements up/down pa...
Definition congrad.c:278
Gauge_force
int Gauge_force(double *dSdpi, Complex_f *ut[2], unsigned int *iu, unsigned int *id, float beta)
Calculates the gauge force due to the Wilson Action at each intermediate time.
Definition force.c:11
Gauge_Update
int Gauge_Update(const double d, double *pp, Complex *ut[2], Complex_f *ut_f[2])
Gauge update for the integration step of the HMC.
Definition integrate.c:37
Force
int Force(double *dSdpi, const bool iflag, double res1, Complex *X0, Complex *X1, Complex *Phi, Complex *ut[2], Complex_f *ut_f[2], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], const Complex_f jqq, const float akappa, const float beta, const float c_sw, double *ancg)
Calculates the force at each intermediate time.
Definition force.c:228
Trial_Exchange
int Trial_Exchange(Complex *ut[2], Complex_f *ut_f[2])
Exchanges the trial fields.
Definition par_mpi.c:1043
Par_dsum
int Par_dsum(double *dval)
Performs a reduction on a double dval to get a sum which is then distributed to all ranks.
Hamilton
int Hamilton(double *h, double *s, double res2, double *pp, Complex *X0, Complex *X1, Complex *Phi, Complex *ud[2], Complex_f *ut[2], unsigned int *iu, unsigned int *id, Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16], Complex *sigval, Complex_f *sigval_f, unsigned short *sigin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float beta, float c_sw, double *ancgh, int traj)
Calculate the Hamiltonian.
Definition su2hmc.c:168
Gauss_c
int Gauss_c(Complex_f *ps, unsigned int n, const Complex_f mu, const float sigma)
Generates a vector of normally distributed random single precision complex numbers using the Box-Mull...
Definition random.c:136
Gauss_d
int Gauss_d(double *ps, unsigned int n, const double mu, const double sigma)
Generates a vector of normally distributed random double precision numbers using the Box-Muller Metho...
Definition random.c:170
matrices.h
Matrix multiplication and related declarations.
rank
int rank
The MPI rank.
Definition par_mpi.c:20
MPI_Finalise
#define MPI_Finalise()
Avoid any accidents with US/UK spelling.
Definition par_mpi.h:32
pcoord
int * pcoord
The processor grid.
Definition par_mpi.c:17
AVX
#define AVX
Alignment of arrays. 64 for AVX-512, 32 for AVX/AVX2. 16 for SSE. Since AVX is standard on modern x86...
Definition sizes.h:279
nc
#define nc
Colours.
Definition sizes.h:182
rescgg
#define rescgg
Conjugate gradient residue for update.
Definition sizes.h:251
nproc
#define nproc
Number of processors for MPI.
Definition sizes.h:138
kmom
#define kmom
sublattice momentum sites
Definition sizes.h:193
ngorkov
#define ngorkov
Gor'kov indices.
Definition sizes.h:190
kferm2Halo
#define kferm2Halo
Dirac lattice and halo.
Definition sizes.h:238
nadj
#define nadj
adjacent spatial indices
Definition sizes.h:184
kvol
#define kvol
Sublattice volume.
Definition sizes.h:163
Complex
#define Complex
Double precision complex number.
Definition sizes.h:64
kferm
#define kferm
sublattice size including Gor'kov indices
Definition sizes.h:195
rescga
#define rescga
Conjugate gradient residue for acceptance.
Definition sizes.h:253
nf
#define nf
Fermion flavours (double it).
Definition sizes.h:160
ndirac
#define ndirac
Dirac indices.
Definition sizes.h:186
cudaDeviceSynchronise
#define cudaDeviceSynchronise()
Get rid of that bastardised yankee English.
Definition sizes.h:53
respbp
#define respbp
Conjugate gradient residue for .
Definition sizes.h:249
halo
#define halo
Total Halo size.
Definition sizes.h:231
Complex_f
#define Complex_f
Single precision complex number.
Definition sizes.h:62
ndim
#define ndim
Dimensions.
Definition sizes.h:188
kferm2
#define kferm2
sublattice size including Dirac indices
Definition sizes.h:197
kvolHalo
#define kvolHalo
Subvolume + halo size.
Definition sizes.h:234
kfermHalo
#define kfermHalo
Gor'kov lattice and halo.
Definition sizes.h:236
Bilinear_a
Structure of arrays for Hermitian bilinear in memory.
Definition clover.h:30
Bilinear_a::offd
Complex_f * offd
Complex valued off-diagonal terms. We only need to store one of these to get the other in .
Definition clover.h:34
Bilinear_a::diag
float * diag
Real valued diagonal terms.
Definition clover.h:32
su2hmc.h
Function declarations for most of the routines.
creal
#define creal(z)
Extract Real Component using C standard notation.
Definition thrust_complex.h:26
cimag
#define cimag(z)
Extract Imaginary Component using C standard notation.
Definition thrust_complex.h:24
I
#define I
Define I in double precision using C standard notation.
Definition thrust_complex.h:28