MPI functions

Functions
int	Par_begin (int argc, char *argv[])
	Initialises the MPI configuration.
int	Par_sread (const int iread, const float beta, const float fmu, const float akappa, const Complex_f ajq, const float c_sw, Complex *u11, Complex *u12, Complex *u11t, Complex *u12t)
	Reads and assigns the gauges from file.
int	Par_swrite (const int itraj, const int icheck, const float beta, const float fmu, const float akappa, const Complex_f ajq, const float c_sw, Complex *u11, Complex *u12)
	Copies u11 and u12 into arrays without halos which then get written to output.
int	Par_isum (int *ival)
	Performs a reduction on an integer ival to get a sum which is then distributed to all ranks.
int	Par_dsum (double *dval)
	Performs a reduction on a double dval to get a sum which is then distributed to all ranks.
int	Par_fsum (float *dval)
	Performs a reduction on a float dval to get a sum which is then distributed to all ranks.
int	Par_csum (Complex_f *cval)
	Performs a reduction on a complex float cval to get a sum which is then distributed to all ranks.
int	Par_zsum (Complex *zval)
	Performs a reduction on a complex double zval to get a sum which is then distributed to all ranks.
int	Par_icopy (int *ival)
	Broadcasts an integer to the other processes.
int	Par_dcopy (double *dval)
	Broadcasts a double to the other processes.
int	Par_fcopy (float *fval)
	Broadcasts a float to the other processes.
int	Par_ccopy (Complex *cval)
	Broadcasts a complex float to the other processes.
int	Par_zcopy (Complex *zval)
	Broadcasts a complex double to the other processes.
int	ZHalo_swap_all (Complex *z, int ncpt)
	Calls the functions to send data to both the up and down halos.
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.
int	CHalo_swap_all (Complex_f *c, int ncpt)
	Calls the functions to send data to both the up and down halos.
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.
int	DHalo_swap_all (double *d, int ncpt)
	Calls the functions to send data to both the up and down halos.
int	DHalo_swap_dir (double *d, int ncpt, int idir, int layer)
	Swaps the halos along the axis given by idir in the direction given by layer.
int	SHalo_swap_all (float *d, int ncpt)
	Calls the functions to send data to both the up and down halos.
int	SHalo_swap_dir (float *d, int ncpt, int idir, int layer)
	Swaps the halos along the axis given by idir in the direction given by layer.
int	Trial_Exchange (Complex *ut[2], Complex_f *ut_f[2])
	Exchanges the trial fields.

Detailed Description

Function Documentation

CHalo_swap_all()

int CHalo_swap_all	(	Complex_f *	c,
		int	ncpt )

Calls the functions to send data to both the up and down halos.

Parameters

[in,out]	c	The data being sent
[in]	ncpt	Number of components being sent

Returns

Zero on success, integer error code otherwise

Postcondition

Halo terms of c updated

References Complex_f.

Here is the caller graph for this function:

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.

Parameters

[in,out]	c	The data being moved about. It should be an array of dimension [kvol+halo][something else]
[in]	ncpt	Number of components being sent
[in]	idir	The axis being moved along in C Indexing
[in]	layer	Either DOWN (0) or UP (1)

Returns

Zero on success, Integer Error code otherwise

Postcondition

Halo terms of c updated in direction idir and layer layer

References Complex_f.

Here is the caller graph for this function:

DHalo_swap_all()

int DHalo_swap_all	(	double *	d,
		int	ncpt )

Calls the functions to send data to both the up and down halos.

Parameters

[in,out]	d	The data being sent
[in]	ncpt	Number of components being sent

Returns

Zero on success, integer error code otherwise

Postcondition

Halo terms of d updated

DHalo_swap_dir()

int DHalo_swap_dir	(	double *	d,
		int	ncpt,
		int	idir,
		int	layer )

Swaps the halos along the axis given by idir in the direction given by layer.

Parameters

[in,out]	d	The data being moved about. It should be an array of dimension [kvol+halo][something else]
[in]	ncpt	Number of components being sent
[in]	idir	The axis being moved along in C Indexing
[in]	layer	Either DOWN (0) or UP (1)

Returns

Zero on success, Integer Error code otherwise

Postcondition

Halo terms of d updated in direction idir and layer layer

Here is the caller graph for this function:

Par_begin()

int Par_begin	(	int	argc,
		char *	argv[] )

Initialises the MPI configuration.

Parameters

[in]	argc	Number of arguments given to the programme
[in]	argv	Array of arguments

Returns

Zero on success, integer error code otherwise.

Postcondition

MPI communicator configured, processor and sublattice topology determined

Definition at line 24 of file par_mpi.c.

                                     {
   //TODO: Remove as much non-MPI stuff from here as possible
   const char funcname[] = "Par_begin";
   int size;
#if(nproc>1)
   if(MPI_Init(&argc, &argv)){
      fprintf(stderr, "Error %i in %s: Failed to initialise MPI\nExiting\n\n", NO_MPI_INIT, funcname);
      MPI_Abort(comm,NO_MPI_INIT);
      exit(NO_MPI_INIT);
   }
 
   if(MPI_Comm_rank(comm, &rank)){
      fprintf(stderr, "Error %i in %s: Failed to find rank.\nExiting...\n\n", NO_MPI_RANK, funcname);
      MPI_Abort(comm,NO_MPI_RANK);
   }
   if(MPI_Comm_size(comm, &size)){
      fprintf(stderr, "Error %i in %s: Failed to find size\nExiting...\n\n", NO_MPI_SIZE, funcname);
      MPI_Abort(comm,NO_MPI_SIZE);
   }
#else
   size=1; rank=0;
#endif
   //If size isn't the same as the max allowed number of processes, then there's a problem somewhere.
   if(size!=nproc){
      fprintf(stderr, "Error %i in %s: For process %i, size %i is not equal to nproc %i.\n"
            "Exiting...\n\n", SIZEPROC, funcname, rank, size, nproc);
#if(nproc>1)
      MPI_Abort(comm,SIZEPROC);
#else
      exit(SIZEPROC);
#endif
   }
   //gsize is the size of the system, lsize is the size of each MPI Grid
   int gsize[4], lsize[4];
   gsize[0]=nx; gsize[1]=ny; gsize[2]=nz; gsize[3]=nt;
   lsize[0]=ksizex; lsize[1]=ksizey; lsize[2]=ksizez; lsize[3]=ksizet;
 
   //Topology layout
   int cartsize[ndim] __attribute__((aligned(AVX)));
   cartsize[0]=npx; cartsize[1]=npy; cartsize[2]=npz; cartsize[3]=npt;
 
   //For the topology, says if each dimension is periodic or not
   //Probably for us everything will be but using the four vector
   //gives the choice at least
   int periods[ndim] __attribute__((aligned(AVX)));
#pragma unroll
   for(int i=0; i<ndim; i++)
      periods[i] = true;
   //Not going to change the rank order
   int reorder = false;
   //Declare the topology
#if(nproc>1)
   MPI_Comm commcart;
   MPI_Cart_create(comm, ndim, cartsize, periods, reorder, &commcart);
#endif
 
   //Get nearest neighbours of processors
#if(nproc>1)
#pragma unroll
   for(int i= 0; i<ndim; i++)
      MPI_Cart_shift(commcart, i, 1, &pd[i], &pu[i]);
#endif
   //Get coordinates of processors in the grid
   pcoord = (int*)malloc(ndim*nproc*sizeof(int));
   memset(pcoord,0,sizeof(int)*ndim*nproc);
#if(nproc>1)
   for(int iproc = 0; iproc<nproc; iproc++){
      MPI_Cart_coords(commcart, iproc, ndim, pcoord+iproc*ndim);
#pragma omp simd aligned(pcoord:AVX)
      for(int idim = 0; idim<ndim; idim++){
         pstart[idim][iproc] = pcoord[idim+ndim*iproc]*lsize[idim];
         pstop[idim][iproc]  = pstart[idim][iproc] + lsize[idim];
      }
   }
#else
   //Set iproc=0 because we only have one proc
   for(int idim = 0; idim<ndim; idim++){
      pstart[idim][0] = 0;
      pstop[idim][0]  = lsize[idim];
   }
#endif
#ifdef _DEBUG
   if(!rank)
      printf("Running on %i processors.\nGrid layout is %ix%ix%ix%i\n",
            nproc, npx,npy,npz,npt);
   printf("Rank: %i pu: %i %i %i %i pd: %i %i %i %i\n", rank, pu[0], pu[1], pu[2], pu[3],
         pd[0], pd[1], pd[2], pd[3]);
#endif
   return 0;
}  

References AVX, ksizet, ksizex, ksizey, ksizez, ndim, NO_MPI_INIT, NO_MPI_RANK, NO_MPI_SIZE, nproc, npt, npx, npy, npz, nt, nx, ny, nz, pcoord, pd, pstart, pstop, pu, rank, size, and SIZEPROC.

Here is the caller graph for this function:

Par_ccopy()

int Par_ccopy

(

Complex *

cval

)

Broadcasts a complex float to the other processes.

Parameters

[in,out]

cval

Complex float being broadcast

Returns

Zero on success, integer error code otherwise

Postcondition

All ranks not broadcasting have their value of cval overwritten

References Complex.

Par_csum()

int Par_csum

(

Complex_f *

cval

)

Performs a reduction on a complex float cval to get a sum which is then distributed to all ranks.

Parameters

[in,out]

cval

The pointer to the element being summed, and the container for said sum.

Postcondition

Reduced sum stored in cval

Returns

Zero on success. Integer error code otherwise.

References Complex_f.

Par_dcopy()

int Par_dcopy

(

double *

dval

)

Broadcasts a double to the other processes.

Parameters

[in,out]

dval

double being broadcast

Returns

Zero on success, integer error code otherwise

Postcondition

All ranks not broadcasting have their value of dval overwritten

