su2hmc
Loading...
Searching...
No Matches
force.c
Go to the documentation of this file.
1
8#include <matrices.h>
9#include <clover.h>
10
11int Gauge_force(double *dSdpi, Complex_f *ut[2],unsigned int *iu,unsigned int *id, float beta){
12 const char funcname[] = "Gauge_force";
13
14 //We define zero halos for debugging
15 // #ifdef _DEBUG
16 // memset(ut[0][kvol], 0, ndim*halo*sizeof(Complex_f));
17 // memset(ut[1][kvol], 0, ndim*halo*sizeof(Complex_f));
18 // #endif
19 //Was a trial field halo exchange here at one point.
20#ifdef USE_GPU
21 cuGauge_force(ut,dSdpi,beta,iu,id,dimGrid,dimBlock);
22 cudaDeviceSynchronise();
23#else
24 Complex_f *Sigma[2], *ush[2];
25 Sigma[0] = (Complex_f *)aligned_alloc(AVX,kvol*sizeof(Complex_f));
26 Sigma[1]= (Complex_f *)aligned_alloc(AVX,kvol*sizeof(Complex_f));
27 ush[0] = (Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f));
28 ush[1] = (Complex_f *)aligned_alloc(AVX,kvolHalo*sizeof(Complex_f));
29 //Holders for directions
30 for(int mu=0; mu<ndim; mu++){
31 memset(Sigma[0],0, kvol*sizeof(Complex_f));
32 memset(Sigma[1],0, kvol*sizeof(Complex_f));
33 for(int nu=0; nu<ndim; nu++)
34 if(nu!=mu){
35 //The +ν Staple
36#pragma omp parallel for simd //aligned(ut[0],ut[1],Sigma[0],Sigma[1],iu:AVX)
37 for(int i=0;i<kvol;i++){
38 int uidm = iu[mu*kvol+i];
39 int uidn = iu[nu*kvol+i];
40 Complex_f a11=ut[0][uidm+kvolHalo*nu]*conj(ut[0][uidn+kvolHalo*mu])+\
41 ut[1][uidm+kvolHalo*nu]*conj(ut[1][uidn+kvolHalo*mu]);
42 Complex_f a12=-ut[0][uidm+kvolHalo*nu]*ut[1][uidn+kvolHalo*mu]+\
43 ut[1][uidm+kvolHalo*nu]*ut[0][uidn+kvolHalo*mu];
44 Sigma[0][i]+=a11*conj(ut[0][i+kvolHalo*nu])+a12*conj(ut[1][i+kvolHalo*nu]);
45 Sigma[1][i]+=-a11*ut[1][i+kvolHalo*nu]+a12*ut[0][i+kvolHalo*nu];
46 }
47 C_gather(ush[0], ut[0], kvol, id, nu);
48 C_gather(ush[1], ut[1], kvol, id, nu);
49#if(nproc>1)
50 CHalo_swap_dir(ush[0], 1, mu, DOWN); CHalo_swap_dir(ush[1], 1, mu, DOWN);
51#endif
52 //Next up, the -ν staple
53#pragma omp parallel for simd //aligned(ut[0],ut[1],ush[0],ush[1],Sigma[0],Sigma[1],iu,id:AVX)
54 for(int i=0;i<kvol;i++){
55 int uidm = iu[mu*kvol+i];
56 int didn = id[nu*kvol+i];
57 //uidm is correct here
58 Complex_f a11=conj(ush[0][uidm])*conj(ut[0][didn+kvolHalo*mu])-\
59 ush[1][uidm]*conj(ut[1][didn+kvolHalo*mu]);
60 Complex_f a12=-conj(ush[0][uidm])*ut[1][didn+kvolHalo*mu]-\
61 ush[1][uidm]*ut[0][didn+kvolHalo*mu];
62 Sigma[0][i]+=a11*ut[0][didn+kvolHalo*nu]-a12*conj(ut[1][didn+kvolHalo*nu]);
63 Sigma[1][i]+=a11*ut[1][didn+kvolHalo*nu]+a12*conj(ut[0][didn+kvolHalo*nu]);
64 }
65 }
66#pragma omp parallel for simd //aligned(ut[0],ut[1],Sigma[0],Sigma[1],dSdpi:AVX)
67 for(int i=0;i<kvol;i++){
68 const unsigned int ind = i+kvolHalo*mu;
69 Complex_f a11 = ut[0][ind]*Sigma[1][i]+ut[1][ind]*conj(Sigma[0][i]);
70 Complex_f a12 = ut[0][ind]*Sigma[0][i]+conj(ut[1][ind])*Sigma[1][i];
71
72 dSdpi[i+kvol*mu]=(double)(beta*cimag(a11));
73 dSdpi[i+kvol*(1*ndim+mu)]=(double)(beta*creal(a11));
74 dSdpi[i+kvol*(2*ndim+mu)]=(double)(beta*cimag(a12));
75 }
76 }
77 free(ush[0]); free(ush[1]); free(Sigma[0]); free(Sigma[1]);
78#endif
79 return 0;
80}
81void Force_s(double *dSdpi, Complex_f *ut[2], Complex_f *X1, Complex_f *X2, Complex_f gamval[20],\
82 unsigned int *iu, const unsigned short gamin[16],const float akappa, const unsigned short mu){
83#pragma omp parallel for simd
84 for(unsigned int i=0;i<kvol;i++){
85 const unsigned int ind=i+kvolHalo*mu;
86 const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
87 const unsigned int uid = iu[i+kvol*mu];
88 //Similarly to Hdslash we always see idirac*nc so we do that here too.
89 for(unsigned short idirac=0;idirac<nc*ndirac;idirac+=nc){
90 Complex_f X1s[nc]; Complex_f X1su[nc];
91 Complex_f X2s[nc]; Complex_f X2su[nc];
92
93 X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
94 X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
95 X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
96 X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
97
98 float dSdpis[3];
99 dSdpis[0]=dSdpi[i+kvol*mu];
100 //Multiplying by i and taking the real component is the same as taking the negative imaginary component
101 //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
102 //generators by hand
103 dSdpis[0]+=-akappa*cimag(
104 conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
105 +conj(X1s[1])*(u11s*X2su[0]+u12s*X2su[1])
106 +conj(X1su[0])*(u12s*X2s[0]-conj(u11s)*X2s[1])
107 +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1]));
108
109 dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
110 dSdpis[1]+=akappa*creal(
111 (conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
112 +conj(X1s[1])*(-u11s*X2su[0]-u12s*X2su[1])
113 +conj(X1su[0])*(-u12s*X2s[0]-conj(u11s)*X2s[1])
114 +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1])));
115
116 dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
117 dSdpis[2]+=-akappa*cimag(
118 conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
119 +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
120 +conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
121 +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1]));
122
123 const unsigned short gindex=mu*ndirac+(idirac>>1);
124 const Complex_f gamval_c=gamval[gindex];
125 //Rescaling gind by nc
126 const unsigned short gind = gamin[gindex]<<1;
127 X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
128 X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
129
130 //If you are asked to rederive the force from Montvay and Munster you'll notice that it should be kappa*gamma
131 //but below is only gamma. We rescaled gamma by kappa already when we defined it so that's where it has gone
132 dSdpis[0]+=-cimag(gamval_c*
133 (conj(X1s[0])* (-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
134 +conj(X1s[1])* (u11s *X2su[0]+u12s *X2su[1])
135 +conj(X1su[0])* (-u12s *X2s[0] +conj(u11s)*X2s[1])
136 +conj(X1su[1])*(u11s *X2s[0] +conj(u12s)*X2s[1])));
137 dSdpi[i+kvol*mu]=dSdpis[0];
138
139 dSdpis[1]+=creal(gamval_c*
140 (conj(X1s[0])* (-conj(u12s)*X2su[0] +conj(u11s)*X2su[1])
141 +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1])
142 +conj(X1su[0])* (u12s *X2s[0]+conj(u11s)*X2s[1])
143 +conj(X1su[1])* (-u11s *X2s[0]+conj(u12s)*X2s[1])));
144 dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
145
146 dSdpis[2]+=-cimag(gamval_c*
147 (conj(X1s[0])*(u11s *X2su[0]+u12s *X2su[1])
148 +conj(X1s[1])*(conj(u12s)*X2su[0]-conj(u11s)*X2su[1])
149 +conj(X1su[0])*(conj(u11s)*X2s[0]+u12s *X2s[1])
150 +conj(X1su[1])*(conj(u12s)*X2s[0]-u11s *X2s[1])));
151 dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
152 }
153 }
154 return;
155}
156void Force_t(double *dSdpi, Complex_f *ut[2],Complex_f *X1, Complex_f *X2, Complex_f gamval[20],\
157 float *dk[2], unsigned int *iu, const unsigned short gamin[16],float akappa){
158
159 const unsigned short mu=3;
160#pragma omp parallel for simd
161 for(unsigned int i=0;i<kvol;i++){
162 const unsigned int ind=i+kvolHalo*mu;
163 const Complex_f u11s=ut[0][ind]; const Complex_f u12s=ut[1][ind];
164 //TODO: The only diffrence with these is that the sign flips for the temporal components
165 // Can we figure out a way of doing this without having to read in a large array.
166 // Will result in a conditional inside a CUDA loop. If i>kvol3
167 const float dks[2] = {dk[0][i],dk[1][i]};
168 //Up indices
169 const unsigned int uid = iu[i+kvol*mu];
170 //Similarly to Hdslash we always see idirac*nc so we do that here too.
171 for(unsigned short idirac=0;idirac<ndirac*nc;idirac+=nc){
172 Complex_f X1s[nc]; Complex_f X1su[nc];
173 Complex_f X2s[nc]; Complex_f X2su[nc];
174
175 X1s[0]=X1[i+kvolHalo*(idirac)]; X1s[1]=X1[i+kvolHalo*(1+idirac)];
176 X1su[0]=X1[uid+kvolHalo*(idirac)]; X1su[1]=X1[uid+kvolHalo*(1+idirac)];
177 X2s[0]=2*X2[i+kvolHalo*(idirac)]; X2s[1]=2*X2[i+kvolHalo*(1+idirac)];
178 X2su[0]=2*X2[uid+kvolHalo*(idirac)]; X2su[1]=2*X2[uid+kvolHalo*(1+idirac)];
179
180 float dSdpis[3];
181 dSdpis[0]=dSdpi[i+kvol*mu];
182 //Multiplying by i and taking the real component is the same as taking the negative imaginary component
183 //The positions of u11 and u12 might look a bit funky here. That's just because we've multiplied by the
184 //generators by hand
185 dSdpis[0]+=-cimag(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
186 +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
187 +dks[1]*(conj(X1su[0])*(+u12s*X2s[0]-conj(u11s)*X2s[1])
188 +conj(X1su[1])*(-u11s*X2s[0]-conj(u12s)*X2s[1])));
189
190 dSdpis[1]=dSdpi[i+kvol*(ndim+mu)];
191 dSdpis[1]+=creal(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
192 +conj(X1s[1])*(-u11s *X2su[0]-u12s *X2su[1]))
193 +dks[1]*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
194 +conj(X1su[1])*( u11s *X2s[0]-conj(u12s)*X2s[1])));
195
196 dSdpis[2]=dSdpi[i+kvol*(2*ndim+mu)];
197 dSdpis[2]+=-cimag(dks[0]* (conj(X1s[0])* (u11s *X2su[0]+u12s *X2su[1])
198 +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
199 +dks[1]*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
200 +conj(X1su[1])* (-conj(u12s)*X2s[0]+u11s *X2s[1])));
201
202 const unsigned short gindex=mu*ndirac+(idirac>>1);
203 //Rescaling gind by nc
204 const unsigned short gind = gamin[gindex]<<1;
205 X2s[0]=2*X2[i+kvolHalo*(gind)]; X2s[1]=2*X2[i+kvolHalo*(1+gind)];
206 X2su[0]=2*X2[uid+kvolHalo*(gind)]; X2su[1]=2*X2[uid+kvolHalo*(1+gind)];
207
208 dSdpis[0]+=-cimag(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
209 +conj(X1s[1])*(u11s *X2su[0]+u12s *X2su[1]))
210 -dks[1]*(conj(X1su[0])* (u12s *X2s[0]-conj(u11s)*X2s[1])
211 +conj(X1su[1])*(-u11s *X2s[0]-conj(u12s)*X2s[1])));
212 dSdpi[i+kvol*mu]=dSdpis[0];
213
214 dSdpis[1]+=creal(dks[0]*(conj(X1s[0])*(-conj(u12s)*X2su[0]+conj(u11s)*X2su[1])
215 +conj(X1s[1])*(-u11s*X2su[0]-u12s *X2su[1]))
216 -dks[1]*(conj(X1su[0])*(-u12s *X2s[0]-conj(u11s)*X2s[1])
217 +conj(X1su[1])*(u11s*X2s[0]-conj(u12s)*X2s[1])));
218 dSdpi[i+kvol*(ndim+mu)]=dSdpis[1];
219
220 dSdpis[2]+=-cimag(dks[0]*(conj(X1s[0])*(u11s*X2su[0] +u12s *X2su[1])
221 +conj(X1s[1])* (conj(u12s)*X2su[0]-conj(u11s)*X2su[1]))
222 -dks[1]*(conj(X1su[0])*(-conj(u11s)*X2s[0]-u12s *X2s[1])
223 +conj(X1su[1])*(-conj(u12s)*X2s[0]+u11s *X2s[1])));
224 dSdpi[i+kvol*(2*ndim+mu)]=dSdpis[2];
225 }
226 }
227}
228int Force(double *dSdpi, const bool iflag, double res1, Complex *X0, Complex *X1, Complex *Phi,\
229 Complex *ut[2], Complex_f *ut_f[2],unsigned int *iu,unsigned int *id,\
230 Complex gamval[20],Complex_f gamval_f[20],const unsigned short gamin[16],Complex *sigval,Complex_f *sigval_f, unsigned short *sigin,\
231 double *dk[2], float *dk_f[2],const Complex_f jqq, const float akappa,const float beta,const float c_sw,double *ancg){
232 const char funcname[] = "Force";
233#ifdef USE_GPU
234 int device=-1;
235 cudaGetDevice(&device);
236#endif
237#ifndef NO_GAUGE
238 Gauge_force(dSdpi,ut_f,iu,id,beta);
239#else
240#endif
243 if(!akappa)
244 return 0;
245 //X1=(M†M)^{1} Phi
246 int itercg=1;
247 Complex_f *clover[2];
248#ifdef USE_GPU
249 Complex_f *X1_f, *X2_f;
250 cudaMallocAsync((void **)&X1_f,kferm2Halo*sizeof(Complex_f),streams[1]);
251 cudaMallocAsync((void **)&X2_f,kferm2Halo*sizeof(Complex_f),streams[0]);
252#else
253 Complex_f *X1_f= (Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
254 Complex_f *X2_f= (Complex_f *)aligned_alloc(AVX,kferm2Halo*sizeof(Complex_f));
255#endif
256 if(c_sw)
257 Clover(clover,ut_f,iu,id);
258
259 for(int na = 0; na<nf; na++){
260#ifdef USE_GPU
261#if(nproc>1) //Strided
262 for(unsigned short j=0;j<nc*idirac;j++)
263 cudaMemcpyAsync(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[j]);
264#else
265 cudaMemcpyAsync(X1,X0+na*kferm2,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,NULL);
266#endif
267#else
268 for(unsigned short j=0;j<nc*ndirac;j++)
269 memcpy(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex));
270#endif
271 if(!iflag){
272 int itercg=1;
273#ifdef USE_GPU
274 Complex *smallPhi;
275 cudaMallocAsync((void **)&smallPhi,kferm2*sizeof(Complex),streams[0]);
276#else
277 Complex *smallPhi = (Complex *)aligned_alloc(AVX,kferm2*sizeof(Complex));
278#endif
279 Fill_Small_Phi(na, smallPhi, Phi);
281 Congradq(na,res1,X1,smallPhi,ut,ut_f,clover,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
282 jqq,akappa,c_sw,&itercg);
283#ifdef USE_GPU
284 cudaFreeAsync(smallPhi,streams[0]);
285#else
286 free(smallPhi);
287#endif
288 *ancg+=itercg;
289#ifdef USE_GPU
290 alignas(16) const Complex blasa=2.0; alignas(16) const double blasb=-1.0;
291 cublasZdscal(cublas_handle,kferm2,&blasb,(cuDoubleComplex *)(X0+na*kferm2),1);
292#if(nproc>1) //strided
293 for(unsigned short j=0;j<nc*ndirac;j++)
294 cublasZaxpy(cublas_handle,kvol,(cuDoubleComplex *)&blasa,(cuDoubleComplex *)X1+j*kvolHalo,1,(cuDoubleComplex *)X0+na*kferm2+j*kvol,1);
295#else
296 cublasZaxpy(cublas_handle,kferm2,(cuDoubleComplex *)&blasa,(cuDoubleComplex *)X1,1,(cuDoubleComplex *)(X0+na*kferm2),1);
297#endif
298#elif (defined __USE_MKL__||defined OPENBLAS||defined AMD_BLAS)
299 const Complex blasa=2.0; const Complex blasb=-1.0;
300 //This is not a general BLAS Routine. BLIS and MKl support it
301 //CUDA and GSL does not support it
302 cblas_zaxpby(kferm2, &blasa, X1, 1, &blasb, X0+na*kferm2, 1);
303#elifdef USE_BLAS
304 const Complex blasa=2.0; const double blasb=-1.0;
305 cblas_zdscal(kferm2,blasb,X0+na*kferm2,1);
306 for(unsigned short j=0;j<nc*ndirac;j++)
307 cblas_zaxpy(kvol,&blasa,X1+j*kvolHalo,1,X0+na*kferm2+j*kvol,1);
308#else
309#pragma omp parallel for simd collapse(2) aligned(X0,X1:AVX)
310 for(int idirac=0;idirac<ndirac;idirac++){
311 for(int i=0;i<kvol;i++)
312 X0[i+kvol*(0+nc*(idirac+ndirac*na))]=
313 2*X1[i+kvolHalo*(0+idirac*c)]-X0[i+kvol*(0+nc*(idirac+ndirac*na))];
314 X0[i+kvol*(1+nc*(idirac+ndirac*na))]=
315 2*X1[i+kvolHalo*(1+idirac*c)]-X0[i+kvol*(1+nc*(idirac+ndirac*na))];
316 }
317#endif
318 }
319#ifdef USE_GPU
320 cudaDeviceSynchronise();
321#endif
322 //Since it has to be stridded in MPI, we have to pass kvol and nc*ndirac instead of kferm2
323 ComplexConvert(X1_f,X1,kvol,true,nc*ndirac);
324 Hdslash_f(X2_f,X1_f,ut_f,iu,id,gamval_f,gamin,dk_f,akappa);
325 if(c_sw)
326 HbyClover_f(X2_f,X1_f,clover,sigval_f,akappa,sigin,false);
327 //NOTE: This was orginally two. But was changed as a test for Claude so the two appears inside Force_s and Force_t
328 //It may need to be reverted back later
329 alignas(8) const float blasd=1.0;
330#ifdef USE_GPU
331 cudaDeviceSynchronise();
332#if(nproc>1)
333 for(unsigned short j=0;j<nc*ndirac;j++)
334 cublasCsscal(cublas_handle,kvol, &blasd, (cuComplex *)X2_f+j*kvolHalo, 1);
335#else
336 cublasCsscal(cublas_handle,kferm2, &blasd, (cuComplex *)X2_f, 1);
337#endif
338#elif defined USE_BLAS
339 for(unsigned short j=0;j<nc*ndirac;j++)
340 cblas_csscal(kvol, blasd, X2_f+j*kvolHalo, 1);
341#else
342#pragma unroll
343#pragma omp parallel for simd collapse(2) aligned(X2_f:AVX)
344 for(unsigned short j=0;j<nc*ndirac;j++)
345 for(unsigned int i=0;i<kvol;i++)
346 X2_f[i+j*kvolHalo]*=1;
347#endif
348#if(npx>1)
349 CHalo_swap_dir(X1_f,nc*ndirac,0,DOWN);
350 CHalo_swap_dir(X2_f,nc*ndirac,0,DOWN);
351#endif
352#if(npy>1)
353 CHalo_swap_dir(X1_f,nc*ndirac,1,DOWN);
354 CHalo_swap_dir(X2_f,nc*ndirac,1,DOWN);
355#endif
356#if(npz>1)
357 CHalo_swap_dir(X1_f,nc*ndirac,2,DOWN);
358 CHalo_swap_dir(X2_f,nc*ndirac,2,DOWN);
359#endif
360#if(npt>1)
361 CHalo_swap_dir(X1_f,nc*ndirac,3,DOWN);
362 CHalo_swap_dir(X2_f,nc*ndirac,3,DOWN);
363#endif
364
365 // The original FORTRAN Comment:
366 // dSdpi=dSdpi-Re(X1*(d(Mdagger)dp)*X2) -- Yikes!
367 // we're gonna need drugs for this one......
368 //
369 // Makes references to X1(.,.,iu(i,mu)) AND X2(.,.,iu(i,mu))
370 // as a result, need to swap the DOWN halos in all dirs for
371 // both these arrays, each of which has 8 cpts
372 //
373#ifdef USE_GPU
374 cuForce(dSdpi,ut_f,X1_f,X2_f,gamval_f,dk_f,iu,gamin,akappa,dimGrid,dimBlock);
375 cudaDeviceSynchronise();
376#else
377 //Thankfully the CUDA version is much neater so we're using that style going forwards
378 for(unsigned short mu=0;mu<ndim-1;mu++)
379 Force_s(dSdpi,ut_f,X1_f,X2_f,gamval_f,iu,gamin,akappa,mu);
380 Force_t(dSdpi,ut_f,X1_f,X2_f,gamval_f,dk_f,iu,gamin,akappa);
381#endif
382 if(c_sw){
383 //Clover_Force(dSdpi,ut_f,X1_f,X2_f,sigval_f,sigin,iu,id,akappa);
384 Clov_Force(dSdpi,ut_f,X1_f,X2_f,sigval_f,sigin,iu,id,akappa);
385 Clover_free(clover);
386 }
387 }
388#ifdef USE_GPU
389 cudaFreeAsync(X1_f,streams[0]); cudaFreeAsync(X2_f,streams[1]);
390#else
391 free(X1_f); free(X2_f);
392#endif
393 return 0;
394}
clover.h
Routines needed for Clover improved wilson fermions.
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:542
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:382
Clover_free
void Clover_free(Complex_f *clover[nc])
Free's memory used for clover terms and leaves.
Definition clover.c:750
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
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
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
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
C_gather
int C_gather(Complex_f *x, Complex_f *y, int n, unsigned int *table, unsigned int mu)
Extracts all the single precision gauge links in the direction only.
Definition su2hmc.c:291
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
cuForce
void cuForce(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, dim3 dimGrid, dim3 dimBlock)
Calculates the force at each intermediate time.
Definition cuforce.cu:373
cuGauge_force
void cuGauge_force(Complex_f *ut[2], double *dSdpi, float beta, unsigned int *iu, unsigned int *id, dim3 dimGrid, dim3 dimBlock)
Calculate the gauge contribution to the force.
Definition cuforce.cu:319
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
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
CHalo_swap_dir
int CHalo_swap_dir(Complex_f *c, int ncpt, int idir, int layer)
Swaps the halos along the axis given by idir in the direction given by layer.
matrices.h
Matrix multiplication and related declarations.
DOWN
#define DOWN
Flag for send down.
Definition par_mpi.h:37
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
kferm2Halo
#define kferm2Halo
Dirac lattice and halo.
Definition sizes.h:238
kvol
#define kvol
Sublattice volume.
Definition sizes.h:163
Complex
#define Complex
Double precision complex number.
Definition sizes.h:64
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
cublas_handle
cublasHandle_t cublas_handle
Handle for cuBLAS.
Definition main.c:47
Complex_f
#define Complex_f
Single precision complex number.
Definition sizes.h:62
dimGrid
dim3 dimGrid
Default grid size. First component is normally nt. Second and third depend whatever is needed to get ...
Definition cusu2hmc.cu:27
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
dimBlock
dim3 dimBlock
Default block size. Usually 128.
Definition cusu2hmc.cu:25
streams
cudaStream_t streams[ndirac *ndim *nadj]
An array of concurrent GPU streams to keep it busy.
Definition cusu2hmc.cu:29
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