Commit f61fbf14 authored by Nandhana Sakhtivel's avatar Nandhana Sakhtivel
Browse files

Revamp of Read Data

parent 677aac52

allvars.c

+7 −4
Original line number Diff line number Diff line
@@ -8,17 +8,20 @@ struct op out = {"grid.txt", "coords.txt", "grid_real.bin", "grid_img.bin", "fft 


struct meta metaData;

struct time timing;

struct parameter param;

struct fileData data;



char filename[1000];

int num_threads;

char datapath_multi[NFILES][900], datapath[900];

char datapath[900];

int xaxis, yaxis;

int ndatasets;

int grid_size_x = 2048;

int grid_size_y = 2048;

int num_w_planes = 8;

int rank;

int size;

long nsectors;

long startrow;

double resolution;



clock_t start, end, start0, startk, endk;

struct timespec begin, finish, begin0, begink, finishk;
 
long nsectors;

allvars.h

+20 −3
Original line number Diff line number Diff line
@@ -91,17 +91,34 @@ extern struct time 


} timing;



extern struct parameter

{

        int num_threads;

        int ndatasets;

        char datapath_multi[NFILES][900];

} param;



extern struct fileData

{

        double * uu;

        double * vv;

        double * ww;

        float * weights;

        float * visreal;

        float * visimg;

}data;



extern char filename[1000];

extern int num_threads;

extern char datapath_multi[NFILES][900],datapath[900];

extern char datapath[900];

extern int xaxis, yaxis;

extern int ndatasets;

extern int grid_size_x;

extern int grid_size_y;

extern int num_w_planes;

extern int rank;

extern int size;

extern long nsectors;

extern long startrow;

extern double resolution;



extern clock_t start, end, start0, startk, endk;

extern struct timespec begin, finish, begin0, begink, finishk;

init.c

+202 −101
Original line number Diff line number Diff line
@@ -4,7 +4,7 @@ 
#include "allvars.h"

#include "proto.h"



void init(int argc, char * argv[])

void init()

{

   // MESH SIZE

   int local_grid_size_x;// = 8;
@@ -12,41 +12,11 @@ void init(int argc, char * argv[]) 
   

   read_parameter_file(in.paramfile);

   

   if ( num_threads == 0 )

   if ( param.num_threads == 0 )

   {

    fprintf(stderr, "Usage: %s number_of_OMP_Threads \n", num_threads);

    fprintf(stderr, "Usage: %s number_of_OMP_Threads \n", param.num_threads);

    exit(1);

   }

   clock_gettime(CLOCK_MONOTONIC, &begin0);

   start0 = clock();

   // Intialize MPI environment

   #ifdef USE_MPI

	MPI_Init(&argc,&argv);

	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	MPI_Comm_size(MPI_COMM_WORLD, &size);

	if(rank == 0)printf("Running with %d MPI tasks\n",size);

  	#ifdef USE_FFTW

		fftw_mpi_init();

  	#endif

   #else

	rank = 0;

	size = 1;

   #endif



   if(rank == 0)printf("Running with %d threads\n",num_threads);



   #ifdef ACCOMP

	if(rank == 0){

  		if (0 == omp_get_num_devices()) {

      			printf("No accelerator found ... exit\n");

      			exit(255);

   		}

   		printf("Number of available GPUs %d\n", omp_get_num_devices());

  		#ifdef NVIDIA

      			prtAccelInfo();

   		#endif

 	}

   #endif



   // set the local size of the image

   local_grid_size_x = grid_size_x;
@@ -63,15 +33,67 @@ void init(int argc, char * argv[]) 
   start = clock();



   // INPUT FILES (only the first ndatasets entries are used)

   strcpy(datapath,datapath_multi[0]);

   strcpy(datapath,param.datapath_multi[0]);

                   

   // Read metadata

   fileName(datapath, in.metafile);

   readMetaData(filename);



   // Local Calculation

   int nsub = 1000;

   printf("Subtracting last %d measurements\n",nsub);

   metaData.Nmeasures = metaData.Nmeasures-nsub;

   metaData.Nvis = metaData.Nmeasures*metaData.freq_per_chan*metaData.polarisations;



   // calculate the coordinates of the center

   double uvshift = metaData.uvmin/(metaData.uvmax-metaData.uvmin);

   

   if (rank == 0)

   {

        printf("N. measurements %ld\n",metaData.Nmeasures);

        printf("N. visibilities %ld\n",metaData.Nvis);

   }



  // Set temporary local size of points

  long nm_pe = (long)(metaData.Nmeasures/size);

  long remaining = metaData.Nmeasures%size;



  startrow = rank*nm_pe;

  if (rank == size-1)nm_pe = nm_pe+remaining;



  long Nmeasures_tot = metaData.Nmeasures;

  metaData.Nmeasures = nm_pe;

  long Nvis_tot = metaData.Nvis;

  metaData.Nvis = metaData.Nmeasures*metaData.freq_per_chan*metaData.polarisations;

  metaData.Nweights = metaData.Nmeasures*metaData.polarisations;



  #ifdef VERBOSE

        printf("N. measurements on %d %ld\n",rank,metaData.Nmeasures);

        printf("N. visibilities on %d %ld\n",rank,metaData.Nvis);

  #endif



  // Allocate Data Buffer

  allocate_memory();

  

  // Reading Data

  readData();

 

  #ifdef USE_MPI

        MPI_Barrier(MPI_COMM_WORLD);

  #endif



  clock_gettime(CLOCK_MONOTONIC, &finish);

  end = clock();

  timing.setup_time = ((double) (end - start)) / CLOCKS_PER_SEC;

  timing.setup_time1 = (finish.tv_sec - begin.tv_sec);

  timing.setup_time1 += (finish.tv_nsec - begin.tv_nsec) / 1000000000.0;

 

}





