su2hmc
Loading...
Searching...
No Matches
Fermionic Observables
Observable Quantities
Collaboration diagram for Fermionic Observables:

Functions

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.

Detailed Description

Function Documentation

◆ 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.

Matrix inversion via conjugate gradient algorithm Solves \(MX=X_1\) (Numerical Recipes section 2.10 pp.70-73)
uses NEW lookup tables ** Implemented in Congradp()

Parameters
[out]pbp\(\langle\bar{\Psi}\Psi\rangle\)
[out]endenfEnergy density
[out]denfNumber Density
[out]qqDiquark condensate
[out]qbqbAntidiquark condensate
[in]resConjugate Gradient Residue
[in]itercgIterations of Conjugate Gradient
[in]ut,ut_fDouble/float precision gauge field
[in]iu,idUp/down Lattice indices
[in]gamval,gamval_fDouble/float precision gamma matrices rescaled by kappa
[in]gaminIndices for Dirac terms
[in]sigval,sigval_fDouble/float Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}/2\)
[in]siginWhat element of the spinor is multiplied by row idirac each sigma matrix?
[in]dk,dk_fDouble/float \(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
[in]jqqDiquark source
[in]akappaHopping parameter
[in]c_swClover parameter
[in]PhiPseudofermion field
Returns
Zero on success, integer error code otherwise
Postcondition
The values of Phi are not used. Since the memory is allocated already it is instead overwritten with the noisy estimator.

Evaluate xi = (M^† M)^-1 R_1

Definition at line 7 of file fermionic.c.

11 {
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}
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
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
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

References AVX, ByClover_f(), Clover(), Complex, Complex_f, ComplexConvert(), Congradp(), conj(), creal, cublas_handle, cudaDeviceSynchronise, DOWN, Dslashd_f(), Gauss_c(), gvol, ITERLIM, kferm, kfermHalo, kvol, kvolHalo, nc, ndirac, ngorkov, Par_dsum(), streams, UP, and ZHalo_swap_dir().

Here is the call graph for this function:
Here is the caller graph for this function: