/github/workspace/main.c

Go to the documentation of this file.

#include <assert.h>
#include <clover.h>
#include <matrices.h>
#ifdef   USE_GPU
#include <cuda_runtime.h>
cublasHandle_t cublas_handle;
cublasStatus_t cublas_status;
cudaMemPool_t mempool;
//Fix this later
#endif
int main(int argc, char *argv[]){
   //Instead of hard coding the function name so the error messages are easier to implement
   const char funcname[] = "main";
 
   Par_begin(argc, argv);
   //Add error catching code...
#if(nproc>1)
   MPI_Comm_rank(comm, &rank);
   MPI_Comm_size(comm, &size);
#endif

   float beta = 1.7f;
   float akappa = 0.1780f;
#ifdef USE_GPU
   __managed__ 
#endif
      Complex_f jqq = 0;
   float fmu = 0.0f; int iread = 0; int istart = 1;
   int iprint = 1; //How often are measurements made
   int icheck = 5; //How often are configurations saved
   int ibound = -1;
 
   float dt=0.004;   float c_sw = 0.0; float delb=0; //Not used?
   float ajq = 0.0;  int stepl = 250;  int ntraj = 10;
   //rank is zero means it must be the "master process"
   if(!rank){
      FILE *midout;
      const char *filename = (argc!=2) ?"midout":argv[1];
      char *fileop = "r";
      if( !(midout = fopen(filename, fileop) ) ){
         fprintf(stderr, "Error %i in %s: Failed to open file %s for %s.\nExiting\n\n",\
               OPENERROR, funcname, filename, fileop);
#if(nproc>1)
         MPI_Abort(comm,OPENERROR);
#else
         exit(OPENERROR);
#endif
      }
      //See the README for what each entry means
      fscanf(midout, "%f %f %f %f %f %f %f %d %d %d %d %d", &dt, &beta, &akappa,\
            &ajq, &c_sw, &fmu, &delb, &stepl, &ntraj, &istart, &icheck, &iread);
      fclose(midout);
      assert(stepl>0);  assert(ntraj>0);    assert(istart>=0);  assert(icheck>0);  assert(iread>=0); 
   }
   //Send inputs to other ranks
#if(nproc>1)
      if(c_sw!=0){
         fprintf(stderr,"Error %i in %s: Multiple MPI Ranks are not currently supported for the clover action.\n"\
               "This is due to the corner halo terms needed to compute the clover force not being implemented.\n"\
               "Please recompile for a single MPI rank or run with c_sw=0.\n"\
               "\nExiting...\n\n",NOIMPL,funcname);
         MPI_Abort(comm,NOIMPL);
         exit(NOIMPL);
      }
   Par_fcopy(&dt); Par_fcopy(&beta); Par_fcopy(&akappa); Par_fcopy(&ajq);
   Par_fcopy(&c_sw); Par_fcopy(&fmu); Par_fcopy(&delb); //Not used?
   Par_icopy(&stepl); Par_icopy(&ntraj); Par_icopy(&istart); Par_icopy(&icheck);
   Par_icopy(&iread); 
#endif
   // Thetaq was never used so depreciated. It's position in midout is now c_sw
   // jqq=ajq*cexp(athq*I);
   jqq=ajq;
   //End of input
#ifdef USE_GPU
   //CUBLAS Handle
   cublasCreate(&cublas_handle);
   //Set up grid and blocks
   blockInit(nx, ny, nz, nt, &dimBlock, &dimGrid);
   //CUDA device
   int device=-1;
   cudaGetDevice(&device);
   //For asynchronous memory, when CUDA syncs any unused memory in the pool is released back to the OS
   //unless a threshold is given. We'll base our threshold off of Congradp
   //12*kvol for the clover and 4*16*kvolHalo for the fermion fields 
   //Factor of 1.5 because we need it in single and double precsion should be plenty without being excessive.
   //Not everything has a halo so larger halos give us more headroom too.
   cudaDeviceGetDefaultMemPool(&mempool, device);
   int threshold=8*kfermHalo*sizeof(Complex);
   cudaMemPoolSetAttribute(mempool, cudaMemPoolAttrReleaseThreshold, &threshold);
#endif
#ifdef _DEBUG
   printf("jqq=%f+(%f)I\n",creal(jqq),cimag(jqq));
#endif
#ifdef _DEBUG
   seed = 967580161;
#else
   seed = time(NULL);
#endif
 
   //Gauge, trial and momentum fields 
   //You'll notice that there are two different allocation/free statements
   //One for CUDA and one for everything else depending on what's
   //being used
   /***  Let's take a quick moment to compare this to the analysis code.
    * The analysis code stores the gauge field as a 4 component real valued vector, whereas the produciton code
    * used two complex numbers.
    *
    * Analysis code: u=(Re(u[0]),Im(u[1]),Re(u[1]),Im(u[0]))
    * Production code: u[0]=u[0]+I*u[3]   u[1]=u[2]+I*u[1]
    *
    */
   Complex *u[2], *ut[2];
   Complex_f *ut_f[2];
   double *dk[2], *pp;
   float *dk_f[2];
   //Halo index arrays
   unsigned int *iu, *id;
   //And clover arrays. These only get assigned if @f$c_\text{SW}>0@f$
   Complex *sigval; Complex_f *sigval_f; unsigned short *sigin;
#ifdef USE_GPU
   //Managed here because it's easier to fill them on CPU
   cudaMallocManaged((void**)&iu,ndim*kvol*sizeof(int),cudaMemAttachGlobal);
   cudaMallocManaged((void**)&id,ndim*kvol*sizeof(int),cudaMemAttachGlobal);
 
   cudaMallocManaged((void **)&dk[0],(kvolHalo)*sizeof(double),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&dk[1],(kvolHalo)*sizeof(double),cudaMemAttachGlobal);
#ifdef _DEBUG
   cudaMallocManaged((void **)&dk_f[0],(kvolHalo)*sizeof(float),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&dk_f[1],(kvolHalo)*sizeof(float),cudaMemAttachGlobal);
#else
   cudaMalloc((void **)&dk_f[0],(kvolHalo)*sizeof(float));
   cudaMalloc((void **)&dk_f[1],(kvolHalo)*sizeof(float));
#endif
 
   unsigned short *gamin; Complex *gamval; Complex_f *gamval_f;
   cudaMallocManaged((void **)&gamval_f,5*4*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&gamval,5*4*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&gamin,4*4*sizeof(short),cudaMemAttachGlobal);
 
   cudaMallocManaged((void **)&u[0],ndim*kvol*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&u[1],ndim*kvol*sizeof(Complex),cudaMemAttachGlobal);
   //Needs to be managed as fermionic.c still used them on CPU
   cudaMallocManaged((void **)&ut[0],ndim*(kvolHalo)*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&ut[1],ndim*(kvolHalo)*sizeof(Complex),cudaMemAttachGlobal);
#ifdef _DEBUG
   cudaMallocManaged((void **)&ut_f[0],ndim*(kvolHalo)*sizeof(Complex_f),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&ut_f[1],ndim*(kvolHalo)*sizeof(Complex_f),cudaMemAttachGlobal);
#else
   cudaMalloc((void **)&ut_f[0],ndim*(kvolHalo)*sizeof(Complex_f));
   cudaMalloc((void **)&ut_f[1],ndim*(kvolHalo)*sizeof(Complex_f));
// cudaMalloc((void **)&ut[0],ndim*(kvolHalo)*sizeof(Complex));
// cudaMalloc((void **)&ut[1],ndim*(kvolHalo)*sizeof(Complex));
#endif
#else
   id = (unsigned int*)aligned_alloc(AVX,ndim*kvol*sizeof(int));
   iu = (unsigned int*)aligned_alloc(AVX,ndim*kvol*sizeof(int));
 
   alignas(AVX) unsigned short gamin[16]; alignas(AVX) Complex gamval[20]; alignas(AVX) Complex_f gamval_f[20];
 
   dk[0] = (double *)aligned_alloc(AVX,(kvolHalo)*sizeof(double));
   dk[1] = (double *)aligned_alloc(AVX,(kvolHalo)*sizeof(double));
   dk_f[0] = (float *)aligned_alloc(AVX,(kvolHalo)*sizeof(float));
   dk_f[1] = (float *)aligned_alloc(AVX,(kvolHalo)*sizeof(float));
 
   u[0] = (Complex *)aligned_alloc(AVX,ndim*kvol*sizeof(Complex));
   u[1] = (Complex *)aligned_alloc(AVX,ndim*kvol*sizeof(Complex));
   ut[0] = (Complex *)aligned_alloc(AVX,ndim*(kvolHalo)*sizeof(Complex));
   ut[1] = (Complex *)aligned_alloc(AVX,ndim*(kvolHalo)*sizeof(Complex));
   ut_f[0] = (Complex_f *)aligned_alloc(AVX,ndim*(kvolHalo)*sizeof(Complex_f));
   ut_f[1] = (Complex_f *)aligned_alloc(AVX,ndim*(kvolHalo)*sizeof(Complex_f));
#endif
   Init(istart,ibound,iread,beta,fmu,akappa,ajq,c_sw,u,ut,ut_f,gamval,gamval_f,gamin,dk,dk_f,iu,id);
   if(c_sw)
      Init_clover(&sigval,&sigval_f,&sigin,c_sw);

#ifdef USE_GPU
   //GPU Initialisation stuff
   Init_CUDA(ut[0],ut[1],gamval,gamval_f,gamin,dk[0],dk[1],iu,id);//&dimBlock,&dimGrid);
#endif
   //Send trials to accelerator for reunitarisation
   Reunitarise(ut);
   //Get trials back
#ifdef USE_GPU
#if(nproc>1) //Memcpy routines need to be strided if there is a halo since the lattice is not contiguous in memory
   for(unsigned short mu=0;mu<ndim;mu++){
      cudaMemcpyAsync(u[0]+kvol*mu, ut[0]+kvolHalo*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
      cudaMemcpyAsync(u[1]+kvol*mu, ut[1]+kvolHalo*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
   }
#else
   cudaMemcpyAsync(u[0], ut[0], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[0]);
   cudaMemcpyAsync(u[1], ut[1], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[1]);
#endif
#else
   for(unsigned short mu=0;mu<ndim;mu++){
      memcpy(u[0]+kvol*mu, ut[0]+kvolHalo*mu, kvol*sizeof(Complex));
      memcpy(u[1]+kvol*mu, ut[1]+kvolHalo*mu, kvol*sizeof(Complex));
   }
#endif
#ifdef USE_GPU
   cudaDeviceSynchronise();
#endif
#ifdef DIAGNOSTIC
   double ancg_diag=0;
   Diagnostics(istart, u, ut, ut_f, iu, id, hu, hd, dk,\
         dk_f, gamin, gamval, gamval_f, sigval, sigval_f, sigin, jqq, akappa, beta, c_sw,ancg_diag);
#endif
 
   //Initial Measurements
   //====================
   Trial_Exchange(ut,ut_f);
   double poly = Polyakov(ut_f);
#ifdef _DEBUG
   if(!rank) printf("Initial Polyakov loop evaluated as %e\n", poly);
#endif
   double hg, avplaqs, avplaqt;
   //Halo exchange of the trial fields
   Average_Plaquette(&hg,&avplaqs,&avplaqt,ut_f,iu,beta);
   //Trajectory length
   double traj=stepl*dt;
   //end trajectory probability
   double proby = 2.5/stepl;
   char suffix[FILELEN]="";
   int buffer; char buff2[7];
   //Add script for extracting correct mu, j etc.
   buffer = (int)round(100*beta);
   sprintf(buff2,"b%03d",buffer);
   strcat(suffix,buff2);
   //κ
   buffer = (int)round(10000*akappa);
   sprintf(buff2,"k%04d",buffer);
   strcat(suffix,buff2);
   //μ
   buffer = (int)round(1000*fmu);
   sprintf(buff2,"mu%04d",buffer);
   strcat(suffix,buff2);
   //J
   buffer = (int)round(1000*ajq);
   sprintf(buff2,"j%03d",buffer);
   strcat(suffix,buff2);
   //c_sw
   if(c_sw){
      buffer = (int)round(100*c_sw);
      sprintf(buff2,"c%03d",buffer);
      strcat(suffix,buff2);
   }
   //nx
   sprintf(buff2,"s%02d",nx);
   strcat(suffix,buff2);
   //nt
   sprintf(buff2,"t%02d",nt);
   strcat(suffix,buff2);
   char outname[FILELEN] = "Output."; char *outop="a";
   strcat(outname,suffix);
   FILE *output;
   if(!rank){
      if(!(output=fopen(outname, outop) )){
         fprintf(stderr,"Error %i in %s: Failed to open file %s for %s.\nExiting\n\n",OPENERROR,funcname,outname,outop);
#if(nproc>1)
         MPI_Abort(comm,OPENERROR);
#else
         exit(OPENERROR);
#endif
      }
      printf("hg = %e, <Ps> = %e, <Pt> = %e, <Poly> = %e\n", hg, avplaqs, avplaqt, poly);
      fprintf(output, "ksize = %i ksizet = %i Nf = %i Halo =%i\nTime step dt = %e Trajectory length = %e\n"\
            "No. of Trajectories = %i β = %e\nκ = %e μ = %e\nDiquark source = %e Clover coefficient = %e\n"\
            "Stopping Residuals: Guidance: %e Acceptance: %e, Estimator: %e\nSeed = %ld\n",
            ksize, ksizet, nf, halo, dt, traj, ntraj, beta, akappa, fmu, ajq, c_sw, rescgg, rescga, respbp, seed);
#ifdef _DEBUG
      //Print to terminal during debugging
      printf("ksize = %i ksizet = %i Nf = %i Halo = %i\nTime step dt = %e Trajectory length = %e\n"\
            "No. of Trajectories = %i β = %e\nκ = %e μ = %e\nDiquark source = %e Clover coefficient = %e\n"\
            "Stopping Residuals: Guidance: %e Acceptance: %e, Estimator: %e\nSeed = %ld\n",
            ksize, ksizet, nf, halo, dt, traj, ntraj, beta, akappa, fmu, ajq, c_sw, rescgg, rescga, respbp, seed);
#endif
   }
   //Initialise for averages
   //======================
   double actiona = 0.0; double vel2a = 0.0; double pbpa = 0.0; double endenfa = 0.0; double denfa = 0.0;
   //Expected canged in Hamiltonian
   double e_dH=0; double e_dH_e=0;
   //Expected Metropolis accept probability. Skewed by cases where the hamiltonian decreases.
   double yav = 0.0; double yyav = 0.0; 
 
   int naccp = 0; int ipbp = 0; int itot = 0;
 
   //Start of classical evolution
   //===========================
   double pbp;
   Complex qq;
   double *dSdpi;
   //Field and related declarations
   Complex *Phi, *X0, *X1;
   //Initialise Arrays. Leaving it late for scoping
   //check the sizes in sizes.h
#ifdef USE_GPU
#ifdef _DEBUG
   cudaMallocManaged((void **)&X0, nf*kferm2*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&Phi, nf*kferm*sizeof(Complex),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&dSdpi, kmom*sizeof(double),cudaMemAttachGlobal);
   cudaMallocManaged((void **)&X1, kferm2Halo*sizeof(Complex),cudaMemAttachGlobal);
#else
   cudaMalloc((void **)&X1, kferm2Halo*sizeof(Complex));
   cudaMalloc((void **)&X0, nf*kferm2*sizeof(Complex));
   cudaMalloc((void **)&Phi, nf*kferm*sizeof(Complex));
   cudaMalloc((void **)&dSdpi, kmom*sizeof(double));
#endif
 
   cudaMallocManaged((void **)&pp, kmom*sizeof(double),cudaMemAttachGlobal);
   cudaDeviceSynchronise();
#else
   Phi= aligned_alloc(AVX,nf*kferm*sizeof(Complex)); 
   X0= aligned_alloc(AVX,nf*kferm2*sizeof(Complex)); 
   X1= aligned_alloc(AVX,kferm2Halo*sizeof(Complex)); 
   dSdpi = aligned_alloc(AVX,kmom*sizeof(double));
   //pp is the momentum field
   pp = aligned_alloc(AVX,kmom*sizeof(double));
#endif
#if (defined SA3AT)
   double start_time=0;
   if(!rank){
#if(nproc>1)
      start_time = MPI_Wtime();
#else
      start_time = omp_get_wtime();
#endif
   }
#endif
   double action;
   //Conjugate Gradient iteration counters
   double ancg,ancgh,totancg,totancgh;
   ancg=ancgh=totancg=totancgh=0;
   for(int itraj = iread+1; itraj <= ntraj+iread; itraj++){
      //Reset conjugate gradient averages
      ancg = 0; ancgh = 0;
#ifdef _DEBUG
      if(!rank)
         printf("Starting itraj %i\n", itraj);
#endif
      Complex_f *clover[2];
      if(c_sw)
         Clover(clover,ut_f,iu,id);
      for(int na=0; na<nf; na++){
         //Probably makes sense to declare this outside the loop
         //but I do like scoping/don't want to break anything else just teat
         //
         //How do we optimise this for use in CUDA? Do we use CUDA's PRNG
         //or stick with MKL and synchronise/copy over the array
#ifdef USE_GPU
         Complex_f *R1_f,*R; Complex *R1;
         cudaMallocManaged((void **)&R,kfermHalo*sizeof(Complex_f),cudaMemAttachGlobal);
#ifdef _DEBUG
         cudaMallocManaged((void **)&R1, kferm*sizeof(Complex),cudaMemAttachGlobal);
         cudaMallocManaged((void **)&R1_f,kferm*sizeof(Complex_f),cudaMemAttachGlobal);
         cudaMemset(R1_f,0,kferm*sizeof(Complex_f));
#else
         cudaMallocAsync((void **)&R1, kferm*sizeof(Complex),streams[1]);
         cudaMallocAsync((void **)&R1_f,kferm*sizeof(Complex_f),streams[0]);
         cudaMemsetAsync(R1_f,0,kferm*sizeof(Complex_f),streams[0]);
#endif
#else
         Complex_f *R=aligned_alloc(AVX,kfermHalo*sizeof(Complex_f));
         Complex *R1= aligned_alloc(AVX,kferm*sizeof(Complex));
         Complex_f *R1_f=aligned_alloc(AVX,kferm*sizeof(Complex_f));
         memset(R1_f,0,kferm*sizeof(Complex_f));
#endif
         //The FORTRAN code had two Gaussian routines.
         //gaussp was the normal Box-Muller and gauss0 didn't have 2 inside the square root
         //Using σ=1/sqrt(2) in these routines has the same effect as gauss0
#if (defined USE_GPU && defined _DEBUG)
         //cudaMemPrefetchAsync(R1_f,kferm*sizeof(Complex_f),device,streams[1]);
#endif
         //Split into chunks to take into account the halos.
         for(unsigned short j=0;j<nc*ngorkov;j++)
            Gauss_c(R+j*kvolHalo,kvol , 0, 1/sqrt(2));
 
         Dslashd_f(R1_f,R,ut_f,iu,id,gamval_f,gamin,dk_f,jqq,akappa);
         if(c_sw)
            ByClover_f(R1_f,R,clover,sigval_f,akappa,sigin,true);
#ifdef USE_GPU
         //Make sure the multiplication is finished before freeing its input!!
         cudaFree(R);//cudaDeviceSynchronise(); 
                     //cudaFree is blocking so don't need to synchronise
         cuComplex_convert(R1_f,R1,kferm,false,dimBlock,dimGrid);
#ifdef _DEBUG
         cudaFree(R1_f);
#else
         //Stream needs to wait for conversion to complete
         cudaDeviceSynchronise();
         cudaFreeAsync(R1_f,streams[0]);
#endif
         cudaMemcpyAsync(Phi+na*kferm,R1, kferm*sizeof(Complex),cudaMemcpyDefault,streams[1]);
#else
         free(R); 
#pragma omp simd aligned(R1_f,R1:AVX)
         for(int i=0;i<kferm;i++)
            R1[i]=(Complex)R1_f[i];
         free(R1_f);
         memcpy(Phi+na*kferm,R1, kferm*sizeof(Complex));
         //Up/down partitioning (using only pseudofermions of flavour 1)
#endif
         UpDownPart(na, X0, R1);
#ifdef USE_GPU
#ifdef _DEBUG
         cudaFree(R1);
#else
         //Stream needs to wait for UpDownPart to complete
         cudaDeviceSynchronise();
         cudaFreeAsync(R1,streams[1]);
#endif
#else
         free(R1);
#endif
      }  
      if(c_sw)
         Clover_free(clover);
      //Heatbath
      //========
      //We're going to make the most of the new Gauss_d routine to send a flattened array
      //and do this all in one step.
#ifdef USE_GPU
#if(nproc>1)//Strided memcpy
      for(unsigned short mu=0;mu<ndim;mu++){
         cudaMemcpyAsync(ut[0]+kvolHalo*mu, u[0]+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
         cudaMemcpyAsync(ut[1]+kvolHalo*mu, u[1]+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
      }
#else 
      cudaMemcpyAsync(ut[0], u[0], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[0]);
      cudaMemcpyAsync(ut[1], u[1], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[1]);
#endif
      cudaDeviceSynchronise();
#else
      for(unsigned short mu=0;mu<ndim;mu++){
         memcpy(ut[0]+kvolHalo*mu, u[0]+kvol*mu, kvol*sizeof(Complex));
         memcpy(ut[1]+kvolHalo*mu, u[1]+kvol*mu, kvol*sizeof(Complex));
      }
#endif
      Trial_Exchange(ut,ut_f);
      Gauss_d(pp, kmom, 0, 1);
 
      //Initialise Trial Fields
      //pp is random at this point so swapping the order isn't really necessary. But it does ensure that it matches
      //previous results 
      // Transpose_d(pp,nadj*ndim,kvol);
      double H0, S0;
      Hamilton(&H0,&S0,rescga,pp,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
            jqq,akappa,beta,c_sw,&ancgh,itraj);
#ifdef _DEBUG
      if(!rank) printf("H0: %e S0: %e\n", H0, S0);
#endif
      if(itraj==1)
         action = S0/gvol;
 
      //Integration 
#if (defined INT_LPFR && defined INT_OMF2) ||(defined INT_LPFR && defined INT_OMF4)||(defined INT_OMF2 && defined INT_OMF4)
#error "Only one integrator may be defined"
#elif defined INT_LPFR
      Leapfrog(ut,ut_f,X0,X1,Phi,dk,dk_f,dSdpi,pp,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,\
            jqq,beta,akappa,c_sw,stepl,dt,&ancg,&itot,proby);
#elif defined INT_OMF2
      OMF2(ut,ut_f,X0,X1,Phi,dk,dk_f,dSdpi,pp,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,\
            jqq,beta,akappa,c_sw,stepl,dt,&ancg,&itot,proby);
#elif defined INT_OMF4
#warning "OMF4 can be less efficient than OMF2 in certain cases. Use with caution. See http://dx.doi.org/10.1103/PhysRevE.73.036706"
      OMF4(ut,ut_f,X0,X1,Phi,dk,dk_f,dSdpi,pp,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,\
            jqq,beta,akappa,c_sw,stepl,dt,&ancg,&itot,proby);
#else
#error "No integrator defined. Please define {INT_LPFR.INT_OMF2,INT_OMF4}"
#endif
 
      totancg+=ancg;
      //Monte Carlo step: Accept new fields with the probability of min(1,exp(H0-X0))
      //Kernel Call needed here?
      Reunitarise(ut);
      double H1, S1;
      Hamilton(&H1,&S1,rescga,pp,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,\
            jqq,akappa,beta,c_sw,&ancgh,itraj);
      ancgh/=2.0; //Hamilton is called at start and end of trajectory
      totancgh+=ancgh;
#ifdef _DEBUG
      printf("H0-H1=%f-%f",H0,H1);
#endif
      double dH = H0 - H1;
#ifdef _DEBUG
      printf("=%f\n",dH);
#endif
      double dS = S0 - S1;
      if(!rank){
         fprintf(output, "dH = %e dS = %e\n", dH, dS);
#ifdef _DEBUG
         printf("dH = %e dS = %e\n", dH, dS);
#endif
      }
      e_dH+=dH; e_dH_e+=dH*dH;
      double y = exp(dH);
      yav+=y;
      yyav+=y*y;
      //The Monte-Carlo
      //Always update  dH is positive (gone from higher to lower energy)
      bool acc;
      if(dH>0 || Par_granf()<=y){
         //Step is accepted. Set s=st
         if(!rank)
            printf("New configuration accepted on trajectory %i.\n", itraj);
         //Original FORTRAN Comment:
         //JIS 20100525: write config here to preempt troubles during measurement!
         //JIS 20100525: remove when all is ok....
#ifdef USE_GPU
#if(nproc>1) //strided Memcpy
         for(unsigned short mu=0;mu<ndim;mu++){
            cudaMemcpyAsync(u[0]+kvol*mu,ut[0]+kvolHalo*mu,kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
            cudaMemcpyAsync(u[1]+kvol*mu,ut[1]+kvolHalo*mu,kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
         }
#else
         cudaMemcpyAsync(u[0],ut[0],ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[0]);
         cudaMemcpyAsync(u[1],ut[1],ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[1]);
#endif
#else
         for(unsigned short mu=0;mu<ndim;mu++){
            memcpy(u[0]+kvol*mu,ut[0]+kvolHalo*mu,kvol*sizeof(Complex));
            memcpy(u[1]+kvol*mu,ut[1]+kvolHalo*mu,kvol*sizeof(Complex));
         }
#endif
         naccp++;
         //Divide by gvol since we've summed over all lattice sites
         action=S1/gvol;
         acc=true;
      }
      else{
         if(!rank)
            printf("New configuration rejected on trajectory %i.\n", itraj);
         acc=false;
      }
      actiona+=action; 
      double vel2=0.0;
#ifdef USE_GPU
      cublasDnrm2(cublas_handle,kmom, pp, 1,&vel2);
      vel2*=vel2;
      cudaDeviceSynchronise();
#elif defined USE_BLAS
      vel2 = cblas_dnrm2(kmom, pp, 1);
      vel2*=vel2;
#else
#pragma unroll
      for(int i=0; i<kmom; i++)
         vel2+=pp[i]*pp[i];
#endif
#if(nproc>1)
      Par_dsum(&vel2);
#endif
      vel2a+=vel2/(ndim*nadj*gvol);
 
      if(itraj%iprint==0){
         //If rejected, copy the previously accepted field in for measurements
         if(!acc){
#ifdef USE_GPU
#if(nproc>1) //Strided Memcpy
            for(unsigned short mu=0;mu<ndim;mu++){
               cudaMemcpyAsync(ut[0]+kvolHalo*mu, u[0]+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
               cudaMemcpyAsync(ut[1]+kvolHalo*mu, u[1]+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault,streams[mu]);
            }
#else
            cudaMemcpyAsync(ut[0], u[0], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[0]);
            cudaMemcpyAsync(ut[1], u[1], ndim*kvol*sizeof(Complex),cudaMemcpyDefault,streams[1]);
#endif
            cudaDeviceSynchronise();
#else
            for(unsigned short mu=0;mu<ndim;mu++){
               memcpy(ut[0]+kvolHalo*mu, u[0]+kvol*mu, kvol*sizeof(Complex));
               memcpy(ut[1]+kvolHalo*mu, u[1]+kvol*mu, kvol*sizeof(Complex));
            }
#endif
            Trial_Exchange(ut,ut_f);
         }
#ifdef _DEBUG
         if(!rank)
            printf("Starting measurements\n");
#endif
         int itercg=0;
         double endenf, denf;
         Complex qbqb;
         //Stop gap for measurement failure on Kay;
         //If the Congrad in Measure fails, don't measure the Diquark or PBP-Density observables for
         //that trajectory
         int measure_check=0;
         measure_check = Measure(&pbp,&endenf,&denf,&qq,&qbqb,respbp,&itercg,ut,ut_f,iu,id,\
               gamval,gamval_f,gamin,sigval,sigval_f,sigin,dk,dk_f,jqq,akappa,c_sw,Phi);
#ifdef _DEBUG
         if(!rank)
            printf("Finished measurements\n");
#endif
         pbpa+=pbp; endenfa+=endenf; denfa+=denf; ipbp++;
         Average_Plaquette(&hg,&avplaqs,&avplaqt,ut_f,iu,beta);
         poly = Polyakov(ut_f);
         //We have four output files, so may as well get the other ranks to help out
         //and abuse scoping rules while we're at it.
         //Can use either OpenMP or MPI to do this
#if (nproc>=4)
         switch(rank)
#else
            if(!rank)
#pragma omp parallel for
               for(int i=0; i<4; i++)
                  switch(i)
#endif
                  {
                     case(0): 
                        //Output code... Some files weren't opened in the main loop of the FORTRAN code 
                        //That will need to be looked into for the C version
                        //It would explain the weird names like fort.1X that looked like they were somehow
                        //FORTRAN related...
                        fprintf(output, "Measure (CG) %i Update (CG) %.3f Hamiltonian (CG) %.3f\n", itercg, ancg, ancgh);
                        fflush(output);
                        break;
                     case(1):
                        {
                           FILE *fortout;
                           char fortname[FILELEN] = "fermi.";
                           strcat(fortname,suffix);
                           const char *fortop= (itraj==1) ? "w" : "a";
                           if(!(fortout=fopen(fortname, fortop) )){
                              fprintf(stderr, "Error %i in %s: Failed to open file %s for %s.\nExiting\n\n",\
                                    OPENERROR, funcname, fortname, fortop);
#if(nproc>1)
                              MPI_Abort(comm,OPENERROR);
#else
                              exit(OPENERROR);
#endif
                           }
                           if(itraj==1)
                              fprintf(fortout, "pbp\tendenf\tdenf\n");
                           if(measure_check)
                              fprintf(fortout, "%e\t%e\t%e\n", NAN, NAN, NAN);
                           else
                              fprintf(fortout, "%e\t%e\t%e\n", pbp, endenf, denf);
                           fclose(fortout);
                           break;
                        }
                     case(2):
                        //The original code implicitly created these files with the name
                        //fort.XX where XX is the file label
                        //from FORTRAN. This was fort.12
                        {
                           FILE *fortout;
                           char fortname[FILELEN] = "bose."; 
                           strcat(fortname,suffix);
                           const char *fortop= (itraj==1) ? "w" : "a";
                           if(!(fortout=fopen(fortname, fortop) )){
                              fprintf(stderr, "Error %i in %s: Failed to open file %s for %s.\nExiting\n\n",\
                                    OPENERROR, funcname, fortname, fortop);
                           }
                           if(itraj==1)
                              fprintf(fortout, "avplaqs\tavplaqt\tpoly\n");
                           fprintf(fortout, "%e\t%e\t%e\n", avplaqs, avplaqt, poly);
                           fclose(fortout);
                           break;
                        }
                     case(3):
                        {
                           FILE *fortout;
                           char fortname[FILELEN] = "diq.";
                           strcat(fortname,suffix);
                           const char *fortop= (itraj==1) ? "w" : "a";
                           if(!(fortout=fopen(fortname, fortop) )){
                              fprintf(stderr, "Error %i in %s: Failed to open file %s for %s.\nExiting\n\n",\
                                    OPENERROR, funcname, fortname, fortop);
#if(nproc>1)
                              MPI_Abort(comm,OPENERROR);
#else
                              exit(OPENERROR);
#endif
                           }
                           if(itraj==1)
                              fprintf(fortout, "Re(qq)\n");
                           if(measure_check)
                              fprintf(fortout, "%e\n", NAN);
                           else
                              fprintf(fortout, "%e\n", creal(qq));
                           fclose(fortout);
                           break;
                        }
                     default: break;
                  }
      }
      if(itraj%icheck==0){
         Par_swrite(itraj,icheck,beta,fmu,akappa,ajq,c_sw,u[0],u[1]);
      }
      if(!rank)
         fflush(output);
   }
#if (defined SA3AT)
   double elapsed = 0;
   if(!rank){
#if(nproc>1)
      elapsed = MPI_Wtime()-start_time;
#else
      elapsed = omp_get_wtime()-start_time;
#endif
   }
#endif
   //End of main loop
   //Free arrays
#ifdef USE_GPU
   //Make a routine that does this for us
   cudaFree(dk[0]); cudaFree(dk[1]); cudaFree(dSdpi); cudaFree(pp);
   cudaFree(Phi); cudaFree(ut[0]); cudaFree(ut[1]);
   cudaFree(X0); cudaFree(X1); cudaFree(u[0]); cudaFree(u[1]);
   cudaFree(id); cudaFree(iu); 
   cudaFree(dk_f[0]); cudaFree(dk_f[1]); cudaFree(ut_f[0]); cudaFree(ut_f[1]);
   cudaFree(gamin); cudaFree(gamval); cudaFree(gamval_f);
   if(c_sw){
      cudaFree(sigval); cudaFree(sigval_f); cudaFree(sigin);
   }
   cublasDestroy(cublas_handle);
#else
   free(dk[0]); free(dk[1]);free(dSdpi); free(pp);
   free(Phi); free(ut[0]); free(ut[1]);
   free(X0); free(X1); free(u[0]); free(u[1]);
   free(id); free(iu);
   free(dk_f[0]); free(dk_f[1]); free(ut_f[0]); free(ut_f[1]);
   if(c_sw){
      free(sigval); free(sigval_f); free(sigin);
   }
#endif
   free(hd); free(hu); free(pcoord);
#ifdef __RANLUX__
   gsl_rng_free(ranlux_instd);
#endif
#if (defined SA3AT)
   if(!rank){
      FILE *sa3at = fopen("Bench_times.csv", "a");
#ifdef USE_GPU
      char version[256];
      int cuversion; cudaRuntimeGetVersion(&cuversion);
      sprintf(version,"CUDA %d\tBlock: (%d,%d,%d)\tGrid: (%d,%d,%d)\n%s\n",cuversion,\
            dimBlock.x,dimBlock.y,dimBlock.z,dimGrid.x,dimGrid.y,dimGrid.z,__VERSION__);
#else
      char *version=__VERSION__;
#endif
      fprintf(sa3at, "%s\nβ%0.3f κ:%0.4f μ:%0.4f j:%0.3f s:%i t:%i kvol:%ld\n"
            "npx:%i npt:%i nthread:%i ncore:%i time:%f traj_time:%f\n\n",\
            version,beta,akappa,fmu,ajq,nx,nt,kvol,npx,npt,nthreads,npx*npy*npz*npt*nthreads,elapsed,elapsed/ntraj);
      fclose(sa3at);
   }
#endif
   //Get averages for final output
   actiona/=ntraj; vel2a/=ntraj; pbpa/=ipbp; endenfa/=ipbp; denfa/=ipbp;
   totancg/=ntraj; totancgh/=ntraj; 
   e_dH/=ntraj; e_dH_e=sqrt((e_dH_e/ntraj-e_dH*e_dH)/(ntraj-1));
   yav/=ntraj; yyav=sqrt((yyav/ntraj - yav*yav)/(ntraj-1));
   float traj_cost=totancg/dt;
   double atraj=dt*itot/ntraj;
 
   if(!rank){
      fprintf(output, "Averages for the last %i trajectories\n"\
            "Number of acceptances: %i\tAverage Trajectory Length = %e\n"\
            "<dH>=%e+/-%e\t<exp(dH)>=%e+/-%e\tTrajectory cost=N_cg/dt =%e\n"\
            "Average number of congrad iter guidance: %.3f acceptance %.3f\n"\
            "psibarpsi = %e\n"\
            "Mean Square Velocity = %e\tAction Per Site = %e\n"\
            "Energy Density = %e\tNumber Density %e\n\n\n",\
            ntraj, naccp, atraj, e_dH,e_dH_e, yav, yyav, traj_cost, totancg, totancgh, pbpa, vel2a, actiona, endenfa, denfa);
      fclose(output);
   }
#if(nproc>1)
   //Ensure writing is done before finalising just in case finalise segfaults and crashes the other ranks mid-write
   MPI_Barrier(comm);
   MPI_Finalise();
#endif
   fflush(stdout);
   return 0;
}