diff --git a/gridding.c b/gridding.c
index 55660cc087e10d71df5c7f231b07c82d02f986ff..758317749388280f2f09acd60aee4924b3998c09 100644
--- a/gridding.c
+++ b/gridding.c
@@ -36,39 +36,44 @@ void initialize_array(){
histo_send = (long*) calloc(nsectors+1,sizeof(long));
int * boundary = (int*) calloc(metaData.Nmeasures,sizeof(int));
- double uuh,vvh;
+ double vvh;
+
for (long iphi = 0; iphi < metaData.Nmeasures; iphi++)
- {
- boundary[iphi] = -1;
- vvh = data.vv[iphi]; //less or equal to 0.6
- int binphi = (int)(vvh*nsectors); //has values expect 0 and nsectors-1. So we use updist and downdist condition
- // 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]++;
- if(updist < w_supporth && updist >= 0.0) {histo_send[binphi+1]++; boundary[iphi] = binphi+1;};
- if(downdist < w_supporth && binphi > 0 && downdist >= 0.0) {histo_send[binphi-1]++; boundary[iphi] = binphi-1;};
- }
+ {
+ boundary[iphi] = -1;
+ vvh = data.vv[iphi]; //less or equal to 0.6
+ int binphi = (int)(vvh*nsectors); //has values expect 0 and nsectors-1. So we use updist and downdist condition
+ // 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]++;
+ if(updist < w_supporth && updist >= 0.0)
+ {histo_send[binphi+1]++; boundary[iphi] = binphi+1;};
+ if(downdist < w_supporth && binphi > 0 && downdist >= 0.0)
+ {histo_send[binphi-1]++; boundary[iphi] = binphi-1;};
+ }
sectorarray = (long**)malloc ((nsectors+1) * sizeof(long*));
for(int sec=0; sec<(nsectors+1); sec++)
- {
- sectorarray[sec] = (long*)malloc(histo_send[sec]*sizeof(long));
- }
-
+ {
+ sectorarray[sec] = (long*)malloc(histo_send[sec]*sizeof(long));
+ }
+
long *counter = (long*) calloc(nsectors+1,sizeof(long));
for (long iphi = 0; iphi < metaData.Nmeasures; iphi++)
- {
- vvh = data.vv[iphi];
- int binphi = (int)(vvh*nsectors);
- double updist = (double)((binphi+1)*yaxis)*dx - vvh;
- double downdist = vvh - (double)(binphi*yaxis)*dx;
- sectorarray[binphi][counter[binphi]] = iphi;
- counter[binphi]++;
- if(updist < w_supporth && updist >= 0.0) { sectorarray[binphi+1][counter[binphi+1]] = iphi; counter[binphi+1]++;};
- if(downdist < w_supporth && binphi > 0 && downdist >= 0.0) { sectorarray[binphi-1][counter[binphi-1]] = iphi; counter[binphi-1]++;};
- }
+ {
+ vvh = data.vv[iphi];
+ int binphi = (int)(vvh*nsectors);
+ double updist = (double)((binphi+1)*yaxis)*dx - vvh;
+ double downdist = vvh - (double)(binphi*yaxis)*dx;
+ sectorarray[binphi][counter[binphi]] = iphi;
+ counter[binphi]++;
+ if(updist < w_supporth && updist >= 0.0)
+ { sectorarray[binphi+1][counter[binphi+1]] = iphi; counter[binphi+1]++;};
+ if(downdist < w_supporth && binphi > 0 && downdist >= 0.0)
+ { sectorarray[binphi-1][counter[binphi-1]] = iphi; counter[binphi-1]++;};
+ }
#ifdef VERBOSE
@@ -235,11 +240,11 @@ void gridding_data(){
int target_rank = (int)(isector % size);
#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_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);
+
// for every task, gridss coincides with the
// that can be avoided if shared window coincides with gridss
memcpy(Me.win.ptr+isector*size_of_grid, gridss, size_of_grid*sizeof(double));
@@ -247,58 +252,59 @@ void gridding_data(){
reduce( isector, target_rank ); // here the reduce is performed within every host
-
- // here thre reduce is performed among hosts
- MPI_Barrier(MPI_COMM_WORLD);
-
- int Im_in_the_new_communicator = MPI_UNDEFINED;
- if(global_rank == target_rank)
- Im_in_the_new_communicator = 1;
- else
- if( Me.Rank[HOSTS] == 0 )
- {
- if( Me.Ranks_to_host[ target_rank ] != Me.myhost )
- Im_in_the_new_communicator = 1;
- }
-
- MPI_Comm Sector_Comm;
- MPI_Comm_split( COMM[WORLD], Im_in_the_new_communicator, global_rank, &Sector_Comm);
-
- if( Sector_Comm != MPI_COMM_NULL )
+ if ( Me.Nhosts > 1 )
{
- int sector_size;
- int sector_rank = 0;
- int sector_target;
+ // here the reduce is performed among hosts
+ MPI_Barrier(MPI_COMM_WORLD);
- MPI_Comm_size( Sector_Comm, §or_size);
- MPI_Comm_rank( Sector_Comm, §or_rank);
- if ( global_rank == target_rank)
- {
- MPI_Send( §or_rank, 1, MPI_INT, 0, 0, Sector_Comm);
- memcpy(grid, Me.swins[Me.Rank[myHOST]].ptr+isector*size_of_grid*sizeof(double), size_of_grid * sizeof(double));
- }
+ int Im_in_the_new_communicator = MPI_UNDEFINED;
+ if(global_rank == target_rank)
+ Im_in_the_new_communicator = 1;
+ else
+ if( Me.Rank[HOSTS] == 0 )
+ {
+ if( Me.Ranks_to_host[ target_rank ] != Me.myhost )
+ Im_in_the_new_communicator = 1;
+ }
- if( sector_rank == 0 )
+ MPI_Comm Sector_Comm;
+ MPI_Comm_split( COMM[WORLD], Im_in_the_new_communicator, global_rank, &Sector_Comm);
+
+ if( Sector_Comm != MPI_COMM_NULL )
{
- MPI_Status status;
- MPI_Recv( §or_target, 1, MPI_INT, MPI_ANY_SOURCE, 0, Sector_Comm, &status);
+ int sector_size;
+ int sector_rank = 0;
+ int sector_target;
+
+ MPI_Comm_size( Sector_Comm, §or_size);
+ MPI_Comm_rank( Sector_Comm, §or_rank);
+ if ( global_rank == target_rank)
+ {
+ MPI_Send( §or_rank, 1, MPI_INT, 0, 0, Sector_Comm);
+ memcpy(grid, Me.swins[Me.Rank[myHOST]].ptr+isector*size_of_grid*sizeof(double), size_of_grid * sizeof(double));
+ }
+
+ if( sector_rank == 0 )
+ {
+ MPI_Status status;
+ MPI_Recv( §or_target, 1, MPI_INT, MPI_ANY_SOURCE, 0, Sector_Comm, &status);
+ }
+
+ MPI_Bcast( §or_target, 1, MPI_INT, 0, Sector_Comm );
+
+ MPI_Reduce(grid, grid, size_of_grid, MPI_DOUBLE,MPI_SUM, sector_target, Sector_Comm);
+
+ MPI_Comm_free( &Sector_Comm );
}
-
- MPI_Bcast( §or_target, 1, MPI_INT, 0, Sector_Comm );
-
- MPI_Reduce(grid, grid, size_of_grid, MPI_DOUBLE,MPI_SUM, sector_target, Sector_Comm);
-
- MPI_Comm_free( &Sector_Comm );
}
-
+ #else // relates to #ifdef ONE_SIDE
- //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
+ #endif // closes #ifdef ONE_SIDE
+ #endif // closes USE_MPI
clock_gettime(CLOCK_MONOTONIC, &finishk);
endk = clock();
@@ -430,6 +436,10 @@ void reduce( int sector, int target_rank )
int local_rank = Me.Rank[myHOST];
if( Me.Ranks_to_host[ target_rank ] == Me.myhost )
+ // exchange rank 0 with target rank
+ // in this way the following log2 alogorithm,
+ // which reduces to rank 0, will work for
+ // every target rank
{
int r = 0;
diff --git a/init.c b/init.c
index 919e6feb6268ff6b33831c6aea1c3e10eb50201c..8312a3b439f98e8d8c6faed57ded03c89e2e9c99 100644
--- a/init.c
+++ b/init.c
@@ -132,120 +132,120 @@ void read_parameter_file(char *fname)
{
if( (file.pFile = fopen (fname,"r")) != NULL )
{
- char buf1[30], buf2[100], buf3[30], num[30];
- int i = 1;
- while(fscanf(file.pFile, "%s" "%s", buf1, buf2) != EOF)
- {
- if(strcmp(buf1, "num_threads") == 0)
- {
- param.num_threads = atoi(buf2);
- }
- if(strcmp(buf1, "Datapath1") == 0)
- {
- strcpy(param.datapath_multi[0], buf2);
- i++;
- }
- if(strcmp(buf1, "ndatasets") == 0)
- {
- param.ndatasets = atoi(buf2);
- }
- if(strcmp(buf1, "w_support") == 0)
- {
- param.w_support = atoi(buf2);
- }
- if(strcmp(buf1, "grid_size_x") == 0)
- {
- param.grid_size_x = atoi(buf2);
- }
- if(strcmp(buf1, "grid_size_y") == 0)
- {
- param.grid_size_y = atoi(buf2);
- }
- if(strcmp(buf1, "num_w_planes") == 0)
- {
- param.num_w_planes = atoi(buf2);
- }
- if(strcmp(buf1, "ufile") == 0)
- {
- strcpy(in.ufile, buf2);
- }
- if(strcmp(buf1, "vfile") == 0)
- {
- strcpy(in.vfile, buf2);
- }
- if(strcmp(buf1, "wfile") == 0)
- {
- strcpy(in.wfile, buf2);
- }
- if(strcmp(buf1, "weightsfile") == 0)
- {
- strcpy(in.weightsfile, buf2);
- }
- if(strcmp(buf1, "visrealfile") == 0)
- {
- strcpy(in.visrealfile, buf2);
- }
- if(strcmp(buf1, "visimgfile") == 0)
- {
- strcpy(in.visimgfile, buf2);
- }
- if(strcmp(buf1, "metafile") == 0)
- {
- strcpy(in.metafile, buf2);
- }
- if(strcmp(buf1, "outfile") == 0)
- {
- strcpy(outparam.outfile, buf2);
- }
- if(strcmp(buf1, "outfile1") == 0)
- {
- strcpy(outparam.outfile1, buf2);
- }
- if(strcmp(buf1, "outfile2") == 0)
- {
- strcpy(outparam.outfile2, buf2);
- }
- if(strcmp(buf1, "outfile3") == 0)
- {
- strcpy(outparam.outfile3, buf2);
- }
- if(strcmp(buf1, "fftfile") == 0)
- {
- strcpy(outparam.fftfile, buf2);
- }
- if(strcmp(buf1, "fftfile2") == 0)
- {
- strcpy(outparam.fftfile2, buf2);
- }
- if(strcmp(buf1, "fftfile3") == 0)
- {
- strcpy(outparam.fftfile3, buf2);
- }
- if(strcmp(buf1, "logfile") == 0)
- {
- strcpy(outparam.logfile, buf2);
- }
- if(strcmp(buf1, "extension") == 0)
- {
- strcpy(outparam.extension, buf2);
- }
- if(strcmp(buf1, "timingfile") == 0)
- {
- strcpy(outparam.timingfile, buf2);
- }
- if(param.ndatasets > 1)
- {
+ char buf1[30], buf2[100], buf3[30], num[30];
+ int i = 1;
+ while(fscanf(file.pFile, "%s" "%s", buf1, buf2) != EOF)
+ {
+ if(strcmp(buf1, "num_threads") == 0)
+ {
+ param.num_threads = atoi(buf2);
+ }
+ if(strcmp(buf1, "Datapath1") == 0)
+ {
+ strcpy(param.datapath_multi[0], buf2);
+ i++;
+ }
+ if(strcmp(buf1, "ndatasets") == 0)
+ {
+ param.ndatasets = atoi(buf2);
+ }
+ if(strcmp(buf1, "w_support") == 0)
+ {
+ param.w_support = atoi(buf2);
+ }
+ if(strcmp(buf1, "grid_size_x") == 0)
+ {
+ param.grid_size_x = atoi(buf2);
+ }
+ if(strcmp(buf1, "grid_size_y") == 0)
+ {
+ param.grid_size_y = atoi(buf2);
+ }
+ if(strcmp(buf1, "num_w_planes") == 0)
+ {
+ param.num_w_planes = atoi(buf2);
+ }
+ if(strcmp(buf1, "ufile") == 0)
+ {
+ strcpy(in.ufile, buf2);
+ }
+ if(strcmp(buf1, "vfile") == 0)
+ {
+ strcpy(in.vfile, buf2);
+ }
+ if(strcmp(buf1, "wfile") == 0)
+ {
+ strcpy(in.wfile, buf2);
+ }
+ if(strcmp(buf1, "weightsfile") == 0)
+ {
+ strcpy(in.weightsfile, buf2);
+ }
+ if(strcmp(buf1, "visrealfile") == 0)
+ {
+ strcpy(in.visrealfile, buf2);
+ }
+ if(strcmp(buf1, "visimgfile") == 0)
+ {
+ strcpy(in.visimgfile, buf2);
+ }
+ if(strcmp(buf1, "metafile") == 0)
+ {
+ strcpy(in.metafile, buf2);
+ }
+ if(strcmp(buf1, "outfile") == 0)
+ {
+ strcpy(outparam.outfile, buf2);
+ }
+ if(strcmp(buf1, "outfile1") == 0)
+ {
+ strcpy(outparam.outfile1, buf2);
+ }
+ if(strcmp(buf1, "outfile2") == 0)
+ {
+ strcpy(outparam.outfile2, buf2);
+ }
+ if(strcmp(buf1, "outfile3") == 0)
+ {
+ strcpy(outparam.outfile3, buf2);
+ }
+ if(strcmp(buf1, "fftfile") == 0)
+ {
+ strcpy(outparam.fftfile, buf2);
+ }
+ if(strcmp(buf1, "fftfile2") == 0)
+ {
+ strcpy(outparam.fftfile2, buf2);
+ }
+ if(strcmp(buf1, "fftfile3") == 0)
+ {
+ strcpy(outparam.fftfile3, buf2);
+ }
+ if(strcmp(buf1, "logfile") == 0)
+ {
+ strcpy(outparam.logfile, buf2);
+ }
+ if(strcmp(buf1, "extension") == 0)
+ {
+ strcpy(outparam.extension, buf2);
+ }
+ if(strcmp(buf1, "timingfile") == 0)
+ {
+ strcpy(outparam.timingfile, buf2);
+ }
+ if(param.ndatasets > 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);
+ 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);
}
else
@@ -278,18 +278,20 @@ void readMetaData(char fileLocal[1000]) {
{
if( (file.pFile = fopen (fileLocal,"r")) != NULL )
{
- 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);
+ int items = 0;
+ items += fscanf(file.pFile,"%ld",&metaData.Nmeasures);
+ items += fscanf(file.pFile,"%ld",&metaData.Nvis);
+ items += fscanf(file.pFile,"%ld",&metaData.freq_per_chan);
+ items += fscanf(file.pFile,"%ld",&metaData.polarisations);
+ items += fscanf(file.pFile,"%ld",&metaData.Ntimes);
+ items += fscanf(file.pFile,"%lf",&metaData.dt);
+ items += fscanf(file.pFile,"%lf",&metaData.thours);
+ items += fscanf(file.pFile,"%ld",&metaData.baselines);
+ items += fscanf(file.pFile,"%lf",&metaData.uvmin);
+ items += fscanf(file.pFile,"%lf",&metaData.uvmax);
+ items += fscanf(file.pFile,"%lf",&metaData.wmin);
+ items += fscanf(file.pFile,"%lf",&metaData.wmax);
+
fclose(file.pFile);
}
else
@@ -315,14 +317,15 @@ void metaData_calculation() {
metaData.Nvis = metaData.Nmeasures*metaData.freq_per_chan*metaData.polarisations;
// calculate the coordinates of the center
+ #warning why it it not used?
double uvshift = metaData.uvmin/(metaData.uvmax-metaData.uvmin);
if (global_rank == 0)
- {
- printf("N. measurements %ld\n",metaData.Nmeasures);
- printf("N. visibilities %ld\n",metaData.Nvis);
- }
-
+ {
+ 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;
@@ -330,16 +333,14 @@ void metaData_calculation() {
startrow = global_rank*nm_pe;
if (global_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",global_rank,metaData.Nmeasures);
- printf("N. visibilities on %d %ld\n",global_rank,metaData.Nvis);
- #endif
+ #ifdef VERBOSE
+ printf("N. measurements on %d %ld\n",global_rank,metaData.Nmeasures);
+ printf("N. visibilities on %d %ld\n",global_rank,metaData.Nvis);
+ #endif
}
@@ -391,7 +392,9 @@ void readData() {
if( (file.pFile = fopen (filename,"rb")) != NULL )
{
fseek (file.pFile,startrow*sizeof(double),SEEK_SET);
- fread(data.uu,metaData.Nmeasures*sizeof(double),1,file.pFile);
+ size_t res = fread(data.uu, sizeof(double), metaData.Nmeasures, file.pFile);
+ if( res != metaData.Nmeasures )
+ printf("an error occurred while reading file %s\n", filename);
fclose(file.pFile);
}
else
@@ -405,7 +408,10 @@ void readData() {
if( (file.pFile = fopen (filename,"rb")) != NULL )
{
fseek (file.pFile,startrow*sizeof(double),SEEK_SET);
- fread(data.vv,metaData.Nmeasures*sizeof(double),1,file.pFile);
+ size_t res = fread(data.vv, sizeof(double), metaData.Nmeasures, file.pFile);
+ if( res != metaData.Nmeasures )
+ printf("an error occurred while reading file %s\n", filename);
+
fclose(file.pFile);
}
else
@@ -419,7 +425,9 @@ void readData() {
if( (file.pFile = fopen (filename,"rb")) != NULL )
{
fseek (file.pFile,startrow*sizeof(double),SEEK_SET);
- fread(data.ww,metaData.Nmeasures*sizeof(double),1,file.pFile);
+ size_t res = fread(data.ww,sizeof(double), metaData.Nmeasures, file.pFile);
+ if( res != metaData.Nmeasures )
+ printf("an error occurred while reading file %s\n", filename);
fclose(file.pFile);
}
else
@@ -433,7 +441,9 @@ void readData() {
if( (file.pFile = fopen (filename,"rb")) != NULL)
{
fseek (file.pFile,startrow*metaData.polarisations*sizeof(float),SEEK_SET);
- fread(data.weights,(metaData.Nweights)*sizeof(float),1,file.pFile);
+ size_t res = fread(data.weights, sizeof(float), metaData.Nweights, file.pFile);
+ if( res != metaData.Nweights )
+ printf("an error occurred while reading file %s\n", filename);
fclose(file.pFile);
}
else
@@ -447,7 +457,9 @@ void readData() {
if((file.pFile = fopen (filename,"rb")) != NULL )
{
fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);
- fread(data.visreal,metaData.Nvis*sizeof(float),1,file.pFile);
+ size_t res = fread(data.visreal, sizeof(float), metaData.Nvis, file.pFile);
+ if( res != metaData.Nvis )
+ printf("an error occurred while reading file %s\n", filename);
fclose(file.pFile);
}
else
@@ -461,7 +473,9 @@ void readData() {
if( (file.pFile = fopen (filename,"rb")) != NULL )
{
fseek (file.pFile,startrow*metaData.freq_per_chan*metaData.polarisations*sizeof(float),SEEK_SET);
- fread(data.visimg,metaData.Nvis*sizeof(float),1,file.pFile);
+ size_t res = fread(data.visimg, sizeof(float), metaData.Nvis, file.pFile);
+ if( res != metaData.Nvis )
+ printf("an error occurred while reading file %s\n", filename);
fclose(file.pFile);
}
else
diff --git a/w-stacking.cu b/w-stacking.cu
index 4fe36254789fff354e28fb887a5d4bb6c3b37704..a74718808658643b56602666ec5ad2ebbc2f855a 100644
--- a/w-stacking.cu
+++ b/w-stacking.cu
@@ -140,7 +140,7 @@ void wstack(
int num_threads)
{
long i;
- long index;
+ //long index;
long visindex;
// initialize the convolution kernel
@@ -242,7 +242,7 @@ void wstack(
visindex = i*freq_per_chan*polarizations;
- double sum = 0.0;
+ //double sum = 0.0;
int j, k;
//if (i%1000 == 0)printf("%ld\n",i);