Here is the caller graph for this function:

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.

Parameters

[in,out]

dval

The pointer to the element being summed, and the container for said sum.

Returns

Zero on success. Integer error code otherwise.

Postcondition

Reduced sum stored in dval

Here is the caller graph for this function:

Par_fcopy()

int Par_fcopy

(

float *

fval

)

Broadcasts a float to the other processes.

Parameters

[in,out]

fval

float being broadcast

Returns

Zero on success, integer error code otherwise

Postcondition

All ranks not broadcasting have their value of fval overwritten

Here is the caller graph for this function:

Par_fsum()

int Par_fsum

(

float *

dval

)

Performs a reduction on a float dval to get a sum which is then distributed to all ranks.

Parameters

[in,out]

dval

The pointer to the element being summed, and the container for said sum.

Returns

Zero on success. Integer error code otherwise.

Postcondition

Reduced sum stored in dval

Here is the caller graph for this function:

Par_icopy()

int Par_icopy

(

int *

ival

)

Broadcasts an integer to the other processes.

Parameters

[in,out]

ival

Integer being broadcast

Returns

Zero on success, integer error code otherwise

Postcondition

All ranks not broadcasting have their value of ival overwritten

Here is the caller graph for this function:

Par_isum()

int Par_isum

(

int *

ival

)

Performs a reduction on an integer ival to get a sum which is then distributed to all ranks.

Parameters

[in,out]

ival

The pointer to the element being summed, and the container for said sum.

Returns

Zero on success. Integer error code otherwise.

Postcondition

Reduced sum stored in ival

Par_sread()

int Par_sread	(	const int	iread,
		const float	beta,
		const float	fmu,
		const float	akappa,
		const Complex_f	ajq,
		const float	c_sw,
		Complex *	u11,
		Complex *	u12,
		Complex *	u11t,
		Complex *	u12t )

Reads and assigns the gauges from file.

Parameters

[in]	iread	Configuration to read in
[in]	beta	Inverse gauge coupling
[in]	fmu	Chemical potential
[in]	akappa	Hopping parameter
[in]	ajq	Diquark source
[in]	c_sw	Clover coefficient
[out]	u11,u12	Gauge fields
[out]	u11t,u12t	Trial fields

Returns

Zero on success, integer error code otherwise

Postcondition

Contents of gauge fields replaced with read in values

