Observable Quantities

Collaboration diagram for Observable Quantities:

Topics
	Bosonic Observables
	Fermionic Observables

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

Detailed Description

Function Documentation

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.

Parameters

[out]	h	Hamiltonian
[out]	s	Action
[in]	res2	Limit for conjugate gradient
[in]	pp	Momentum field
[in]	X0	Up/down partitioned pseudofermion field
[in]	X1	Holder for the partitioned fermion field, then the conjugate gradient output
[in]	Phi	Pseudofermion field
[in]	ut,ud	Gauge fields (single/double precision)
[in]	iu,id	Lattice indices
[in]	gamval,gamval_f	Gamma matrices rescaled by kappa
[in]	gamin	Gamma indices
[in]	sigval,sigval_f	Commutators of gamma matrices scaled by \(\frac{c_\text{SW}}{2}\)
[in]	sigin	What element of the spinor is multiplied by row idirac each sigma matrix?
[in]	dk,dk_f	\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
[in]	jqq	Diquark source
[in]	akappa	Hopping parameter
[in]	beta	Inverse gauge coupling
[in]	c_sw	Clover coefficient. If non-zero calculate the clover contribution
[in]	ancgh	Conjugate gradient iterations counter
[in]	traj	Calling trajectory for error reporting

Returns

Zero on success. Integer Error code otherwise.

Postcondition

h and s overwritten with output

Definition at line 168 of file su2hmc.c.

               {
   const char funcname[] = "Hamilton";
   //Iterate over momentum terms.
#ifdef USE_GPU
   double hp;
   int device=-1;
   cudaGetDevice(&device);
   //cudaMemPrefetchAsync(pp,kmom*sizeof(double),device,NULL);
   cublasDnrm2(cublas_handle, kmom, pp, 1,&hp);
   hp*=hp;
#elif defined USE_BLAS
   double hp = cblas_dnrm2(kmom, pp, 1);
   hp*=hp;
#else
   double hp=0;
   for(unsigned int i = 0; i<kmom; i++)
      hp+=pp[i]*pp[i]; 
#endif
   hp*=0.5;
   double avplaqs, avplaqt;
   double hg = 0;
   //avplaq? isn't seen again here.
   Average_Plaquette(&hg,&avplaqs,&avplaqt,ut,iu,beta);
 
   alignas(8) double hf = 0; int itercg = 0;
#ifdef USE_GPU
   Complex *smallPhi;
#ifdef _DEBUG
   cudaMallocManaged((void **)&smallPhi,kferm2*sizeof(Complex),cudaMemAttachGlobal);
#else
   cudaMallocAsync((void **)&smallPhi,kferm2*sizeof(Complex),NULL);
#endif
#else
   Complex *smallPhi = aligned_alloc(AVX,kferm2*sizeof(Complex));
#endif
   Complex_f *clover[nc];
   if(c_sw)
      Clover(clover,ut,iu,id);
   //Iterating over flavours
   for(unsigned short na=0;na<nf;na++){
#ifdef USE_GPU
#if (nproc>1) //strided for multi-GPU
      for(unsigned short j=0;j<nc*ndirac;j++)
         cudaMemcpyAsync(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[j]);
#else
      cudaMemcpyAsync(X1,X0+na*kferm2,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[0]);
#endif
#else
      for(unsigned short j=0;j<nc*ndirac;j++)
         memcpy(X1+j*kvolHalo,X0+na*kferm2+j*kvol,kvol*sizeof(Complex));
#endif
      Fill_Small_Phi(na, smallPhi, Phi);
      if(Congradq(na,res2,X1,smallPhi,ud,ut,clover,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
               jqq,akappa,c_sw,&itercg))
         fprintf(stderr,"Trajectory %d\n", traj);
 
      *ancgh+=itercg;
#ifdef USE_GPU
#if (nproc>1) //strided for multi-GPU
      for(unsigned short j=0;j<nc*ndirac;j++)
         cudaMemcpyAsync(X0+na*kferm2+j*kvol,X1+j*kvolHalo,kvol*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[j]);
#else
      cudaMemcpyAsync(X0+na*kferm2,X1,kferm2*sizeof(Complex),cudaMemcpyDeviceToDevice,streams[0]);
#endif
#else
      for(unsigned short j=0;j<nc*ndirac;j++)
         memcpy(X0+na*kferm2+j*kvol,X1+j*kvolHalo,kvol*sizeof(Complex));
#endif
      Fill_Small_Phi(na, smallPhi,Phi);
#ifdef USE_GPU
      alignas(16) Complex dot=0;
#if (nproc>1)
      for(unsigned short j=0;j<nc*ndirac;j++){
         alignas(16) Complex buff;
         cublasZdotc(cublas_handle,kvol,(cuDoubleComplex *)smallPhi+j*kvol,1,(cuDoubleComplex *) X1+j*kvolHalo,1,(cuDoubleComplex *) &buff);
         dot+=buff;
      }
#else
      cublasZdotc(cublas_handle,kferm2,(cuDoubleComplex *)smallPhi,1,(cuDoubleComplex *) X1,1,(cuDoubleComplex *) &dot);
#endif
      hf+=creal(dot);
#elif defined USE_BLAS
      Complex dot=0;
      for(unsigned short j=0;j<nc*ndirac;j++){
         alignas(16) Complex buff=0;
         cblas_zdotc_sub(kvol, smallPhi+j*kvol, 1, X1+j*kvolHalo, 1, &buff);
         dot+=buff;
      }
      hf+=creal(dot);
#else
      //It is a dot product of the flattened arrays, could use
      //a module to convert index to coordinate array...
#pragma omp parallel for simd collapse(2) aligned(smallPhi,X1:AVX)
      for(unsigned short j=0;j<nc*ndirac;j++)
         for(unsigned int i=0;i<kvol;i++)
            hf+=creal(conj(smallPhi[i+j*kvol])*X1[i+j*kvolHalo]);
#endif
   }
   if(c_sw)
      Clover_free(clover);
#ifdef USE_GPU
#ifdef _DEBUG
   cudaFree(smallPhi);
#else
   cudaFreeAsync(smallPhi,NULL);
#endif
#else
   free(smallPhi);
#endif
   //hg was summed over inside of Average_Plaquette.
#if(nproc>1)
   Par_dsum(&hp); Par_dsum(&hf);
#endif
   *s=hg+hf; *h=(*s)+hp;
#ifdef _DEBUG
   if(!rank)
      printf("hg=%.5e; hf=%.5e; hp=%.5e; h=%.5e\n", hg, hf, hp, *h);
#endif
   return 0;
}

References Average_Plaquette(), AVX, Clover(), Clover_free(), Complex, Complex_f, Congradq(), conj(), creal, cublas_handle, Fill_Small_Phi(), kferm2, kmom, kvol, kvolHalo, nc, ndirac, nf, Par_dsum(), rank, and streams.

Here is the call graph for this function:

Here is the caller graph for this function: