su2hmc
Loading...
Searching...
No Matches
fermionic.c
Go to the documentation of this file.
1
5#include <matrices.h>
6#include <clover.h>
7int Measure(double *pbp, double *endenf, double *denf, Complex *qq, Complex *qbqb, double res, int *itercg,\
8 Complex *ut[2], Complex_f *ut_f[2], unsigned int *iu, unsigned int *id,\
9 Complex gamval[20], Complex_f gamval_f[20], const unsigned short gamin[16],\
10 Complex *sigval,Complex_f *sigval_f, unsigned short *sigin, double *dk[2],float *dk_f[2],\
11 Complex_f jqq, float akappa, float c_sw,Complex *Phi){
12 const char funcname[] = "Measure";
13 //This x is just a storage container
14
15#ifdef USE_GPU
16 int device=-1;
17 cudaGetDevice(&device);
18 Complex *x, *xi, *R1 ; Complex_f *xi_f, *R1_f, *clover[nc];
19#ifdef _DEBUG
20 cudaMallocManaged((void **)&R1,kferm*sizeof(Complex), cudaMemAttachGlobal);
21 cudaMallocManaged((void **)&R1_f,kferm*sizeof(Complex_f), cudaMemAttachGlobal);
22#else
23 cudaMallocAsync((void **)&R1,kferm*sizeof(Complex),streams[0]);
24 cudaMallocAsync((void **)&R1_f,kferm*sizeof(Complex_f),streams[0]);
25#endif
26 cudaMallocManaged((void **)&x,kfermHalo*sizeof(Complex), cudaMemAttachGlobal);
27 cudaMallocManaged((void **)&xi,kferm*sizeof(Complex), cudaMemAttachGlobal);
28 cudaMallocManaged((void **)&xi_f,kfermHalo*sizeof(Complex_f), cudaMemAttachGlobal);
29#else
30 Complex_f *clover[nc];
31 Complex *x =(Complex *)aligned_alloc(AVX,kfermHalo*sizeof(Complex));
32 Complex *xi =(Complex *)aligned_alloc(AVX,kferm*sizeof(Complex));
33 Complex_f *xi_f =(Complex_f *)aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
34 Complex_f *R1_f = (Complex_f *)aligned_alloc(AVX,kferm*sizeof(Complex_f));
35 Complex *R1 = (Complex *)aligned_alloc(AVX,kferm*sizeof(Complex));
36#endif
37 //Setting up noise. Again need that annoying stride
38 for(unsigned short j=0;j<nc*ngorkov;j++){
39 Gauss_c(xi_f+j*kvolHalo, kvol, 0, (float)(1/sqrt(2)));
40 ComplexConvert(xi_f+j*kvolHalo,xi+j*kvol,kvol,false,1);
41 }
42#ifdef USE_GPU
43#if (nproc>1) //strided
44 for(unsigned short j=0;j<nc*ngorkov;j++){
45 cudaMemcpyAsync(x+j*kvolHalo, xi+j*kvol, kvol*sizeof(Complex),cudaMemcpyDefault,0);
46 }
47#else
48 cudaMemcpyAsync(x, xi, kferm*sizeof(Complex),cudaMemcpyDefault,0);
49#endif
50#else
51#pragma omp parallel for
52 for(unsigned short j=0;j<nc*ngorkov;j++){
53 memcpy(x+j*kvolHalo, xi+j*kvol, kvol*sizeof(Complex));
54 }
55#endif
56 //R_1= @f$M^\dagger\Xi@f$
57 //global
58 Dslashd_f(R1_f,xi_f,ut_f,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
59 if(c_sw){
60 Clover(clover,ut_f,iu,id);
61 ByClover_f(R1_f,xi_f,clover,sigval_f,akappa,sigin,true);
62 }
63 ComplexConvert(R1_f,R1,kferm,false,1);
64#ifdef USE_GPU
65 cudaFree(xi_f);
66 cudaDeviceSynchronise();
67#ifdef _DEBUG
68 cudaFree(R1_f);
69#else
70 cudaFreeAsync(R1_f,streams[0]);
71#endif
72 cudaMemcpyAsync(Phi, R1, kferm*sizeof(Complex),cudaMemcpyDefault,streams[0]);
73 cudaDeviceSynchronise();
74#else
75 free(xi_f); free(R1_f);
76 memcpy(Phi, R1, kferm*sizeof(Complex));
77#endif
79 if(Congradp(0, res, Phi,R1,ut,ut_f,clover,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,c_sw,itercg)==ITERLIM){
80 //Clean exit
81#ifdef USE_GPU
82#ifdef _DEBUG
83 if(c_sw){
84 cudaFree(clover[0]); cudaFree(clover[1]);
85 }
86 cudaFree(R1);
87#else
88 if(c_sw){
89 cudaFreeAsync(clover[0],streams[1]); cudaFreeAsync(clover[1],streams[2]);
90 }
91 cudaFreeAsync(R1,streams[0]);
92#endif
93 cudaFree(x); cudaFree(xi);
94#else
95 if(c_sw){
96 free(clover[0]); free(clover[1]);
97 }
98 free(x); free(xi); free(R1);
99 return ITERLIM;
100#endif
101 }
102#ifdef USE_GPU
103 cudaMemcpyAsync(xi,R1,kferm*sizeof(Complex),cudaMemcpyDefault,streams[0]);
104#ifdef _DEBUG
105 if(c_sw){
106 cudaFree(clover[0]); cudaFree(clover[1]);
107 }
108 cudaFree(R1);
109#else
110 if(c_sw){
111 cudaFreeAsync(clover[0],streams[1]); cudaFreeAsync(clover[1],streams[2]);
112 }
113 cudaDeviceSynchronise();
114 cudaFreeAsync(R1,streams[0]);
115#endif
116 cudaDeviceSynchronise();
117#else
118 memcpy(xi,R1,kferm*sizeof(Complex));
119 if(c_sw){
120 free(clover[0]); free(clover[1]);
121 }
122 free(R1);
123#endif
124 *pbp = 0;
125#ifdef USE_BLAS
126 alignas(16) Complex buff;
127#ifdef USE_GPU
128#if(nproc>1)
129 for(unsigned short j=0;j<ngorkov*nc;j++){
130 buff=0;
131 cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)x+j*kvolHalo,1,(cuDoubleComplex *)xi+j*kvol,1,(cuDoubleComplex *)&buff);
132 *pbp+=creal(buff);
133 }
134#else
135 cublasZdotc(cublas_handle,kferm,(cuDoubleComplex *)x,1,(cuDoubleComplex *)xi,1,(cuDoubleComplex *)&buff);
136 *pbp+=creal(buff);
137#endif
138 cudaDeviceSynchronise();
139#elif defined USE_BLAS
140 for(unsigned short j=0;j<ngorkov*nc;j++){
141 buff=0;
142 cblas_zdotc_sub(kvol, x+j*kvolHalo, 1, xi+j*kvol, 1, &buff);
143 *pbp+=creal(buff);
144 }
145#endif
146#else
147#pragma unroll
148 for(int i=0;i<kferm;i++)
149 *pbp+=creal(conj(x[i])*xi[i]);
150#endif
151#if(nproc>1)
152 Par_dsum(pbp);
153#endif
154 *pbp/=4*gvol;
155
156 *qbqb=*qq=0;
157#if defined USE_BLAS
158 for(int idirac = 0; idirac<ndirac; idirac++){
159 int igork=idirac+4;
160 //Unrolling the colour indices, Then its just (γ_5*x)*Ξ or (γ_5*Ξ)*x
161#pragma unroll
162 for(int ic = 0; ic<nc; ic++){
163 alignas(16) Complex dot=0;
164 //Because we have kvol on the outer index and are summing over it, we set the
165 //step for BLAS to be ngorkov*nc=16.
166 //Does this make sense to do on the GPU?
167#ifdef USE_GPU
168 cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)x+kvolHalo*(idirac*nc+ic),1,(cuDoubleComplex *)xi+kvol*(igork*nc+ic), 1,(cuDoubleComplex *)&dot);
169#else
170 cblas_zdotc_sub(kvol, x+kvolHalo*(idirac*nc+ic), 1, xi+kvol*(igork*nc+ic), 1, &dot);
171#endif
172 *qbqb+=gamval[4*ndirac+idirac]*dot;
173#ifdef USE_GPU
174 cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)x+kvolHalo*(igork*nc+ic),1,(cuDoubleComplex *)xi+kvol*(idirac*nc+ic), 1,(cuDoubleComplex *)&dot);
175#else
176 cblas_zdotc_sub(kvol, x+kvolHalo*(igork*nc+ic), 1, xi+kvol*(idirac*nc+ic), 1, &dot);
177#endif
178 *qq-=gamval[4*ndirac+idirac]*dot;
179 }
180 }
181#else
182 //What is the optimal order to evaluate these in?
183#pragma omp parallel for simd collapse(2) aligned(x,xi:AVX) reduction(+:*qq,*qbqb)
184 for(int idirac = 0; idirac<ndirac; idirac++)
185 for(int i=0; i<kvol; i++){
186 int igork=idirac+4;
187 *qbqb+=gamval[4*ndirac+idirac]*conj(x[i+kvolHalo*(idirac*nc)])*xi[i+kvol*(igork*nc)];
188 *qbqb+=gamval[4*ndirac+idirac]*conj(x[i+kvolHalo*(idirac*nc+1)])*xi[i+kvol*(igork*nc+1)];
189 *qq-=gamval[4*ndirac+idirac]*conj(x[i+kvolHalo*(igork*nc)])*xi[i+kvol*(idirac*nc)];
190 *qq-=gamval[4*ndirac+idirac]*conj(x[i+kvolHalo*(igork*nc+1)])*xi[i+kvol*(idirac*nc+1)];
191 }
192#endif
193 //In the FORTRAN Code dsum was used instead despite qq and qbqb being complex
194 //Since we only care about the real part this shouldn't cause (m)any serious issues
195#if(nproc>1)
196 Par_dsum((double *)qq); Par_dsum((double *)qbqb);
197#endif
198 *qq=(*qq+*qbqb)/(2*gvol);
199 Complex xu, xd, xuu, xdd;
200 xu=xd=xuu=xdd=0;
201
202 //Halos
203#if(npt>1)
204 ZHalo_swap_dir(x,16,3,DOWN); ZHalo_swap_dir(x,16,3,UP);
205#endif
206 //Pesky halo exchange indices again
207 //The halo exchange for the trial fields was done already at the end of the trajectory
208 //No point doing it again
209
210 //Instead of typing id[i+kvol*3] a lot, we'll just assign them to variables.
211 //Idea. One loop instead of two loops but for xuu and xdd just use ngorkov-(igorkov+1) instead
212 //Dirty CUDA work around since it won't convert thrust<complex> to double
213
214 //TODO: Make the code below CUDA friendly.
215 for(unsigned short igorkov=0; igorkov<4; igorkov++){
216 const unsigned short igork1=gamin[3*ndirac+igorkov];
217 //For the C Version I'll try and factorise where possible
218#pragma omp parallel for simd aligned(dk,x,xi:AVX) reduction(+:xu)
219 for(unsigned int i = 0; i<kvol; i++){
220 unsigned int did=id[3*kvol+i];
221 xu+=dk[1][did]*(conj(x[did+kvolHalo*(igorkov*nc)])*(\
222 ut[0][did+kvol*3]*(xi[i+kvol*(igork1)*nc]-xi[i+kvol*(igorkov)*nc])+\
223 ut[1][did+kvol*3]*(xi[i+kvol*(igork1)*nc+1]-xi[i+kvol*(igorkov)*nc+1]) )+\
224 conj(x[did+kvolHalo*(igorkov*nc+1)])*(\
225 conj(ut[0][did+kvol*3])*(xi[i+kvol*(igork1)*nc+1]-xi[i+kvol*(igorkov)*nc+1])+\
226 conj(ut[1][did+kvol*3])*(xi[i+kvol*(igorkov)*nc]-xi[i+kvol*(igork1)*nc])));
227 }
228 }
229 for(unsigned short igorkov=0; igorkov<4; igorkov++){
230 const unsigned short igork1=gamin[3*ndirac+igorkov];
231#pragma omp parallel for simd aligned(dk,x,xi:AVX) reduction(+:xd)
232 for(unsigned int i = 0; i<kvol; i++){
233 unsigned int uid=iu[3*kvol+i];
234 xd+=dk[0][i]*(conj(x[uid+kvolHalo*(igorkov*nc)])*(\
235 conj(ut[0][i+kvol*3])*(xi[i+kvol*(igork1*nc)]+xi[i+kvol*(igorkov*nc)])-\
236 ut[1][i+kvol*3]*(xi[i+kvol*(igork1*nc+1)]+xi[i+kvol*(igorkov*nc+1)]) )+\
237 conj(x[uid+kvolHalo*(igorkov*nc+1)])*(\
238 ut[0][i+kvol*3]*(xi[i+kvol*(igork1*nc+1)]+xi[i+kvol*(igorkov*nc+1)])+\
239 conj(ut[1][i+kvol*3])*(xi[i+kvol*(igorkov*nc)]+xi[i+kvol*(igork1*nc)]) ) );
240 }
241 }
242 for(unsigned short igorkovPP=4; igorkovPP<8; igorkovPP++){
243 const unsigned short igork1PP=4+gamin[3*ndirac+igorkovPP-4];
244#pragma omp parallel for simd aligned(dk,x,xi:AVX) reduction(+:xuu)
245 for(unsigned int i = 0; i<kvol; i++){
246 unsigned int did=id[3*kvol+i];
247 xuu-=dk[0][did]*(conj(x[did+kvolHalo*(igorkovPP*nc)])*(\
248 ut[0][did+kvol*3]*(xi[i+kvol*(igork1PP*nc)]-xi[i+kvol*(igorkovPP*nc)])+\
249 ut[1][did+kvol*3]*(xi[i+kvol*(igork1PP*nc+1)]-xi[i+kvol*(igorkovPP*nc+1)]) )+\
250 conj(x[did+kvolHalo*(igorkovPP*nc+1)])*(\
251 conj(ut[0][did+kvol*3])*(xi[i+kvol*(igork1PP)*nc+1]-xi[i+kvol*(igorkovPP)*nc+1])+\
252 conj(ut[1][did+kvol*3])*(xi[i+kvol*(igorkovPP)*nc]-xi[i+kvol*(igork1PP)*nc]) ) );
253 }
254 }
255 for(unsigned short igorkovPP=4; igorkovPP<8; igorkovPP++){
256 const unsigned short igork1PP=4+gamin[3*ndirac+igorkovPP-4];
257#pragma omp parallel for simd aligned(dk,x,xi:AVX) reduction(+:xdd)
258 for(unsigned int i = 0; i<kvol; i++){
259 unsigned int uid=iu[3*kvol+i];
260 xdd-=dk[1][i]*(conj(x[uid+kvolHalo*(igorkovPP*nc)])*(\
261 conj(ut[0][i+kvol*3])*(xi[i+kvol*(igork1PP*nc)]+xi[i+kvol*(igorkovPP*nc)])-\
262 ut[1][i+kvol*3]*(xi[i+kvol*(igork1PP*nc+1)]+xi[i+kvol*(igorkovPP*nc+1)]) )+\
263 conj(x[uid+kvolHalo*(igorkovPP*nc+1)])*(\
264 ut[0][i+kvol*3]*(xi[i+kvol*(igork1PP*nc+1)]+xi[i+kvol*(igorkovPP*nc+1)])+\
265 conj(ut[1][i+kvol*3])*(xi[i+kvol*(igorkovPP*nc)]+xi[i+kvol*(igork1PP*nc)]) ) );
266 }
267 }
268 *endenf=creal(xu-xd-xuu+xdd);
269 *denf=creal(xu+xd+xuu+xdd);
270
271#if(nproc>1)
272 Par_dsum(endenf); Par_dsum(denf);
273#endif
274 *endenf/=2*gvol; *denf/=2*gvol;
275 //Future task. Chiral susceptibility measurements
276#ifdef USE_GPU
277 cudaFree(x); cudaFree(xi);
278 //Revert index and gauge arrays
279 // Transpose_z(ut[0],ndim,kvol);
280 // Transpose_z(ut[1],ndim,kvol);
281 //Transpose_U(iu,ndim,kvol);
282 //Transpose_U(id,ndim,kvol);
283#else
284 free(x); free(xi);
285#endif
286 return 0;
287}
clover.h
Routines needed for Clover improved wilson fermions.
ITERLIM
#define ITERLIM
Exceeded max number of iterations.
Definition errorcodes.h:137
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:336
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
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
Measure
int Measure(double *pbp, double *endenf, double *denf, Complex *qq, Complex *qbqb, double res, int *itercg, 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], Complex_f jqq, float akappa, float c_sw, Complex *Phi)
Calculate fermion expectation values via a noisy estimator.
Definition fermionic.c:7
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
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
ZHalo_swap_dir
int ZHalo_swap_dir(Complex *z, int ncpt, int idir, int layer)
Swaps the halos along the axis given by idir in the direction given by layer.
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.
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
matrices.h
Matrix multiplication and related declarations.
UP
#define UP
Flag for send up.
Definition par_mpi.h:39
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
ngorkov
#define ngorkov
Gor'kov indices.
Definition sizes.h:190
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
ndirac
#define ndirac
Dirac indices.
Definition sizes.h:186
cudaDeviceSynchronise
#define cudaDeviceSynchronise()
Get rid of that bastardised yankee English.
Definition sizes.h:53
gvol
#define gvol
Lattice volume.
Definition sizes.h:98
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
kvolHalo
#define kvolHalo
Subvolume + halo size.
Definition sizes.h:234
kfermHalo
#define kfermHalo
Gor'kov lattice and halo.
Definition sizes.h:236
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