void read_parameter_file(char *fname)

{

   if(rank == 0)

   {

   	if(file.pFile = fopen (fname,"r"))

   	{
@@ -80,15 +102,15 @@ void read_parameter_file(char *fname) 
      		{

       			if(strcmp(buf1, "num_threads") == 0)

        		{

           num_threads = atoi(buf2);

           			param.num_threads = atoi(buf2);

        		}

        		if(strcmp(buf1, "Datapath") == 0)

       			{

            strcpy(datapath_multi[0], buf2);

            			strcpy(param.datapath_multi[0], buf2);

        		}

        		if(strcmp(buf1, "ndatasets") == 0)

        		{

           ndatasets = atoi(buf2);

           			param.ndatasets = atoi(buf2);

        		}

        		if(strcmp(buf1, "ufile") == 0)

        		{
@@ -129,12 +151,19 @@ void read_parameter_file(char *fname) 
       }

    }

   

    /* Communicating the relevent parameters to the other process */



    MPI_Bcast(&in, sizeof(struct ip), MPI_BYTE, 0, MPI_COMM_WORLD);

    MPI_Bcast(&param, sizeof(struct parameter), MPI_BYTE, 0, MPI_COMM_WORLD);

}

void fileName(char datapath[900], char file[30]) {

     strcpy(filename,datapath);

     strcat(filename,file);

}



void readMetaData(char fileLocal[1000]) {



     if(rank == 0) {

     	file.pFile = fopen (fileLocal,"r");

     	fscanf(file.pFile,"%ld",&metaData.Nmeasures);

     	fscanf(file.pFile,"%ld",&metaData.Nvis);
@@ -150,3 +179,75 @@ void readMetaData(char fileLocal[1000]) { 
     	fscanf(file.pFile,"%lf",&metaData.wmax);

     	fclose(file.pFile);

      }

      

      /* Communicating the relevent parameters to the other process */



      MPI_Bcast(&metaData, sizeof(struct meta), MPI_BYTE, 0, MPI_COMM_WORLD);

}



void allocate_memory() {





     // DAV: all these arrays can be allocatate statically for the sake of optimization. However be careful that if MPI is used

     // all the sizes are rescaled by the number of MPI tasks

     //  Allocate arrays



     data.uu = (double*) calloc(metaData.Nmeasures,sizeof(double));

     data.vv = (double*) calloc(metaData.Nmeasures,sizeof(double));

     data.ww = (double*) calloc(metaData.Nmeasures,sizeof(double));

     data.weights = (float*) calloc(metaData.Nweights,sizeof(float));

     data.visreal = (float*) calloc(metaData.Nvis,sizeof(float));

     data.visimg = (float*) calloc(metaData.Nvis,sizeof(float));

}



void readData() {



  //  if(rank == 0) {

         printf("READING DATA\n");



         fileName(datapath, in.ufile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

         fread(data.uu,metaData.Nmeasures*sizeof(double),1,file.pFile);

         fclose(file.pFile);



         fileName(datapath, in.vfile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

         fread(data.vv,metaData.Nmeasures*sizeof(double),1,file.pFile);

         fclose(file.pFile);



         fileName(datapath, in.wfile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

         fread(data.ww,metaData.Nmeasures*sizeof(double),1,file.pFile);

         fclose(file.pFile);



    /*     fileName(datapath, in.weightsfile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*metaData.polarisations*sizeof(float),SEEK_SET);

         fread(data.weights,(metaData.Nweights)*sizeof(float),1,file.pFile);

         fclose(file.pFile);



         fileName(datapath, in.visrealfile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);

         fread(data.visreal,metaData.Nvis*sizeof(float),1,file.pFile);

         fclose(file.pFile);



         fileName(datapath, in.visimgfile);

         file.pFile = fopen (filename,"rb");

         fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);

         fread(data.visimg,metaData.Nvis*sizeof(float),1,file.pFile);

         fclose(file.pFile);

     */

//     }



     #ifdef USE_MPI

        MPI_Barrier(MPI_COMM_WORLD);

     #endif



     /* Communicating the relevent parameters to the other process */



//     MPI_Bcast(&data, sizeof(struct fileData), MPI_BYTE, 0, MPI_COMM_WORLD);

}

main.c

+50 −170
Original line number Diff line number Diff line
@@ -23,19 +23,6 @@ int main(int argc, char * argv[]) 


	char srank[4];



	double * uu;

	double * vv;

	double * ww;

	float * weights;

	float * visreal;

	float * visimg;

	double resolution;



  // MESH SIZE

//	int local_grid_size_x;// = 8;

//	int local_grid_size_y;// = 8;

//	int xaxis;

//	int yaxis;



	// DAV: the corresponding KernelLen is calculated within the wstack function. It can be anyway hardcoded for optimization

	int w_support = 7;
@@ -44,7 +31,6 @@ int main(int argc, char * argv[]) 
	double w_supporth = (double)((w_support-1)/2)*dx;



	double elapsed;

//	long nsectors;

  

        if(argc > 1)

        {
@@ -56,8 +42,11 @@ int main(int argc, char * argv[]) 
          exit(1);

        }

 

        init(argc, argv); 

  /*       

        clock_gettime(CLOCK_MONOTONIC, &begin0);

   	start0 = clock();

 

        /* Initializing MPI Environment */



        #ifdef USE_MPI

		MPI_Init(&argc,&argv);

		MPI_Comm_rank(MPI_COMM_WORLD, &rank);
@@ -70,9 +59,6 @@ int main(int argc, char * argv[]) 
		rank = 0;

		size = 1;

	#endif



  if(rank == 0)printf("Running with %d threads\n",num_threads);



	#ifdef ACCOMP

	if(rank == 0){

  		if (0 == omp_get_num_devices()) {
@@ -86,114 +72,8 @@ if(rank == 0){ 
 	}

	#endif













	// set the local size of the image

	local_grid_size_x = grid_size_x;

	nsectors = NUM_OF_SECTORS;

	if (nsectors < 0) nsectors = size;

	local_grid_size_y = grid_size_y/nsectors;

	//nsectors = size;



	// LOCAL grid size

	xaxis = local_grid_size_x;

	yaxis = local_grid_size_y;



	clock_gettime(CLOCK_MONOTONIC, &begin);

	start = clock();



	// INPUT FILES (only the first ndatasets entries are used)

    	strcpy(datapath,datapath_multi[0]);

	

        // Read metadata

        fileName(datapath, in.metafile);

        readMetaData(filename);

*/

	// WATCH THIS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	int nsub = 1000;

	//int nsub = 10;

	printf("Subtracting last %d measurements\n",nsub);

  	metaData.Nmeasures = metaData.Nmeasures-nsub;

	metaData.Nvis = metaData.Nmeasures*metaData.freq_per_chan*metaData.polarisations;



  // calculate the coordinates of the center

	double uvshift = metaData.uvmin/(metaData.uvmax-metaData.uvmin);

	//printf("UVSHIFT %f %f %f %f %f\n",metaData.uvmin, metaData.uvmax, metaData.wmin, metaData.wmax, uvshift);



	if (rank == 0)

	{

	   printf("N. measurements %ld\n",metaData.Nmeasures);

	   printf("N. visibilities %ld\n",metaData.Nvis);

	}



  // Set temporary local size of points

	long nm_pe = (long)(metaData.Nmeasures/size);

	long remaining = metaData.Nmeasures%size;



  	long startrow = rank*nm_pe;

  	if (rank == size-1)nm_pe = nm_pe+remaining;



	long Nmeasures_tot = metaData.Nmeasures;

	metaData.Nmeasures = nm_pe;

	long Nvis_tot = metaData.Nvis;

	metaData.Nvis = metaData.Nmeasures*metaData.freq_per_chan*metaData.polarisations;

	metaData.Nweights = metaData.Nmeasures*metaData.polarisations;



#ifdef VERBOSE

	printf("N. measurements on %d %ld\n",rank,metaData.Nmeasures);

	printf("N. visibilities on %d %ld\n",rank,metaData.Nvis);

#endif





	// DAV: all these arrays can be allocatate statically for the sake of optimization. However be careful that if MPI is used

	// all the sizes are rescaled by the number of MPI tasks

	// Allocate arrays

	uu = (double*) calloc(metaData.Nmeasures,sizeof(double));

	vv = (double*) calloc(metaData.Nmeasures,sizeof(double));

	ww = (double*) calloc(metaData.Nmeasures,sizeof(double));

	weights = (float*) calloc(metaData.Nweights,sizeof(float));

	visreal = (float*) calloc(metaData.Nvis,sizeof(float));

	visimg = (float*) calloc(metaData.Nvis,sizeof(float));



  if(rank == 0)printf("READING DATA\n");

	// Read data

	fileName(datapath, in.ufile);

	//printf("Reading %s\n",filename);



	file.pFile = fopen (filename,"rb");

	fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

	fread(uu,metaData.Nmeasures*sizeof(double),1,file.pFile);

	fclose(file.pFile);



        fileName(datapath, in.vfile);

	//printf("Reading %s\n",filename);



  file.pFile = fopen (filename,"rb");

	fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

	fread(vv,metaData.Nmeasures*sizeof(double),1,file.pFile);

	fclose(file.pFile);



        fileName(datapath, in.wfile);

	//printf("Reading %s\n",filename);



	file.pFile = fopen (filename,"rb");

	fseek (file.pFile,startrow*sizeof(double),SEEK_SET);

	fread(ww,metaData.Nmeasures*sizeof(double),1,file.pFile);

	fclose(file.pFile);



#ifdef USE_MPI

	MPI_Barrier(MPI_COMM_WORLD);

#endif



	clock_gettime(CLOCK_MONOTONIC, &finish);

	end = clock();

	timing.setup_time = ((double) (end - start)) / CLOCKS_PER_SEC;

	timing.setup_time1 = (finish.tv_sec - begin.tv_sec);

	timing.setup_time1 += (finish.tv_nsec - begin.tv_nsec) / 1000000000.0;



        /*INIT function */

        init();



  if(rank == 0)printf("GRIDDING DATA\n");
@@ -218,8 +98,8 @@ if(rank == 0){ 
	for (long iphi = 0; iphi < metaData.Nmeasures; iphi++)

	{

     boundary[iphi] = -1;

	   uuh = uu[iphi];

	   vvh = vv[iphi];

	   uuh = data.uu[iphi];

	   vvh = data.vv[iphi];

	   int binphi = (int)(vvh*nsectors);

	   // check if the point influence also neighboring slabs

	   double updist = (double)((binphi+1)*yaxis)*dx - vvh;
@@ -292,13 +172,13 @@ if(rank == 0){ 


	// MAIN LOOP OVER FILES

	//

	for (int ifiles=0; ifiles<ndatasets; ifiles++)

	for (int ifiles=0; ifiles<param.ndatasets; ifiles++)

	{

	strcpy(filename,datapath_multi[ifiles]);

        printf("Processing %s, %d of %d\n",filename,ifiles+1,ndatasets);

	strcpy(filename,param.datapath_multi[ifiles]);

        printf("Processing %s, %d of %d\n",filename,ifiles+1,param.ndatasets);



        // Read metadata

        fileName(datapath_multi[ifiles], in.metafile);

        fileName(param.datapath_multi[ifiles], in.metafile);

        //readMetaData(filename);

        file.pFile = fopen (filename,"r");

        fscanf(file.pFile,"%ld",&metaData.Nmeasures0);
@@ -325,13 +205,13 @@ if(rank == 0){ 
        fileName(datapath, in.weightsfile);

        file.pFile = fopen (filename,"rb");

        fseek (file.pFile,startrow*metaData.polarisations*sizeof(float),SEEK_SET);

        fread(weights,(metaData.Nweights)*sizeof(float),1,file.pFile);

        fread(data.weights,(metaData.Nweights)*sizeof(float),1,file.pFile);

        fclose(file.pFile);

  

  fileName(datapath, in.visrealfile);

  file.pFile = fopen (filename,"rb");

  fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);

  fread(visreal,metaData.Nvis*sizeof(float),1,file.pFile);

  fread(data.visreal,metaData.Nvis*sizeof(float),1,file.pFile);

  fclose(file.pFile);

  fileName(datapath, in.visimgfile);
@@ -340,7 +220,7 @@ if(rank == 0){ 
#endif

  file.pFile = fopen (filename,"rb");

  fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);

  fread(visimg,metaData.Nvis*sizeof(float),1,file.pFile);

  fread(data.visimg,metaData.Nvis*sizeof(float),1,file.pFile);

  fclose(file.pFile);



#ifdef USE_MPI
@@ -386,21 +266,21 @@ if(rank == 0){ 
	      while (current->index != -1)

              {

                   long ilocal = current->index;

      	           //double vvh = vv[ilocal];

      	           //double vvh = data.vv[ilocal];

             	   //int binphi = (int)(vvh*nsectors);

	           //if (binphi == isector || boundary[ilocal] == isector) {

	           uus[icount] = uu[ilocal];

	           vvs[icount] = vv[ilocal]-isector*shift;

	           wws[icount] = ww[ilocal];

	           uus[icount] = data.uu[ilocal];

	           vvs[icount] = data.vv[ilocal]-isector*shift;

	           wws[icount] = data.ww[ilocal];

             	  for (long ipol=0; ipol<metaData.polarisations; ipol++)

                  {

		        weightss[ip] = weights[ilocal*metaData.polarisations+ipol];

		        weightss[ip] = data.weights[ilocal*metaData.polarisations+ipol];

                        ip++;

                  }

                 for (long ifreq=0; ifreq<metaData.polarisations*metaData.freq_per_chan; ifreq++)

                 {

	                 visreals[inu] = visreal[ilocal*metaData.polarisations*metaData.freq_per_chan+ifreq];

	                 visimgs[inu] = visimg[ilocal*metaData.polarisations*metaData.freq_per_chan+ifreq];

	                 visreals[inu] = data.visreal[ilocal*metaData.polarisations*metaData.freq_per_chan+ifreq];

	                 visimgs[inu] = data.visimg[ilocal*metaData.polarisations*metaData.freq_per_chan+ifreq];

		               //if(visimgs[inu]>1e10 || visimgs[inu]<-1e10)printf("%f %f %ld %ld %d %ld %ld\n",visreals[inu],visimgs[inu],inu,Nvissec,rank,ilocal*metaData.polarisations*metaData.freq_per_chan+ifreq,metaData.Nvis);

	                 inu++;

	         }
@@ -458,7 +338,7 @@ if(rank == 0){ 
	             xaxis,

	             yaxis,

               gridss,

               num_threads);

               param.num_threads);



/* int z =0 ;

#pragma omp target map(to:test_i_gpu) map(from:z)
@@ -798,7 +678,7 @@ if(rank == 0){ 
        double* image_real = (double*) calloc(xaxis*yaxis,sizeof(double));	

        double* image_imag = (double*) calloc(xaxis*yaxis,sizeof(double));	



        phase_correction(gridss,image_real,image_imag,xaxis,yaxis,num_w_planes,grid_size_x,grid_size_y,resolution,metaData.wmin,metaData.wmax,num_threads);

        phase_correction(gridss,image_real,image_imag,xaxis,yaxis,num_w_planes,grid_size_x,grid_size_y,resolution,metaData.wmin,metaData.wmax,param.num_threads);



        end = clock();

        clock_gettime(CLOCK_MONOTONIC, &finish);
@@ -865,7 +745,7 @@ if(rank == 0){ 
          printf("Phase time:    %f sec\n",timing.phase_time);

#endif

          printf("TOT time:      %f sec\n",timing.tot_time);

	  if(num_threads > 1)

	  if(param.num_threads > 1)

          {

            printf("PSetup time:   %f sec\n",timing.setup_time1);

            printf("PProcess time: %f sec\n",timing.process_time1);
@@ -881,7 +761,7 @@ if(rank == 0){ 
	if (rank == 0)

	{

	 file.pFile = fopen (out.timingfile,"w");

	 if (num_threads == 1)

	 if (param.num_threads == 1)

         {

	   fprintf(file.pFile, "%f %f %f %f %f %f %f\n",timing.setup_time,timing.kernel_time,timing.compose_time,timing.reduce_time,timing.fftw_time,timing.phase_time,timing.tot_time);

	 } else {

proto.h

+3 −1
Original line number Diff line number Diff line
@@ -3,7 +3,9 @@ 


/* init.c */



void init(int argc, char * argv[]);

void init();

void read_parameter_file(char *);

void fileName(char datapath[900], char file[30]);

void readMetaData(char fileLocal[1000]);

void allocate_memory();

void readData();