From 0894606b48dc47182e2fd1734e18fea95c301c53 Mon Sep 17 00:00:00 2001 From: nandhanas Date: Wed, 2 Mar 2022 04:35:51 +0100 Subject: [PATCH] Restructing of i/o --- Makefile | 14 +- allvars.c | 17 +- allvars.h | 21 +- data/paramfile.txt | 13 + fourier_transform.c | 213 +++++++++++++ gridding.c | 398 +++++++++++++++++++++++ init.c | 327 +++++++++++++------ main.c | 753 +++----------------------------------------- proto.h | 20 +- result.c | 52 +++ 10 files changed, 994 insertions(+), 834 deletions(-) create mode 100644 data/paramfile.txt create mode 100644 fourier_transform.c create mode 100644 gridding.c create mode 100644 result.c diff --git a/Makefile b/Makefile index 2a3f21c..6b529a5 100644 --- a/Makefile +++ b/Makefile @@ -39,8 +39,8 @@ OPT += -DWRITE_IMAGE OPT += -DPHASE_ON -DEPS = w-stacking.h main.c w-stacking.cu phase_correction.cu allvars.h init.c -COBJ = w-stacking.o main.o phase_correction.o allvars.o init.o +DEPS = w-stacking.h main.c w-stacking.cu phase_correction.cu allvars.h init.c gridding.c fourier_transform.c result.c +COBJ = w-stacking.o main.o phase_correction.o allvars.o init.o gridding.o fourier_transform.o result.o w-stacking.c: w-stacking.cu cp w-stacking.cu w-stacking.c @@ -60,17 +60,17 @@ serial: $(COBJ) $(CC) $(OPTIMIZE) $(OPT) -o w-stackingCfftw_serial $^ $(LIBS) serial_omp: phase_correction.c - $(CC) $(OPTIMIZE) $(OPT) -o w-stackingOMP_serial main.c init.c w-stacking_omp.c $(CFLAGS) $(LIBS) + $(CC) $(OPTIMIZE) $(OPT) -o w-stackingOMP_serial main.c init.c gridding.c fourier_transform.c result.c w-stacking_omp.c $(CFLAGS) $(LIBS) simple_mpi: phase_correction.c - $(MPICC) $(OPTIMIZE) $(OPT) -o w-stackingMPI_simple w-stacking_omp.c main.c init.c phase_correction.c $(CFLAGS) $(LIBS) + $(MPICC) $(OPTIMIZE) $(OPT) -o w-stackingMPI_simple w-stacking_omp.c main.c init.c gridding.c fourier_transform.c result.c phase_correction.c $(CFLAGS) $(LIBS) mpi_omp: phase_correction.c - $(MPICC) $(OPTIMIZE) $(OPT) -o w-stackingMPI_omp w-stacking_omp.c main.c init.c phase_correction.c $(CFLAGS) $(LIBS) + $(MPICC) $(OPTIMIZE) $(OPT) -o w-stackingMPI_omp w-stacking_omp.c main.c init.c gridding.c fourier_transform.c result.c phase_correction.c $(CFLAGS) $(LIBS) serial_cuda: $(NVCC) $(NVFLAGS) -c w-stacking.cu phase_correction.cu $(NVLIB) - $(CC) $(OPTIMIZE) $(OPT) -c main.c init.c $(CFLAGS) $(LIBS) + $(CC) $(OPTIMIZE) $(OPT) -c main.c init.c gridding.c fourier_transform.c result.c $(CFLAGS) $(LIBS) $(CXX) $(OPTIMIZE) $(OPT) -o w-stackingfftw_serial w-stacking-fftw.o w-stacking.o phase_correction.o $(CFLAGS) $(NVLIB) -lm mpi: $(COBJ) @@ -78,7 +78,7 @@ mpi: $(COBJ) mpi_cuda: $(NVCC) $(NVFLAGS) -c w-stacking.cu phase_correction.cu $(NVLIB) - $(MPICC) $(OPTIMIZE) $(OPT) -c main.c init.c$(CFLAGS) $(LIBS) + $(MPICC) $(OPTIMIZE) $(OPT) -c main.c init.c fourier_transform.c result.c $(CFLAGS) $(LIBS) $(MPIC++) $(OPTIMIZE) $(OPT) -o w-stackingfftw w-stacking-fftw.o w-stacking.o phase_correction.o $(NVLIB) $(CFLAGS) $(LIBS) clean: diff --git a/allvars.c b/allvars.c index cb7dad6..4015037 100644 --- a/allvars.c +++ b/allvars.c @@ -4,14 +4,16 @@ struct io file; struct ip in; -struct op out = {"grid.txt", "coords.txt", "grid_real.bin", "grid_img.bin", "fft.txt", "fft_real.bin", "fft_img.bin", "run.log", ".txt", "timings.dat"}; +struct op out; +struct op outparam = {"grid.txt", "coords.txt", "grid_real.bin", "grid_img.bin", "fft.txt", "fft_real.bin", "fft_img.bin", "run.log", ".txt", "timings.dat"}; struct meta metaData; struct time timing; struct parameter param; struct fileData data; -char filename[1000]; +int w_support = 7; +char filename[1000], buf[30], num_buf[30]; char datapath[900]; int xaxis, yaxis; int grid_size_x = 2048; @@ -21,7 +23,16 @@ int rank; int size; long nsectors; long startrow; -double resolution; +double resolution, dx, dw, w_supporth; clock_t start, end, start0, startk, endk; struct timespec begin, finish, begin0, begink, finishk; + +struct sectorlist ** sectorhead; +long * histo_send; +double * grid, *gridss; +fftw_complex *fftwgrid; + +#ifdef USE_MPI + MPI_Win slabwin; +#endif diff --git a/allvars.h b/allvars.h index dcbf6ee..31528e3 100644 --- a/allvars.h +++ b/allvars.h @@ -61,7 +61,7 @@ extern struct op char extension[30]; char timingfile[30]; -} out; +} out, outparam; extern struct meta { @@ -108,17 +108,30 @@ extern struct fileData float * visimg; }data; -extern char filename[1000]; +extern struct sectorlist { + long index; + struct sectorlist * next; +}** sectorhead; + + +extern char filename[1000], buf[30], num_buf[30]; extern char datapath[900]; extern int xaxis, yaxis; extern int grid_size_x; extern int grid_size_y; -extern int num_w_planes; +extern int num_w_planes, w_support; extern int rank; extern int size; extern long nsectors; extern long startrow; -extern double resolution; +extern double resolution, dx, dw, w_supporth; extern clock_t start, end, start0, startk, endk; extern struct timespec begin, finish, begin0, begink, finishk; +extern long * histo_send; +extern double * grid, *gridss; +extern fftw_complex *fftwgrid; + +#ifdef USE_MPI + extern MPI_Win slabwin; +#endif diff --git a/data/paramfile.txt b/data/paramfile.txt new file mode 100644 index 0000000..4127543 --- /dev/null +++ b/data/paramfile.txt @@ -0,0 +1,13 @@ +ndatasets 3 +Datapath1 /u/nsakthivel/LOFAR/data/newgauss2noconj_t201806301100_SBL180.binMS/ +Datapath2 /u/nsakthivel/LOFAR/data/newgauss2noconj_t201806301100_SBL180.binMS/ +Datapath3 /u/nsakthivel/LOFAR/data/newgauss2noconj_t201806301100_SBL180.binMS/ +num_threads 2 +ufile ucoord.bin +vfile vcoord.bin +wfile wcoord.bin +weightsfile weights.bin +visrealfile visibilities_real.bin +visimgfile visibilities_img.bin +metafile meta.txt + diff --git a/fourier_transform.c b/fourier_transform.c new file mode 100644 index 0000000..130cd9d --- /dev/null +++ b/fourier_transform.c @@ -0,0 +1,213 @@ +#include +#include "allvars.h" +#include "proto.h" + +void fftw_data(){ + + #ifdef USE_FFTW + // FFT transform the data (using distributed FFTW) + if(rank == 0)printf("PERFORMING FFT\n"); + clock_gettime(CLOCK_MONOTONIC, &begin); + start = clock(); + fftw_plan plan; + ptrdiff_t alloc_local, local_n0, local_0_start; + double norm = 1.0/(double)(grid_size_x*grid_size_y); + + // map the 1D array of complex visibilities to a 2D array required by FFTW (complex[*][2]) + // x is the direction of contiguous data and maps to the second parameter + // y is the parallelized direction and corresponds to the first parameter (--> n0) + // and perform the FFT per w plane + alloc_local = fftw_mpi_local_size_2d(grid_size_y, grid_size_x, MPI_COMM_WORLD,&local_n0, &local_0_start); + fftwgrid = fftw_alloc_complex(alloc_local); + plan = fftw_mpi_plan_dft_2d(grid_size_y, grid_size_x, fftwgrid, fftwgrid, MPI_COMM_WORLD, FFTW_BACKWARD, FFTW_ESTIMATE); + + long fftwindex = 0; + long fftwindex2D = 0; + for (int iw=0; iw 1) + { + for (int iw=0; iw +#include "allvars.h" +#include "proto.h" + +void Push(struct sectorlist** headRef, long data) { + struct sectorlist* newNode = malloc(sizeof(struct sectorlist)); + newNode->index = data; + newNode->next = *headRef; + *headRef = newNode; +} + +void gridding(){ + + if(rank == 0)printf("GRIDDING DATA\n"); + + // Create histograms and linked lists + + clock_gettime(CLOCK_MONOTONIC, &begin); + start = clock(); + + // Initialize linked list + initialize_list(); + + //Sector and Gridding data + gridding_data(); + + #ifdef USE_MPI + MPI_Barrier(MPI_COMM_WORLD); + #endif + + end = clock(); + clock_gettime(CLOCK_MONOTONIC, &finish); + timing.process_time = ((double) (end - start)) / CLOCKS_PER_SEC; + timing.process_time1 = (finish.tv_sec - begin.tv_sec); + timing.process_time1 += (finish.tv_nsec - begin.tv_nsec) / 1000000000.0; + clock_gettime(CLOCK_MONOTONIC, &begin); + +} + +void initialize_list(){ + + sectorhead = (struct sectorlist **) malloc((nsectors+1) * sizeof(struct sectorlist)); + for (int isec=0; isec<=nsectors; isec++) + { + sectorhead[isec] = malloc(sizeof(struct sectorlist)); + sectorhead[isec]->index = -1; + sectorhead[isec]->next = NULL; + } + + + histo_send = (long*) calloc(nsectors+1,sizeof(long)); + int * boundary = (int*) calloc(metaData.Nmeasures,sizeof(int)); + double uuh,vvh; + for (long iphi = 0; iphi < metaData.Nmeasures; iphi++) + { + boundary[iphi] = -1; + 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; + double downdist = vvh - (double)(binphi*yaxis)*dx; + // + histo_send[binphi]++; + Push(§orhead[binphi],iphi); + if(updist < w_supporth && updist >= 0.0) {histo_send[binphi+1]++; boundary[iphi] = binphi+1; Push(§orhead[binphi+1],iphi);}; + if(downdist < w_supporth && binphi > 0 && downdist >= 0.0) {histo_send[binphi-1]++; boundary[iphi] = binphi-1; Push(§orhead[binphi-1],iphi);}; + } + #ifdef PIPPO + struct sectorlist * current; + long iiii = 0; + for (int j=0; jindex != -1) + { + printf("%d %d %ld %ld %ld\n",rank,j,iiii,histo_send[j],current->index); + current = current->next; + iiii++; + } + } + #endif + + #ifdef VERBOSE + for (int iii=0; iiiindex != -1) + { + long ilocal = current->index; + //double vvh = data.vv[ilocal]; + //int binphi = (int)(vvh*nsectors); + //if (binphi == isector || boundary[ilocal] == isector) { + uus[icount] = data.uu[ilocal]; + vvs[icount] = data.vv[ilocal]-isector*shift; + wws[icount] = data.ww[ilocal]; + for (long ipol=0; ipol1e10 || 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++; + } + icount++; + current = current->next; + } + + clock_gettime(CLOCK_MONOTONIC, &finishk); + endk = clock(); + timing.compose_time += ((double) (endk - startk)) / CLOCKS_PER_SEC; + timing.compose_time1 += (finishk.tv_sec - begink.tv_sec); + timing.compose_time1 += (finishk.tv_sec - begink.tv_sec); + timing.compose_time1 += (finishk.tv_nsec - begink.tv_nsec) / 1000000000.0; + + #ifndef USE_MPI + double vvmin = 1e20; + double uumax = -1e20; + double vvmax = -1e20; + + for (long ipart=0; ipart %f\n",gridss[iii]); + + #ifndef USE_MPI + long stride = isector*2*xaxis*yaxis*num_w_planes; + for (long iii=0; iii<2*xaxis*yaxis*num_w_planes; iii++)gridtot[stride+iii] = gridss[iii]; + #endif + + // Write grid in the corresponding remote slab + #ifdef USE_MPI + int target_rank = (int)isector; + //int target_rank = (int)(size-isector-1); + #ifdef ONE_SIDE + printf("One Side communication active\n"); + MPI_Win_lock(MPI_LOCK_SHARED,target_rank,0,slabwin); + MPI_Accumulate(gridss,size_of_grid,MPI_DOUBLE,target_rank,0,size_of_grid,MPI_DOUBLE,MPI_SUM,slabwin); + MPI_Win_unlock(target_rank,slabwin); + //MPI_Put(gridss,size_of_grid,MPI_DOUBLE,target_rank,0,size_of_grid,MPI_DOUBLE,slabwin); + #else + MPI_Reduce(gridss,grid,size_of_grid,MPI_DOUBLE,MPI_SUM,target_rank,MPI_COMM_WORLD); + #endif //ONE_SIDE + #endif //USE_MPI + + clock_gettime(CLOCK_MONOTONIC, &finishk); + endk = clock(); + timing.reduce_time += ((double) (endk - startk)) / CLOCKS_PER_SEC; + timing.reduce_time1 += (finishk.tv_sec - begink.tv_sec); + timing.reduce_time1 += (finishk.tv_nsec - begink.tv_nsec) / 1000000000.0; + // Go to next sector + for (long inull=0; inull<2*num_w_planes*xaxis*yaxis; inull++)gridss[inull] = 0.0; + + // Deallocate all sector arrays + free(uus); + free(vvs); + free(wws); + free(weightss); + free(visreals); + free(visimgs); + // End of loop over sector + } + // Finalize MPI communication + // Finalize MPI communication + #ifdef USE_MPI + MPI_Win_fence(0,slabwin); + #endif + + #ifndef USE_MPI + fclose(file.pFile1); + #endif + + #ifdef USE_MPI + MPI_Barrier(MPI_COMM_WORLD); + #endif +} + +void write_grided_data() +{ + + #ifdef WRITE_DATA + // Write results + if (rank == 0) + { + printf("WRITING GRIDDED DATA\n"); + file.pFilereal = fopen (out.outfile2,"wb"); + file.pFileimg = fopen (out.outfile3,"wb"); + #ifdef USE_MPI + for (int isector=0; isector 1) + { + for (int iw=0; iw 1) + { + + sprintf(num, "%d", i); + strcat(strcpy(buf3,"Datapath"),num); + if(strcmp(buf1,buf3) == 0) + { + strcpy(param.datapath_multi[i-1], buf2); + i++; + } + } } fclose(file.pFile); @@ -163,91 +186,185 @@ void fileName(char datapath[900], char file[30]) { 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); - fscanf(file.pFile,"%ld",&metaData.freq_per_chan); - fscanf(file.pFile,"%ld",&metaData.polarisations); - fscanf(file.pFile,"%ld",&metaData.Ntimes); - fscanf(file.pFile,"%lf",&metaData.dt); - fscanf(file.pFile,"%lf",&metaData.thours); - fscanf(file.pFile,"%ld",&metaData.baselines); - fscanf(file.pFile,"%lf",&metaData.uvmin); - fscanf(file.pFile,"%lf",&metaData.uvmax); - fscanf(file.pFile,"%lf",&metaData.wmin); - fscanf(file.pFile,"%lf",&metaData.wmax); - fclose(file.pFile); - } + if(rank == 0) + { + if(file.pFile = fopen (fileLocal,"r")) + { + fscanf(file.pFile,"%ld",&metaData.Nmeasures); + fscanf(file.pFile,"%ld",&metaData.Nvis); + fscanf(file.pFile,"%ld",&metaData.freq_per_chan); + fscanf(file.pFile,"%ld",&metaData.polarisations); + fscanf(file.pFile,"%ld",&metaData.Ntimes); + fscanf(file.pFile,"%lf",&metaData.dt); + fscanf(file.pFile,"%lf",&metaData.thours); + fscanf(file.pFile,"%ld",&metaData.baselines); + fscanf(file.pFile,"%lf",&metaData.uvmin); + fscanf(file.pFile,"%lf",&metaData.uvmax); + fscanf(file.pFile,"%lf",&metaData.wmin); + fscanf(file.pFile,"%lf",&metaData.wmax); + fclose(file.pFile); + } + else + { + printf("error opening meta file"); + exit(1); + } + } /* Communicating the relevent parameters to the other process */ MPI_Bcast(&metaData, sizeof(struct meta), MPI_BYTE, 0, MPI_COMM_WORLD); } +void metaData_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 + +} + 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) { + 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); + if(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); + } + else + { + printf("error opening ucoord file"); + exit(1); + } 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); + if(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); + } + else + { + printf("error opening vcoord file"); + exit(1); + } 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); + if(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); + } + else + { + printf("error opening wcoord file"); + exit(1); + } + + fileName(datapath, in.weightsfile); + if(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); + } + else + { + printf("error opening weights file"); + exit(1); + } 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); + if(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); + } + else + { + printf("error opening visibilities_real file"); + exit(1); + } 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); - */ -// } + if(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); + } + else + { + printf("error opening visibilities_img file"); + exit(1); + } + + } #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); + + MPI_Bcast(data.uu, metaData.Nmeasures * sizeof(double), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(data.vv, metaData.Nmeasures * sizeof(double), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(data.ww, metaData.Nmeasures * sizeof(double), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(data.weights, metaData.Nweights * sizeof(float), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(data.visreal, metaData.Nvis * sizeof(float), MPI_BYTE, 0, MPI_COMM_WORLD); + MPI_Bcast(data.visimg, metaData.Nvis * sizeof(float), MPI_BYTE, 0, MPI_COMM_WORLD); } diff --git a/main.c b/main.c index 8a8ccb8..925444a 100644 --- a/main.c +++ b/main.c @@ -3,34 +3,17 @@ #include "proto.h" -// Linked List set-up -struct sectorlist { - long index; - struct sectorlist * next; -}; - -void Push(struct sectorlist** headRef, long data) { - struct sectorlist* newNode = malloc(sizeof(struct sectorlist)); - newNode->index = data; - newNode->next = *headRef; - *headRef = newNode; -} - - // Main Code int main(int argc, char * argv[]) { char srank[4]; - - - // DAV: the corresponding KernelLen is calculated within the wstack function. It can be anyway hardcoded for optimization - int w_support = 7; - double dx = 1.0/(double)grid_size_x; - double dw = 1.0/(double)num_w_planes; - double w_supporth = (double)((w_support-1)/2)*dx; - double elapsed; + + // DAV: the corresponding KernelLen is calculated within the wstack function. It can be anyway hardcoded for optimization + dx = 1.0/(double)grid_size_x; + dw = 1.0/(double)num_w_planes; + w_supporth = (double)((w_support-1)/2)*dx; if(argc > 1) { @@ -44,7 +27,7 @@ int main(int argc, char * argv[]) clock_gettime(CLOCK_MONOTONIC, &begin0); start0 = clock(); - + /* Initializing MPI Environment */ #ifdef USE_MPI @@ -60,7 +43,7 @@ int main(int argc, char * argv[]) size = 1; #endif #ifdef ACCOMP - if(rank == 0){ + if(rank == 0){ if (0 == omp_get_num_devices()) { printf("No accelerator found ... exit\n"); exit(255); @@ -69,711 +52,53 @@ int main(int argc, char * argv[]) #ifdef NVIDIA prtAccelInfo(); #endif - } + } #endif - /*INIT function */ - init(); - - if(rank == 0)printf("GRIDDING DATA\n"); - - // Create histograms and linked lists - clock_gettime(CLOCK_MONOTONIC, &begin); - start = clock(); - - // Initialize linked list - struct sectorlist ** sectorhead; - sectorhead = (struct sectorlist **) malloc((nsectors+1) * sizeof(struct sectorlist)); - for (int isec=0; isec<=nsectors; isec++) - { - sectorhead[isec] = malloc(sizeof(struct sectorlist)); - sectorhead[isec]->index = -1; - sectorhead[isec]->next = NULL; - } - - - long * histo_send = (long*) calloc(nsectors+1,sizeof(long)); - int * boundary = (int*) calloc(metaData.Nmeasures,sizeof(int)); - double uuh,vvh; - for (long iphi = 0; iphi < metaData.Nmeasures; iphi++) - { - boundary[iphi] = -1; - 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; - double downdist = vvh - (double)(binphi*yaxis)*dx; - // - histo_send[binphi]++; - Push(§orhead[binphi],iphi); - if(updist < w_supporth && updist >= 0.0) {histo_send[binphi+1]++; boundary[iphi] = binphi+1; Push(§orhead[binphi+1],iphi);}; - if(downdist < w_supporth && binphi > 0 && downdist >= 0.0) {histo_send[binphi-1]++; boundary[iphi] = binphi-1; Push(§orhead[binphi-1],iphi);}; - } - -#ifdef PIPPO - struct sectorlist * current; - long iiii = 0; - for (int j=0; jindex != -1) - { - printf("%d %d %ld %ld %ld\n",rank,j,iiii,histo_send[j],current->index); - current = current->next; - iiii++; - } - } -#endif - -#ifdef VERBOSE - for (int iii=0; iiiindex != -1) - { - long ilocal = current->index; - //double vvh = data.vv[ilocal]; - //int binphi = (int)(vvh*nsectors); - //if (binphi == isector || boundary[ilocal] == isector) { - uus[icount] = data.uu[ilocal]; - vvs[icount] = data.vv[ilocal]-isector*shift; - wws[icount] = data.ww[ilocal]; - for (long ipol=0; ipol1e10 || 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++; - } - icount++; - //} - current = current->next; - } - - clock_gettime(CLOCK_MONOTONIC, &finishk); - endk = clock(); - timing.compose_time += ((double) (endk - startk)) / CLOCKS_PER_SEC; - timing.compose_time1 += (finishk.tv_sec - begink.tv_sec); - timing.compose_time1 += (finishk.tv_nsec - begink.tv_nsec) / 1000000000.0; - - #ifndef USE_MPI - double uumin = 1e20; - double vvmin = 1e20; - double uumax = -1e20; - double vvmax = -1e20; - - for (long ipart=0; ipart %f\n",gridss[iii]); + /* Reading Parameter file */ - #ifndef USE_MPI - long stride = isector*2*xaxis*yaxis*num_w_planes; - for (long iii=0; iii<2*xaxis*yaxis*num_w_planes; iii++)gridtot[stride+iii] = gridss[iii]; - #endif + read_parameter_file(in.paramfile); - // Write grid in the corresponding remote slab - #ifdef USE_MPI - int target_rank = (int)isector; - //int target_rank = (int)(size-isector-1); - #ifdef ONE_SIDE - printf("One Side communication active\n"); - MPI_Win_lock(MPI_LOCK_SHARED,target_rank,0,slabwin); - MPI_Accumulate(gridss,size_of_grid,MPI_DOUBLE,target_rank,0,size_of_grid,MPI_DOUBLE,MPI_SUM,slabwin); - MPI_Win_unlock(target_rank,slabwin); - //MPI_Put(gridss,size_of_grid,MPI_DOUBLE,target_rank,0,size_of_grid,MPI_DOUBLE,slabwin); - #else - MPI_Reduce(gridss,grid,size_of_grid,MPI_DOUBLE,MPI_SUM,target_rank,MPI_COMM_WORLD); - #endif //ONE_SIDE - #endif //USE_MPI - - clock_gettime(CLOCK_MONOTONIC, &finishk); - endk = clock(); - timing.reduce_time += ((double) (endk - startk)) / CLOCKS_PER_SEC; - timing.reduce_time1 += (finishk.tv_sec - begink.tv_sec); - timing.reduce_time1 += (finishk.tv_nsec - begink.tv_nsec) / 1000000000.0; - // Go to next sector - for (long inull=0; inull<2*num_w_planes*xaxis*yaxis; inull++)gridss[inull] = 0.0; - - // Deallocate all sector arrays - free(uus); - free(vvs); - free(wws); - free(weightss); - free(visreals); - free(visimgs); - // End of loop over sectors - } - // End of loop over input files - } - - // Finalize MPI communication - #ifdef USE_MPI - MPI_Win_fence(0,slabwin); - #endif - - // Swap left and right parts APPARENTLY NOT NECESSARY - /* - for (long kswap=0; kswap 1) - { - for (int iw=0; iw n0) - // and perform the FFT per w plane - alloc_local = fftw_mpi_local_size_2d(grid_size_y, grid_size_x, MPI_COMM_WORLD,&local_n0, &local_0_start); - fftwgrid = fftw_alloc_complex(alloc_local); - plan = fftw_mpi_plan_dft_2d(grid_size_y, grid_size_x, fftwgrid, fftwgrid, MPI_COMM_WORLD, FFTW_BACKWARD, FFTW_ESTIMATE); - - long fftwindex = 0; - long fftwindex2D = 0; - for (int iw=0; iw 1) - { - for (int iw=0; iw 1) - { - printf("PSetup time: %f sec\n",timing.setup_time1); - printf("PProcess time: %f sec\n",timing.process_time1); - printf("PKernel time = %f, PArray Composition time %f, PReduce time: %f sec\n",timing.kernel_time1,timing.compose_time1,timing.reduce_time1); -#ifdef USE_FFTW - printf("PFFTW time: %f sec\n",timing.fftw_time1); - printf("PPhase time: %f sec\n",timing.phase_time1); -#endif - printf("PTOT time: %f sec\n",timing.tot_time1); - } - } + } + + //Close MPI Environment + + #ifdef USE_MPI + MPI_Finalize(); + #endif - if (rank == 0) - { - file.pFile = fopen (out.timingfile,"w"); - 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 { - fprintf(file.pFile, "%f %f %f %f %f %f %f\n",timing.setup_time1,timing.kernel_time1,timing.compose_time1,timing.reduce_time1,timing.fftw_time1,timing.phase_time1,timing.tot_time1); - } - fclose(file.pFile); - } - - // Close MPI environment - #ifdef USE_MPI - MPI_Win_fence(0,slabwin); - MPI_Win_free(&slabwin); - MPI_Finalize(); - #endif } diff --git a/proto.h b/proto.h index f408c2e..9fb9b06 100644 --- a/proto.h +++ b/proto.h @@ -3,9 +3,27 @@ /* init.c */ -void init(); +void init(int i); +void op_filename(); void read_parameter_file(char *); void fileName(char datapath[900], char file[30]); void readMetaData(char fileLocal[1000]); +void metaData_calculation(); void allocate_memory(); void readData(); + + +/* gridding.c */ + +void gridding(); +void Push(struct sectorlist** headRef, long data); +void initialize_list(); +void gridding_data(); +void write_grided_data(); + + +/* fourier_transform.c */ + +void fftw_data(); +void write_fftw_data(); +void write_result(); diff --git a/result.c b/result.c new file mode 100644 index 0000000..7a7607d --- /dev/null +++ b/result.c @@ -0,0 +1,52 @@ +#include +#include "allvars.h" +#include "proto.h" + +void write_result() { + + end = clock(); + clock_gettime(CLOCK_MONOTONIC, &finish); + timing.tot_time = ((double) (end - start0)) / CLOCKS_PER_SEC; + timing.tot_time1 = (finish.tv_sec - begin0.tv_sec); + timing.tot_time1 += (finish.tv_nsec - begin0.tv_nsec) / 1000000000.0; + + if (rank == 0) + { + printf("Setup time: %f sec\n",timing.setup_time); + printf("Process time: %f sec\n",timing.process_time); + printf("Kernel time = %f, Array Composition time %f, Reduce time: %f sec\n",timing.kernel_time,timing.compose_time,timing.reduce_time); + #ifdef USE_FFTW + printf("FFTW time: %f sec\n",timing.fftw_time); + printf("Phase time: %f sec\n",timing.phase_time); + #endif + printf("TOT time: %f sec\n",timing.tot_time); + if(param.num_threads > 1) + { + printf("PSetup time: %f sec\n",timing.setup_time1); + printf("PProcess time: %f sec\n",timing.process_time1); + printf("PKernel time = %f, PArray Composition time %f, PReduce time: %f sec\n",timing.kernel_time1,timing.compose_time1,timing.reduce_time1); + #ifdef USE_FFTW + printf("PFFTW time: %f sec\n",timing.fftw_time1); + printf("PPhase time: %f sec\n",timing.phase_time1); + #endif + printf("PTOT time: %f sec\n",timing.tot_time1); + } + } + + if (rank == 0) + { + file.pFile = fopen (out.timingfile,"w"); + 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 { + fprintf(file.pFile, "%f %f %f %f %f %f %f\n",timing.setup_time1,timing.kernel_time1,timing.compose_time1,timing.reduce_time1,timing.fftw_time1,timing.phase_time1,timing.tot_time1); + } + fclose(file.pFile); + } + + #ifdef USE_MPI + MPI_Win_fence(0,slabwin); + MPI_Win_free(&slabwin); + #endif +} -- GitLab