Full working CFITSIO parallel fits writing implementation

w-stacking-fftw.c

@@ -122,8 +122,8 @@ int main(int argc, char * argv[])
 	int num_threads;
 
         // Resolution 
-	double dx = 0.5/(double)grid_size_x;
-	double dw = 0.5/(double)num_w_planes;
+	double dx = 1.0/(double)grid_size_x;
+	double dw = 1.0/(double)num_w_planes;
 	// Half support size
 	double w_supporth = (double)((w_support-1)/2)*dx;
 
@@ -133,9 +133,9 @@ int main(int argc, char * argv[])
 	fitsfile *fptrimg;
 	int status;
 	long nelements;
-	long fpixel, lpixel;
-	char testfitsreal[FILENAMELENGTH] = "parallel_np2_real.fits";
-	char testfitsimag[FILENAMELENGTH] = "parallel_np2_img.fits";
+//	long fpixel, lpixel;
+	char testfitsreal[FILENAMELENGTH] = "parallel_np4_real.fits";
+	char testfitsimag[FILENAMELENGTH] = "parallel_np4_img.fits";
 
 	long naxis = 2;
 	long naxes[2] = { grid_size_x, grid_size_y };
@@ -217,7 +217,6 @@ if(rank == 0){
 	// LOCAL grid size
 	xaxis = local_grid_size_x;
 	yaxis = local_grid_size_y;
-
 	clock_gettime(CLOCK_MONOTONIC, &begin);
 	start = clock();
 
@@ -975,33 +974,38 @@ if(rank == 0){
 	    fits_close_file(fptreal, &status);
 	    #endif
 
-//            pFilereal = fopen (fftfile2,"wb");
-//            pFileimg = fopen (fftfile3,"wb");
-
-//            fclose(pFilereal);
-//            fclose(pFileimg);
+            pFilereal = fopen (fftfile2, "wb");
+	    pFileimg = fopen (fftfile3, "wb");
+	    fclose(pFilereal);
+	    fclose(pFileimg);
         }
 	#ifdef USE_MPI
 	MPI_Barrier(MPI_COMM_WORLD);
         #endif
         if(rank == 0)printf("WRITING IMAGE\n");
+        long * fpixel = (long *) malloc(sizeof(long)*naxis);
+        long * lpixel = (long *) malloc(sizeof(long)*naxis);
+        
 	for (int isector=0; isector<size; isector++)
 	{
 	    #ifdef USE_MPI
 	    MPI_Barrier(MPI_COMM_WORLD);
             #endif
 	    if(isector == rank)
+            //if(rank == 0)
 	    {
                #ifdef FITSIO
 
 	       printf("%d writing\n",isector);
-               long * fpixel = (long *) malloc(sizeof(long)*naxis);
-               long * lpixel = (long *) malloc(sizeof(long)*naxis);
+               //long * fpixel = (long *) malloc(sizeof(long)*naxis);
+               //long * lpixel = (long *) malloc(sizeof(long)*naxis);
                
                fpixel[0] = 1;
-               lpixel[0] = xaxis;
                fpixel[1] = isector*yaxis+1;
+               lpixel[0] = xaxis;
                lpixel[1] = (isector+1)*yaxis;
+               //printf("fpixel %d, %d\n", fpixel[0], fpixel[1]);
+	       //printf("lpixel %d, %d\n", lpixel[0], lpixel[1]);
 	       
                status = 0;
 	       fits_open_image(&fptreal, testfitsreal, READWRITE, &status);
@@ -1015,18 +1019,18 @@ if(rank == 0){
 
                #endif
 
-//               pFilereal = fopen (fftfile2,"ab");
-//               pFileimg = fopen (fftfile3,"ab");
+               pFilereal = fopen (fftfile2,"ab");
+               pFileimg = fopen (fftfile3,"ab");
 
-//	       long global_index = isector*(xaxis*yaxis)*sizeof(double);
+	       long global_index = isector*(xaxis*yaxis)*sizeof(double);
 
-//               fseek(pFilereal, global_index, SEEK_SET);
-//               fwrite(image_real, xaxis*yaxis, sizeof(double), pFilereal);
-//               fseek(pFileimg, global_index, SEEK_SET);
-//               fwrite(image_imag, xaxis*yaxis, sizeof(double), pFileimg);
+               fseek(pFilereal, global_index, SEEK_SET);
+               fwrite(image_real, xaxis*yaxis, sizeof(double), pFilereal);
+               fseek(pFileimg, global_index, SEEK_SET);
+               fwrite(image_imag, xaxis*yaxis, sizeof(double), pFileimg);
 
-//               fclose(pFilereal);
-//               fclose(pFileimg);
+               fclose(pFilereal);
+               fclose(pFileimg);
 	    }
 	}
 	#ifdef USE_MPI