Definition at line 114 of file par_mpi.c.

                                                                                 {
   const char funcname[] = "Par_sread";
#if(nproc>1)
   MPI_Status status;
   //For sending the seeds later
   MPI_Datatype MPI_SEED_TYPE = (sizeof(seed)==sizeof(int)) ? MPI_INT:MPI_LONG;
#endif
   //We shall allow the almighty master thread to open the file
   Complex *u1buff = (Complex *)aligned_alloc(AVX,kvol*sizeof(Complex));
   Complex *u2buff = (Complex *)aligned_alloc(AVX,kvol*sizeof(Complex));
   if(!rank){
      //Containers for input. Only needed by the master rank
      Complex *u11Read = (Complex *)aligned_alloc(AVX,ndim*gvol*sizeof(Complex));
      Complex *u12Read = (Complex *)aligned_alloc(AVX,ndim*gvol*sizeof(Complex));
      char gauge_file[FILELEN]="config.";
      int buffer; char buff2[7];
      //Add script for extracting correct mu, j etc.
      buffer = (int)round(100*beta);
      sprintf(buff2,"b%03d",buffer);
      strcat(gauge_file,buff2);
      //κ
      buffer = (int)round(10000*akappa);
      sprintf(buff2,"k%04d",buffer);
      strcat(gauge_file,buff2);
      //μ
      buffer = (int)round(1000*fmu);
      sprintf(buff2,"mu%04d",buffer);
      strcat(gauge_file,buff2);
      //J
      buffer = (int)round(1000*creal(ajq));
      sprintf(buff2,"j%03d",buffer);
      strcat(gauge_file,buff2);
      //c_sw
      if(c_sw!=0){
         buffer = (int)round(100*c_sw);
         sprintf(buff2,"c%02d",buffer);
         strcat(gauge_file,buff2);
      }
      //nx
      sprintf(buff2,"s%02d",nx);
      strcat(gauge_file,buff2);
      //nt
      sprintf(buff2,"t%02d",nt);
      strcat(gauge_file,buff2);
      //nconfig
      char c[8];
      sprintf(c,".%06d", iread);
      strcat(gauge_file, c);
 
      char *fileop = "rb";
      printf("Opening gauge file on processor: %i\n",rank); 
      FILE *con;
      if(!(con = fopen(gauge_file, fileop))){
         fprintf(stderr, "Error %i in %s: Failed to open %s for %s.\
               \nExiting...\n\n", OPENERROR, funcname, gauge_file, fileop);
#if(nproc>1)
         MPI_Abort(comm,OPENERROR);
#endif
         exit(OPENERROR);
      }
      //TODO: SAFETY CHECKS FOR EACH READ OPERATION
      int old_nproc;
      //What was previously the FORTRAN integer is now used to store the number of processors used to
      //generate the configuration
      fread(&old_nproc, sizeof(int), 1, con);
      if(old_nproc!=nproc)
         fprintf(stderr, "Warning %i in %s: Previous run was done on %i processors, current run uses %i.\n",\
               DIFNPROC,funcname,old_nproc,nproc);
      fread(u11Read, ndim*gvol*sizeof(Complex), 1, con);
      fread(u12Read, ndim*gvol*sizeof(Complex), 1, con);
      //The seed array will be used to gather and sort the seeds from each rank so they can be in a continuation run
      //If less processors are used then only nproc seeds are used (breaking the Markov Chain)
      //If more processors are used then we use the first seed to generate the rest as in Par_ranset
#ifdef __RANLUX__
      unsigned long *seed_array=(unsigned long*)calloc(nproc,sizeof(seed));
#elif defined __USE_MKL__ && !defined USE_RAN2
      int *seed_array=(int *)calloc(nproc,sizeof(seed));
#else
      long *seed_array=(long*)calloc(nproc,sizeof(seed));
#endif
      for(int i=0; i<fmin(old_nproc,nproc);i++)
         fread(seed_array+i, sizeof(seed), 1, con);
      fclose(con);
      //Any remaining processors get their initial value set as is done in Par_ranset
      for(int i=old_nproc; i<nproc; i++)
         seed_array[i] = seed_array[0]*(1.0f+8.0f*(float)i/(float)(size-1));
      if(!rank)
         seed=seed_array[0];
#if(nproc>1)
      for(int iproc = 1; iproc<nproc; iproc++)
         if(MPI_Send(&seed_array[iproc], 1, MPI_SEED_TYPE,iproc, 1, comm)){
            fprintf(stderr, "Error %i in %s: Failed to send seed to process %i.\nExiting...\n\n",
                  CANTSEND, funcname, iproc);
            MPI_Abort(comm,CANTSEND);
         }
#endif
 
      for(int iproc = 0; iproc < nproc; iproc++)
         for(int idim = 0; idim < ndim; idim++){
            int i = 0;
            //Index order is reversed from FORTRAN for performance
            //Going to split up assigning icoord[i] to reduce the
            //number of assignments.
            //We're weaving our way through the memory here, converting
            //between lattice and memory coordinates
            for(int it=pstart[3][iproc]; it<pstop[3][iproc]; it++)
               for(int iz=pstart[2][iproc]; iz<pstop[2][iproc]; iz++)
                  for(int iy=pstart[1][iproc]; iy<pstop[1][iproc]; iy++)
                     for(int ix=pstart[0][iproc]; ix<pstop[0][iproc]; ix++){
                        //j is the relative memory index of icoord
                        int j = Coord2gindex(ix,iy,iz,it);
                        //TODO: Don't change this for SoA. We can use this to flip the memory layout without calling a
                        //transpose
                        u1buff[i]=u11Read[idim*gvol+j];
                        u2buff[i]=u12Read[idim*gvol+j];
                        //C starts counting from zero, not 1 so increment afterwards or start at int i=-1
                        i++;
                     }
            if(i!=kvol){
               fprintf(stderr, "Error %i in %s: Number of elements %i is not equal to\
                     kvol %i.\nExiting...\n\n", NUMELEM, funcname, i, kvol);
#if(nproc>1)
               MPI_Abort(comm,NUMELEM);
#else
               exit(NUMELEM);
#endif
            }
            if(!iproc){
#ifdef USE_GPU
               cudaMemcpy(u11+idim*kvol,u1buff,kvol*sizeof(Complex),cudaMemcpyDefault);
               cudaMemcpy(u12+idim*kvol,u2buff,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
               memcpy(u11+idim*kvol,u1buff,kvol*sizeof(Complex));
               memcpy(u12+idim*kvol,u2buff,kvol*sizeof(Complex));
#endif
            }     
#if(nproc>1)
            else{
               //The master thread did all the hard work, the minions just need to receive their
               //data and go.
               if(MPI_Send(u1buff, kvol, MPI_C_DOUBLE_COMPLEX,iproc, 2*idim, comm)){
                  fprintf(stderr, "Error %i in %s: Failed to send ubuff to process %i.\nExiting...\n\n",
                        CANTSEND, funcname, iproc);
#if(nproc>1)
                  MPI_Abort(comm,CANTSEND);
#else
                  exit(CANTSEND);
#endif
               }
               if(MPI_Send(u2buff, kvol, MPI_C_DOUBLE_COMPLEX,iproc, 2*idim+1, comm)){
                  fprintf(stderr, "Error %i in %s: Failed to send ubuff to process %i.\nExiting...\n\n",
                        CANTSEND, funcname, iproc);
#if(nproc>1)
                  MPI_Abort(comm,CANTSEND);
#else
                  exit(CANTSEND);
#endif
               }
            }
#endif
         }
      free(u11Read); free(u12Read);
      free(seed_array);
   }
#if(nproc>1)
   else{
      if(MPI_Recv(&seed, 1, MPI_SEED_TYPE, masterproc, 1, comm, &status)){
         fprintf(stderr, "Error %i in %s: Falied to receive seed on process %i.\nExiting...\n\n",
               CANTRECV, funcname, rank);
#if(nproc>1)
         MPI_Abort(comm,CANTRECV);
#else
         exit(CANTRECV);
#endif
      }
      for(int idim = 0; idim<ndim; idim++){
         //Receiving the data from the master threads.
         if(MPI_Recv(u1buff, kvol, MPI_C_DOUBLE_COMPLEX, masterproc, 2*idim, comm, &status)){
            fprintf(stderr, "Error %i in %s: Falied to receive u11 on process %i.\nExiting...\n\n",
                  CANTRECV, funcname, rank);
            MPI_Abort(comm,CANTRECV);
         }
         if(MPI_Recv(u2buff, kvol, MPI_C_DOUBLE_COMPLEX, masterproc, 2*idim+1, comm, &status)){
            fprintf(stderr, "Error %i in %s: Falied to receive u12 on process %i.\nExiting...\n\n",
                  CANTRECV, funcname, rank);
            MPI_Abort(comm,CANTRECV);
         }
#ifdef USE_GPU
         cudaMemcpy(u11+idim*kvol,u1buff,kvol*sizeof(Complex),cudaMemcpyDefault);
         cudaMemcpy(u12+idim*kvol,u2buff,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
         memcpy(u11+idim*kvol,u1buff,kvol*sizeof(Complex));
         memcpy(u12+idim*kvol,u2buff,kvol*sizeof(Complex));
#endif
      }
   }
#endif
   free(u1buff); free(u2buff);
   for(unsigned short mu=0;mu<ndim;mu++){
#ifdef USE_GPU
      cudaMemcpy(u11t+kvolHalo*mu, u11+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault);
      cudaMemcpy(u12t+kvolHalo*mu, u12+kvol*mu, kvol*sizeof(Complex),cudaMemcpyDefault);
#else
      memcpy(u11t+kvolHalo*mu, u11+kvol*mu, kvol*sizeof(Complex));
      memcpy(u12t+kvolHalo*mu, u12+kvol*mu, kvol*sizeof(Complex));
#endif
   }
   return 0;
}

References AVX, CANTRECV, CANTSEND, Complex, Complex_f, Coord2gindex(), creal, DIFNPROC, FILELEN, gvol, kvol, kvolHalo, masterproc, ndim, nproc, nt, NUMELEM, nx, OPENERROR, pstart, pstop, rank, seed, and size.

Here is the call graph for this function:

Here is the caller graph for this function:

Par_swrite()

int Par_swrite	(	const int	itraj,
		const int	icheck,
		const float	beta,
		const float	fmu,
		const float	akappa,
		const Complex_f	ajq,
		const float	c_sw,
		Complex *	u11,
		Complex *	u12 )

Copies u11 and u12 into arrays without halos which then get written to output.

Modified from an original version of swrite in FORTRAN

Parameters

[in]	itraj	Trajectory to write
[in]	icheck	Not currently used but haven't gotten around to removing it
[in]	beta	Inverse gauge coupling
[in]	fmu	Chemical potential
[in]	akappa	Hopping parameter
[in]	ajq	Diquark source
[in]	c_sw	Clover coefficient
[in]	u11,u12	Gauge fields

Returns

Zero on success, integer error code otherwise

Postcondition

Gauge fields saved to file

Definition at line 324 of file par_mpi.c.

                                                                        {
   const char funcname[] = "par_swrite";
#if (nproc>1)
   MPI_Status status;
   //Used for seed array later on
   MPI_Datatype MPI_SEED_TYPE = (sizeof(seed)==sizeof(int)) ? MPI_INT:MPI_LONG;
#endif
   Complex *u1buff = (Complex *)aligned_alloc(AVX,kvol*sizeof(Complex));
   Complex *u2buff = (Complex *)aligned_alloc(AVX,kvol*sizeof(Complex));
#ifdef _DEBUG
   char dump_prefix[FILELEN]="u11.";
   char dump_buff[32];
   sprintf(dump_buff,"r%01d_c%06d",rank,itraj);
   strcat(dump_prefix,dump_buff);
   FILE *gauge_dump=fopen(dump_prefix,"wb");
   //Print the local trial field in the order it is stored in memory.
   //This is not the same order as it is stored in secondary storage
   fwrite(u11,ndim*kvol*sizeof(Complex),1,gauge_dump);
   fclose(gauge_dump);
#endif
#ifdef __RANLUX__
   seed=gsl_rng_get(ranlux_instd);
#endif
   if(!rank){
      //Array to store the seeds. nth index is the nth processor
#ifdef __RANLUX__
      unsigned long *seed_array=(unsigned long*)calloc(nproc,sizeof(seed));
#elif defined __USE_MKL__ && !defined USE_RAN2
      int *seed_array=(int *)calloc(nproc,sizeof(seed));
#else
      long *seed_array=(long*)calloc(nproc,sizeof(seed));
#endif
      seed_array[0]=seed;
#if(nproc>1)
      for(int iproc = 1; iproc<nproc; iproc++)
         if(MPI_Recv(&seed_array[iproc], 1, MPI_SEED_TYPE,iproc, 1, comm, &status)){
            fprintf(stderr, "Error %i in %s: Failed to receive seed from process %i.\nExiting...\n\n",
                  CANTRECV, funcname, iproc);
            MPI_Abort(comm,CANTRECV);
         }
#endif
      Complex *u11Write = (Complex *)aligned_alloc(AVX,ndim*gvol*sizeof(Complex));
      Complex *u12Write = (Complex *)aligned_alloc(AVX,ndim*gvol*sizeof(Complex));
      //Get correct parts of u11read etc from remote processors
      for(int iproc=0;iproc<nproc;iproc++)
         for(int idim=0;idim<ndim;idim++){
#if(nproc>1)
            if(iproc){
               if(MPI_Recv(u1buff, kvol, MPI_C_DOUBLE_COMPLEX, iproc, 2*idim, comm, &status)){
                  fprintf(stderr, "Error %i in %s: Falied to receive u11 from process %i.\nExiting...\n\n",
                        CANTRECV, funcname, iproc);
                  MPI_Abort(comm,CANTRECV);
               }
               if(MPI_Recv(u2buff, kvol, MPI_C_DOUBLE_COMPLEX, iproc, 2*idim+1, comm, &status)){
                  fprintf(stderr, "Error %i in %s: Falied to receive u12 from process %i.\nExiting...\n\n",
                        CANTRECV, funcname, iproc);
                  MPI_Abort(comm,CANTRECV);
               }
            }
            else{
#endif
               //No need to do MPI Send/Receive on the master rank
               //Array looping is slow so we use memcpy instead
#ifdef USE_GPU
               cudaMemcpy(u1buff,u11+idim*kvol,kvol*sizeof(Complex),cudaMemcpyDefault);
               cudaMemcpy(u2buff,u12+idim*kvol,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
               memcpy(u1buff,u11+idim*kvol,kvol*sizeof(Complex));
               memcpy(u2buff,u12+idim*kvol,kvol*sizeof(Complex));
#endif
#ifdef _DEBUG
               char part_dump[FILELEN]="";
               strcat(part_dump,dump_prefix);
               sprintf(dump_buff,"_d%d",idim);
               strcat(part_dump,dump_buff);
               FILE *pdump=fopen(part_dump,"wb");
               fwrite(u1buff,kvol*sizeof(Complex),1,pdump);
               fclose(pdump);
#endif
#if(nproc>1)
            }
#endif
            int i=0;
            for(int it=pstart[3][iproc]; it<pstop[3][iproc]; it++)
               for(int iz=pstart[2][iproc]; iz<pstop[2][iproc]; iz++)
                  for(int iy=pstart[1][iproc]; iy<pstop[1][iproc]; iy++)
                     for(int ix=pstart[0][iproc]; ix<pstop[0][iproc]; ix++){
                        //j is the relative memory index of icoord
                        int j = Coord2gindex(ix, iy, iz, it);
                        u11Write[idim*gvol+j] = u1buff[i];  
                        u12Write[idim*gvol+j] = u2buff[i];  
                        //C starts counting from zero, not 1 so increment afterwards or start at int i=-1
                        i++;
                     }
            if(i!=kvol){
               fprintf(stderr, "Error %i in %s: Number of elements %i is not equal to\
                     kvol %i.\nExiting...\n\n", NUMELEM, funcname, i, kvol);
#if(nproc>1)
               MPI_Abort(comm,NUMELEM);
#else
               exit(NUMELEM);
#endif
            }
         }
      free(u1buff); free(u2buff);
 
      char gauge_title[FILELEN]="config.";
      int buffer; char buff2[7];
      //Add script for extracting correct mu, j etc.
      buffer = (int)round(100*beta);
      sprintf(buff2,"b%03d",buffer);
      strcat(gauge_title,buff2);
      //κ
      buffer = (int)round(10000*akappa);
      sprintf(buff2,"k%04d",buffer);
      strcat(gauge_title,buff2);
      //μ
      buffer = (int)round(1000*fmu);
      sprintf(buff2,"mu%04d",buffer);
      strcat(gauge_title,buff2);
      //J
      buffer = (int)round(1000*creal(ajq));
      sprintf(buff2,"j%03d",buffer);
      strcat(gauge_title,buff2);
      //c_sw
      if(c_sw!=0){
         buffer = (int)round(100*c_sw);
         sprintf(buff2,"c%02d",buffer);
         strcat(gauge_title,buff2);
      }
      //nx
      sprintf(buff2,"s%02d",nx);
      strcat(gauge_title,buff2);
      //nt
      sprintf(buff2,"t%02d",nt);
      strcat(gauge_title,buff2);
 
      char gauge_file[FILELEN];
      strcpy(gauge_file,gauge_title);
      char c[8];
      sprintf(c,".%06d", itraj);
      strcat(gauge_file, c);
      printf("Gauge file name is %s\n", gauge_file);
      printf("Writing the gauge file on processor %i.\n", rank);
      FILE *con;
      char *fileop = "wb";
      if(!(con=fopen(gauge_file, fileop))){
         fprintf(stderr, "Error %i in %s: Failed to open %s for %s.\
               \nExiting...\n\n", OPENERROR, funcname, gauge_file, fileop);
#if(nproc>1)
         MPI_Abort(comm,OPENERROR);
#else
         exit(OPENERROR);
#endif
      }
      //TODO: SAFETY CHECKS FOR EACH WRITE OPERATION
      //Write the number of processors used in the previous run. This takes the place of the FORTRAN integer rather nicely
#if(nproc==1)
      int size=nproc;
#endif
      fwrite(&size,sizeof(int),1,con);
      fwrite(u11Write, ndim*gvol*sizeof(Complex), 1, con);
      fwrite(u12Write, ndim*gvol*sizeof(Complex), 1, con);
      //TODO
      //Make a seed array, where the nth component is the seed on the nth rank for continuation runs.
      fwrite(seed_array, nproc*sizeof(seed), 1, con);
      fclose(con);
      free(u11Write); free(u12Write);
      free(seed_array);
   }
#if(nproc>1)
   else{
      if(MPI_Send(&seed, 1, MPI_SEED_TYPE, masterproc, 1, comm)){
         fprintf(stderr, "Error %i in %s: Falied to send u11 from process %i.\nExiting...\n\n",
               CANTSEND, funcname, rank);
         MPI_Abort(comm,CANTSEND);
      }
      for(int idim = 0; idim<ndim; idim++){
#ifdef USE_GPU
         cudaMemcpy(u1buff,u11+idim*kvol,kvol*sizeof(Complex),cudaMemcpyDefault);
         cudaMemcpy(u2buff,u12+idim*kvol,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
         memcpy(u1buff,u11+idim*kvol,kvol*sizeof(Complex));
         memcpy(u2buff,u12+idim*kvol,kvol*sizeof(Complex));
#endif
#ifdef _DEBUG
         char part_dump[FILELEN]="";
         strcat(part_dump,dump_prefix);
         sprintf(dump_buff,"_d%d",idim);
         strcat(part_dump,dump_buff);
         FILE *pdump=fopen(part_dump,"wb");
         fwrite(u1buff,kvol*sizeof(Complex),1,pdump);
         fclose(pdump);
#endif
         int i=0;
         if(MPI_Send(u1buff, kvol, MPI_C_DOUBLE_COMPLEX, masterproc, 2*idim, comm)){
            fprintf(stderr, "Error %i in %s: Falied to send u11 from process %i.\nExiting...\n\n",
                  CANTSEND, funcname, rank);
            MPI_Abort(comm,CANTSEND);
         }
         if(MPI_Send(u2buff, kvol, MPI_C_DOUBLE_COMPLEX, masterproc, 2*idim+1, comm)){
            fprintf(stderr, "Error %i in %s: Falied to send u12 from process %i.\nExiting...\n\n",
                  CANTSEND, funcname, rank);
            MPI_Abort(comm,CANTSEND);
         }
      }
      free(u1buff); free(u2buff);
   }
#endif
   return 0;
}

References AVX, CANTRECV, CANTSEND, Complex, Complex_f, Coord2gindex(), creal, FILELEN, gvol, kvol, masterproc, ndim, nproc, nt, NUMELEM, nx, OPENERROR, pstart, pstop, rank, ranlux_instd, seed, and size.

Here is the call graph for this function:

Here is the caller graph for this function:

Par_zcopy()

int Par_zcopy

(

Complex *

zval

)

Broadcasts a complex double to the other processes.

Parameters

[in,out]

zval

Complex double being broadcast

Returns

Zero on success, integer error code otherwise

Postcondition

All ranks not broadcasting have their value of zval overwritten

References Complex.

Par_zsum()

int Par_zsum

(

Complex *

zval

)

Performs a reduction on a complex double zval to get a sum which is then distributed to all ranks.

Parameters

[in,out]

zval

The pointer to the element being summed, and the container for said sum.

Returns

Zero on success. Integer error code otherwise.

Postcondition

Reduced sum stored in zval

References Complex.

SHalo_swap_all()

int SHalo_swap_all	(	float *	d,
		int	ncpt )

Calls the functions to send data to both the up and down halos.

Parameters

[in,out]	d	The data being sent
[in]	ncpt	Number of components being sent

Returns

Zero on success, integer error code otherwise

Postcondition

Halo terms of d updated

Here is the caller graph for this function:

SHalo_swap_dir()

int SHalo_swap_dir	(	float *	d,
		int	ncpt,
		int	idir,
		int	layer )

Swaps the halos along the axis given by idir in the direction given by layer.

Parameters

[in,out]	d	The data being moved about. It should be an array of dimension [kvol+halo][something else]
[in]	ncpt	Number of components being sent
[in]	idir	The axis being moved along in C Indexing
[in]	layer	Either DOWN (0) or UP (1)

Returns

Zero on success, Integer Error code otherwise

Postcondition

Halo terms of d updated in direction idir and layer layer

References Complex, and Complex_f.

Trial_Exchange()

int Trial_Exchange	(	Complex *	ut[2],
		Complex_f *	ut_f[2] )

Exchanges the trial fields.

I noticed that this halo exchange was happening even though the trial fields hadn't been updated. To get around this I'm making a function that does the halo exchange and only calling it after the trial fields get updated.

Parameters

[in,out]	ut	Double precision trial fields
[out]	ut_f	Single precision trial fields

Returns

Zero on success, Integer Error code otherwise

Postcondition

Halos of ut updated. ut_f overwritten with single precision values of ut

Definition at line 1043 of file par_mpi.c.

                                                     {
   /*
    * Exchanges the trial fields. I noticed that this halo exchange was happening
    * even though the trial fields hadn't been updated. To get around this
    * I'm making a function that does the halo exchange and only calling it after
    * the trial fields get updated.
    */
   const char *funchame = "Trial_Exchange";
   //Prefetch the trial fields from the GPU, halos come later
#if(nproc>1)
#ifdef USE_GPU
   int device=-1;
   cudaGetDevice(&device);
   Complex *z;
#ifdef _DEBUG
   cudaMallocManaged((void **)&z,kvolHalo*sizeof(Complex),cudaMemAttachGlobal);
#else
   cudaMallocAsync((void **)&z,kvolHalo*sizeof(Complex),streams[0]);
#endif
   cudaDeviceSynchronise();
#else
   Complex *z = (Complex *)aligned_alloc(AVX,kvolHalo*sizeof(Complex));
#endif
   // ZHalo_swap_all(ut[0],ndim); ZHalo_swap_all(ut[1],ndim);
   // 
   for(int mu=0;mu<ndim;mu++){
      //Copy the column from ut[0]
#ifdef USE_GPU
      cudaMemcpy(z,ut[0]+kvolHalo*mu,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
      memcpy(z,ut[0]+kvolHalo*mu,kvol*sizeof(Complex));
#endif
 
      //Halo exchange on that column
      ZHalo_swap_all(z, 1);
      //And the swap back/getting the next halo
#ifdef USE_GPU
      cudaMemcpy(ut[0]+kvolHalo*mu,z,kvolHalo*sizeof(Complex),cudaMemcpyDefault);
      cudaMemcpy(z,ut[1]+kvolHalo*mu,kvol*sizeof(Complex),cudaMemcpyDefault);
#else
      memcpy(ut[0]+kvolHalo*mu,z,kvolHalo*sizeof(Complex));
      memcpy(z,ut[1]+kvolHalo*mu,kvol*sizeof(Complex));
#endif
 
      //Repeat
      ZHalo_swap_all(z, 1);
#ifdef USE_GPU
      cudaMemcpy(ut[1]+kvolHalo*mu,z,kvolHalo*sizeof(Complex),cudaMemcpyDefault);
#else
      memcpy(ut[1]+kvolHalo*mu,z,kvolHalo*sizeof(Complex));
#endif
   }
   //Now we prefetch the halo
   //And get the single precision gauge fields preppeed
   //Since we want the halos converted too set the stride to one
#ifdef USE_GPU
#ifdef _DEBUG
   cudaFree(z);
#else
   cudaFreeAsync(z,streams[0]);
#endif
   cudaDeviceSynchronise();
#else
   free(z);
#endif
#endif
   ComplexConvert(ut_f[0],ut[0],ndim*kvolHalo,true,1);
   ComplexConvert(ut_f[1],ut[1],ndim*kvolHalo,true,1);
   return 0;
}

References AVX, Complex, Complex_f, ComplexConvert(), cudaDeviceSynchronise, kvol, kvolHalo, ndim, streams, and ZHalo_swap_all().

Here is the call graph for this function:

Here is the caller graph for this function:

ZHalo_swap_all()

int ZHalo_swap_all	(	Complex *	z,
		int	ncpt )

Calls the functions to send data to both the up and down halos.

Parameters

[in,out]	z	The data being sent
[in]	ncpt	Number of components being sent

Returns

Zero on success, integer error code otherwise

Postcondition

Halo terms of z updated

References Complex.

Here is the caller graph for this function:

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.

Parameters

[in,out]	z	The data being moved about. It should be an array of dimension [kvol+halo][something else]
[in]	ncpt	Number of components being sent
[in]	idir	The axis being moved along in C Indexing
[in]	layer	Either DOWN (0) or UP (1)

Returns

Zero on success, Integer Error code otherwise

Postcondition

Halo terms of z updated in direction idir and layer layer

References Complex.

Here is the caller graph for this function: