Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

utilcomm_Cart.cpp

Go to the documentation of this file.
00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 #include "emdata.h"
00004 
00005 #include "ali3d_unified_mpi.h"
00006 #include "utilcomm_Cart.h"
00007 #include "utilcomm.h"
00008 
00009 using namespace EMAN;
00010 
00011 int ReadStackandDist_Cart(MPI_Comm comm_2d, EMData ***images2D, char *stackfname, int *nloc)
00012 {
00013     int ncpus, ierr, mpierr=0;
00014     MPI_Status mpistatus;
00015     EMUtil *my_util;
00016     int nima; 
00017     FILE *fp=NULL;
00018     int ROW = 0, COL = 1;
00019     int my2dpid, mycoords[2], dims[2], periods[2];
00020     int srpid, srcoords[2]; // Send/receive info
00021 
00022     // Get dims and my coordinates
00023     MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00024     MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
00025     ncpus = dims[ROW]*dims[COL];
00026 
00027     ierr = 0;
00028     if (mycoords[ROW] == 0 && mycoords[COL] == 0) {
00029         fp = fopen(stackfname,"r");
00030         if (!fp) {
00031             ierr = 1;
00032             printf("failed to open %s\n", stackfname);
00033         }
00034         else {
00035             fclose(fp);
00036             nima = my_util->get_image_count(stackfname);
00037 /**********************************************/
00038 //    nima = 84;
00039 /****************************/
00040         }
00041     }
00042     mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm_2d);
00043     if (ierr == 0) {
00044         int *psize = new int[dims[ROW]];
00045         int *nbase = new int[dims[ROW]];
00046     
00047         mpierr = MPI_Bcast(&nima, 1, MPI_INT, 0, comm_2d);
00048     
00049         *nloc = setpart_gc1(comm_2d, nima, psize, nbase);
00050         *images2D = new EMData*[*nloc]; // NB!: whoever calls ReadStackandDist must delete this!
00051 
00052         EMData *img_ptr;
00053         int img_index;
00054         float *imgdata;
00055 
00056         // read the first image to get size
00057         img_ptr = new EMData();
00058         img_ptr->read_image(stackfname, 0);
00059         int nx = img_ptr->get_xsize();
00060         int ny = img_ptr->get_ysize();
00061 
00062         float s2x, s2y;
00063 
00064         if (mycoords[COL] == 0 && mycoords[ROW] == 0) {
00065             printf("Master node reading and distributing %d images along the first column of processors\n", nima);
00066             for ( int ip = 0 ; ip < dims[ROW] ; ++ip ) {
00067                 for ( int i = 0 ; i < psize[ip] ; ++i ) {
00068                     img_index = nbase[ip] + i;
00069 
00070                     if (ip != 0) {
00071                         img_ptr->read_image(stackfname, img_index);
00072                         // get a pointer to the image's data
00073                         imgdata = img_ptr->get_data();
00074                         // find the x/y shift values if it has them, otherwise set them to 0.0
00075                         try {
00076                             s2x = (*images2D)[i]->get_attr("s2x");
00077                         } catch ( std::exception& e ) {
00078                             s2x = 0.0;
00079                         }
00080                         try {
00081                             s2y = (*images2D)[i]->get_attr("s2y");
00082                         } catch ( std::exception& e ) {
00083                             s2y = 0.0;
00084                         }
00085                         // send these to processor (ip,0)
00086                         srcoords[COL] = 0;
00087                         srcoords[ROW] = ip;
00088                         MPI_Cart_rank(comm_2d, srcoords, &srpid);
00089 
00090                         MPI_Send(imgdata, nx*ny, MPI_FLOAT, srpid, srpid, comm_2d);
00091                         MPI_Send(&s2x, 1, MPI_FLOAT, srpid, srpid, comm_2d);
00092                         MPI_Send(&s2y, 1, MPI_FLOAT, srpid, srpid, comm_2d);
00093                     } else { // ip == 0                             
00094                         (*images2D)[i] = new EMData();
00095                         (*images2D)[i]->read_image(stackfname, img_index);
00096                         try {
00097                             s2x = (*images2D)[i]->get_attr("s2x");
00098                         } catch ( std::exception& e ) {
00099                             (*images2D)[i]->set_attr("s2x",0.0);
00100                         }
00101                         try {
00102                             s2y = (*images2D)[i]->get_attr("s2y");
00103                         } catch ( std::exception& e ) {
00104                             (*images2D)[i]->set_attr("s2y",0.0);
00105                         }
00106                     }
00107                 }
00108                 //printf("finished reading data for processor %d\n", ip);
00109             }
00110 
00111         } else if (mycoords[COL] == 0 && mycoords[ROW] != 0) {
00112             for ( int i = 0 ; i < psize[mycoords[ROW]] ; ++i ) {
00113                 (*images2D)[i] = new EMData();
00114                 (*images2D)[i]->set_size(nx, ny, 1);
00115                 imgdata = (*images2D)[i]->get_data();
00116                 
00117                 MPI_Recv(imgdata, nx*ny, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
00118                 MPI_Recv(&s2x, 1, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
00119                 MPI_Recv(&s2y, 1, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
00120                 (*images2D)[i]->set_attr("s2x",s2x);
00121                 (*images2D)[i]->set_attr("s2y",s2y);
00122             }
00123             //printf("received %d images for processor (%d,%d)\n", *nloc, mycoords[ROW], mycoords[COL]);
00124         }
00125 
00126 // Now have all the processors in group g_c_0 broadcast the images and angles along the row communicator
00127     if (mycoords[COL] == 0 && mycoords[ROW] == 0)
00128       printf("First column of processors distributing images to other columns..\n");
00129 
00130     if (mycoords[COL] == 0) {
00131       for (int img_index = 0; img_index < *nloc; img_index++){
00132         imgdata = (*images2D)[img_index]->get_data();
00133         
00134         // find the x/y shift values if it has them, otherwise set them to 0.0
00135         try {
00136             s2x = (*images2D)[img_index]->get_attr("s2x");
00137         } catch ( std::exception& e ) {
00138             s2x = 0.0;
00139         }
00140         try {
00141             s2y = (*images2D)[img_index]->get_attr("s2y");
00142         } catch ( std::exception& e ) {
00143             s2y = 0.0;
00144         }
00145 
00146         for(int ip=1; ip< dims[COL]; ip++){
00147 // printf("Proc (%d, %d) sending image %d of %d to Proc (%d,%d) \n", mycoords[ROW], mycoords[COL], img_index, *nloc, mycoords[ROW], ip); 
00148           srcoords[ROW] = mycoords[ROW];
00149           srcoords[COL] = ip;
00150           MPI_Cart_rank(comm_2d, srcoords, &srpid);
00151 
00152           MPI_Send(imgdata, nx*ny, MPI_FLOAT, srpid, img_index, comm_2d);
00153           MPI_Send(&s2x, 1, MPI_FLOAT, srpid, img_index, comm_2d);
00154           MPI_Send(&s2y, 1, MPI_FLOAT, srpid, img_index, comm_2d);
00155         }
00156       }
00157     }
00158     if (mycoords[COL] != 0) {
00159   //      printf("Proc (%d, %d) receiving...", mycoords[ROW], mycoords[COL]);
00160         for(int img_index = 0; img_index < *nloc; img_index++){
00161           (*images2D)[img_index] = new EMData();
00162           (*images2D)[img_index]->set_size(nx, ny, 1);
00163           imgdata = (*images2D)[img_index]->get_data();
00164           
00165           srcoords[ROW] = mycoords[ROW];
00166           srcoords[COL] = 0;
00167           MPI_Cart_rank(comm_2d, srcoords, &srpid);
00168 
00169           MPI_Recv(imgdata, nx*ny, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);
00170           MPI_Recv(&s2x, 1, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);
00171           MPI_Recv(&s2y, 1, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);
00172 
00173           (*images2D)[img_index]->set_attr("s2x",s2x);
00174           (*images2D)[img_index]->set_attr("s2y",s2y);
00175          }
00176     }
00177 
00178         ierr = mpierr; 
00179 
00180         EMDeletePtr(img_ptr);
00181         EMDeleteArray(psize);
00182         EMDeleteArray(nbase);
00183     }
00184     return ierr;
00185 }
00186 
00187 int ReadAngTrandDist_Cart(MPI_Comm comm_2d, MPI_Comm comm_row, int *dims, float *angleshift, char *angfname, int nloc)
00188 {
00189    int ierr = 0;
00190    int ROW = 0, COL = 1;
00191    int my2dpid, mycoords[2];
00192    int srpid, srcoords[2], keep_dims[2];
00193    MPI_Status mpistatus;
00194    FILE *fp = NULL;
00195 
00196    int nimgs=0;
00197 
00198    MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
00199    MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords); // Get my coordinates
00200 
00201    float * iobuffer   = new float[5*nloc];
00202    if (!iobuffer) {
00203       fprintf(stderr,"failed to allocate buffer to read angles shifts\n");
00204       ierr = -1;
00205       goto EXIT;
00206    }
00207 
00208    if (mycoords[COL] == 0 && mycoords[ROW] == 0) { //I am Proc (0,0)
00209       fp = fopen(angfname,"r");
00210       if (!fp)  ierr = 1;
00211    }
00212    MPI_Bcast(&ierr, 1, MPI_INT, 0, comm_2d);
00213 
00214    if ( ierr ) {
00215       if (mycoords[COL] == 0 && mycoords[ROW] == 0) 
00216           fprintf(stderr,"failed to open %s\n", angfname);
00217       ierr = MPI_Finalize();
00218       goto EXIT;
00219    }
00220    else {
00221        if (mycoords[COL] == 0 && mycoords[ROW] == 0) { //I am Proc (0,0)
00222           for (int iproc = 0; iproc < dims[ROW]; iproc++) {
00223              // figure out the number of images assigned to processor (iproc,0)
00224              if (iproc > 0) {
00225                 srcoords[COL] = 0;
00226                 srcoords[ROW] = iproc;
00227                 MPI_Cart_rank(comm_2d, srcoords, &srpid);
00228 
00229                 MPI_Recv(&nimgs, 1, MPI_INT, srpid, srpid, comm_2d, &mpistatus);
00230 
00231                 // Read the next nimgs set of angles and shifts
00232                 for (int i = 0; i < nimgs; i++) {
00233                    fscanf(fp,"%f %f %f %f %f", 
00234                           &iobuffer[5*i+0],
00235                           &iobuffer[5*i+1],
00236                           &iobuffer[5*i+2],
00237                           &iobuffer[5*i+3],
00238                           &iobuffer[5*i+4]);
00239                 }
00240                 MPI_Send(iobuffer,5*nimgs,MPI_FLOAT,srpid,srpid,comm_2d);
00241              }
00242              else {
00243                 for (int i = 0; i < nloc; i++) {
00244                    fscanf(fp,"%f %f %f %f %f", 
00245                           &angleshift[5*i+0],
00246                           &angleshift[5*i+1],
00247                           &angleshift[5*i+2],
00248                           &angleshift[5*i+3],
00249                           &angleshift[5*i+4]);
00250                 }
00251              }
00252           }
00253           fclose(fp);
00254        }
00255        else if (mycoords[COL] == 0 && mycoords[ROW] != 0) { //I am in the first column
00256           // send image count to the master processor (mypid = 0)
00257 
00258           MPI_Send(&nloc, 1, MPI_INT, 0, my2dpid, comm_2d);
00259           // Receive angleshifts
00260           MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
00261        }
00262   }
00263 
00264   // Now have all the processors in group g_c_0 broadcast the angles along the row communicator
00265   srcoords[ROW] = 0;
00266   MPI_Cart_rank(comm_row, srcoords, &srpid);
00267   MPI_Bcast(angleshift, 5*nloc, MPI_FLOAT, srpid, comm_row);
00268 
00269   EMDeleteArray(iobuffer);
00270 
00271 EXIT:
00272    return ierr;
00273 }
00274 
00275 int CleanStack_Cart(MPI_Comm comm_col, EMData ** image_stack, int nloc, int ri, Vec3i volsize, Vec3i origin)
00276 {
00277     int nx = volsize[0];
00278     int ny = volsize[1];
00279         
00280     float * rhs = new float[nx*ny];
00281     float * imgdata;
00282 
00283     // Calculate average "background" from all pixels strictly outside of radius
00284     double aba, abaloc; // average background
00285     aba = 0.0;
00286     abaloc = 0.0;
00287     int klp, klploc; // number of pixels in background
00288     klp = 0;
00289     klploc = 0;
00290 
00291     // calculate avg background in parallel
00292     for ( int i = 0 ; i < nloc ; ++i ) {
00293         imgdata = image_stack[i]->get_data();
00294         asta2(imgdata, nx, ny, ri, &abaloc, &klploc);
00295     }
00296         
00297     // All reduce using the column-based sub-topology
00298 
00299     MPI_Allreduce(&abaloc, &aba, 1, MPI_DOUBLE, MPI_SUM, comm_col);
00300     MPI_Allreduce(&klploc, &klp, 1, MPI_INT, MPI_SUM, comm_col);
00301 
00302     aba /= klp;
00303 
00304     // subtract off the average background from pixels weakly inside of radius
00305     int x_summand, y_summand;
00306     int r_squared = ri * ri;
00307     for ( int i = 0 ; i < nloc ; ++i ) {
00308         imgdata = image_stack[i]->get_data();
00309         for ( int j = 0 ; j < nx ; ++j) {
00310             x_summand = (j - origin[0]) *  (j - origin[0]);
00311             for ( int k = 0 ; k < ny ; ++k ) {
00312                 y_summand = (k - origin[1]) *  (k - origin[1]);
00313                 if ( x_summand + y_summand <= r_squared) {
00314                     imgdata[j*ny + k] -= aba;
00315                 }
00316             }
00317         }
00318     }
00319     EMDeleteArray(rhs);
00320     return 0;
00321 }
00322 
00323 int setpart_gc1(MPI_Comm comm_2d, int nangs, int *psize, int *nbase)
00324 /* This function provides the partition for the set of 2D images along the first column of processors
00325 
00326         Input:
00327             comm_2d - Cartesian communicator
00328               nangs - number of 2D images
00329         Output:
00330               psize - vector containing the partition distribution
00331               nbase - vector containing the base of the partitions
00332            nangsloc - number of local angles
00333 */
00334 {
00335         int ROW = 0, COL = 1;
00336         int dims[2], periods[2], mycoords[2];
00337         int nangsloc, nrem;
00338 
00339         // Get information associated with comm_2d
00340         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00341  
00342         nangsloc = nangs/dims[ROW];
00343         nrem = nangs - dims[ROW]*nangsloc;
00344         if (mycoords[ROW] < nrem) nangsloc++;
00345 
00346         for (int i = 0; i < dims[ROW]; i++) {
00347                 psize[i] = nangs/dims[ROW];
00348                 if (i < nrem) psize[i]++;
00349         }
00350  
00351         for (int i = 0; i < dims[ROW]; i++) {
00352           if(i==0)
00353             nbase[i] = 0;
00354           else
00355             nbase[i] = nbase[i-1] + psize[i-1];
00356         }
00357    
00358 //printf("I am [%d,%d], nloc = %d, psize = [%d, %d], nbase = [%d, %d]", mycoords[ROW], mycoords[COL],nangsloc,psize[0],psize[1], nbase[0], nbase[1]);
00359 
00360         return nangsloc;
00361 }
00362 
00363 
00364 int setpart_gr1(MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase)
00365 /* This function provides the ideal partition for nnz in the 3D volume along the first row of processors.
00366 
00367         Input:
00368             comm_2d - Cartesian communicator
00369                 nnz - total number of non-zeros in the volume
00370         Output:
00371             nnzpart - vector containing the partition distribution
00372             nnzbase - vector containing the base of the partitions
00373              nnzloc - initial local nnz
00374 */
00375 {
00376         int ROW = 0, COL = 1;
00377         int dims[2], periods[2], mycoords[2];
00378         int nnzloc, nrem;
00379         
00380         // Get information associated with comm_2d
00381         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00382  
00383         nnzloc = nnz/dims[COL];
00384         nrem = nnz - dims[COL]*nnzloc;
00385         if (mycoords[COL] < nrem) nnzloc++;
00386 
00387         for (int i = 0; i < dims[COL]; i++) {
00388                 nnzpart[i] = nnz/dims[COL];
00389                 if (i < nrem) nnzpart[i]++;
00390         }
00391  
00392         nnzbase[0] = 0; 
00393         for (int i = 1; i < dims[COL]; i++) {
00394                 nnzbase[i] = nnzbase[i-1] + nnzpart[i-1];
00395         }
00396         nnzbase[dims[COL]] = nnz;
00397         return nnzloc;
00398 }
00399 
00400 

Generated on Tue Jun 11 13:40:45 2013 for EMAN2 by  doxygen 1.3.9.1