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

utilcomm_Cart.h File Reference

#include "mpi.h"
#include "stdlib.h"
#include "emdata.h"
#include "alignoptions.h"

Include dependency graph for utilcomm_Cart.h:

Include dependency graph

This graph shows which files directly or indirectly include this file:

Included by dependency graph

Go to the source code of this file.

Functions

int ReadStackandDist_Cart (MPI_Comm comm, EMData ***images2D, char *stackfname, int *nloc)
int CleanStack_Cart (MPI_Comm comm, EMData **image_stack, int nloc, int ri, Vec3i volsize, Vec3i origin)
int ReadAngTrandDist_Cart (MPI_Comm comm_2d, MPI_Comm comm_row, int *dim, float *angleshift, char *angfname, int nloc)
int setpart_gc1 (MPI_Comm comm_2d, int nangs, int *psize, int *nbase)
int setpart_gr1 (MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase)
int sphpart (MPI_Comm comm_2d, int nrays, int *ptrs, int *nnzbase, int *ptrstart)
int getcb2sph (Vec3i volsize, int ri, Vec3i origin, int nnz0, int *ptrs, int *cord)
int fwdpj3_Cart (Vec3i volsize, int nraysloc, int nnzloc, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)
int bckpj3_Cart (Vec3i volsize, int nraysloc, int nnzloc, float *dm, Vec3i origin, int ri, int *ptrs, int *cord, int myptrstart, float *x, float *y)


Function Documentation

int bckpj3_Cart Vec3i  volsize,
int  nraysloc,
int  nnzloc,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
int  myptrstart,
float *  x,
float *  y
 

Definition at line 237 of file project3d_Cart.cpp.

References cord, dm, ifix(), myptrstart, nx, ny, ptrs, status, EMAN::Vec3i, x, and y.

00240              :
00241           volsize - size of 3D cube volume
00242           nraysloc - local number of rays
00243           nnzloc - local number of voxels
00244           dm - projection matrix
00245           origin, ri - origin and radius of sphere
00246           ptrs, cord - pointers and coordinates for first ray
00247           myptrstart - determines starting rays
00248           x - 2D image
00249         Output:
00250           y - portion of the 3D volume
00251 */
00252 {
00253     int       i, j, iqx,iqy, xc, yc, zc, jj;
00254     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
00255     int       status = 0; 
00256 
00257     int xcent = origin[0];
00258     int ycent = origin[1];
00259     int zcent = origin[2];
00260 
00261     int nx = volsize[0];
00262     int ny = volsize[1];
00263 
00264     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00265     float sx, sy;
00266 
00267     sx = dm(7);
00268     sy = dm(8);
00269 
00270     if ( nx > 2*ri) {
00271         for (i = 1; i <= nraysloc; i++) {
00272             zc = cord(1,myptrstart+i) - zcent;
00273             yc = cord(2,myptrstart+i) - ycent;
00274             xc = cord(3,myptrstart+i) - xcent;
00275 
00276             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
00277             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;
00278 
00279             for (j = ptrs(myptrstart+i); j <ptrs(myptrstart+i+1); j++) {
00280                 jj = j-ptrs[myptrstart]+1;
00281 
00282                 // jj is the index on the local volume
00283                 iqx = ifix(xb);
00284                 iqy = ifix(yb);
00285 
00286                 dx = xb - iqx;
00287                 dy = yb - iqy;
00288                 dx1m = 1.0 - dx;
00289                 dy1m = 1.0 - dy;
00290                 dxdy = dx*dy;
00291 /*
00292 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
00293 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
00294 c     &                     + dx  *dy1m*x(iqx+1, iqy)
00295 c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
00296 c
00297 c              --- faster version of the above commented out
00298 c                  code (derived by summing the following table 
00299 c                  of coefficients along  the colunms) ---
00300 c
00301 c                        1         dx        dy      dxdy
00302 c                     ------   --------  --------  -------
00303 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
00304 c                                        x(i,j+1) -x(i,j+1)
00305 c                              x(i+1,j)           -x(i+1,j)
00306 c                                                x(i+1,j+1) 
00307 c
00308 */
00309                 // Phi: add index checking, now that shifts are being used
00310                 if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
00311                     y(jj) += x(iqx,iqy);
00312                     if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
00313                         y(jj) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
00314                     }
00315                     if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
00316                         y(jj) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
00317                     }
00318                     if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
00319                         y(jj) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00320                     }
00321                 }
00322 
00323 //                y(j) += x(iqx,iqy)
00324 //                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
00325 //                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
00326 //                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
00327 //                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
00328 
00329                xb += dm(1);
00330                yb += dm(4);
00331             } // end for j
00332         } // end for i
00333      }
00334     else {
00335         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
00336     }
00337 
00338     return status;
00339 }

int CleanStack_Cart MPI_Comm  comm,
EMData **  image_stack,
int  nloc,
int  ri,
Vec3i  volsize,
Vec3i  origin
 

Definition at line 275 of file utilcomm_Cart.cpp.

References asta2(), EMDeleteArray(), EMAN::EMData::get_data(), imgdata, nx, ny, rhs, and EMAN::Vec3i.

Referenced by main().

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 }

int fwdpj3_Cart Vec3i  volsize,
int  nraysloc,
int  nnzloc,
float *  dm,
Vec3i  origin,
int  ri,
int *  ptrs,
int *  cord,
int  myptrstart,
float *  x,
float *  y
 

Definition at line 147 of file project3d_Cart.cpp.

References cord, dm, ifix(), myptrstart, nx, ny, ptrs, status, EMAN::Vec3i, x, and y.

00150              :
00151           volsize - size (nx,ny,nz) of 3D cube volume
00152           nraysloc - local number of rays within the compact spherical representation
00153           nnzloc - local number of voxels
00154           dm - an array of size 9 storing transformation 
00155           origin, ri - origin and radius of sphere
00156           ptrs, cord - pointers and coordinates for first ray
00157           myptrstart - determines starting rays
00158           x - portion of the 3D input volume
00159         Output:
00160           y - 2D output image
00161 */
00162 {
00163     int    iqx, iqy, i, j, xc, yc, zc, jj;
00164     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
00165     int    status = 0;
00166     
00167     // Phi: adding the shift parameters that get passed in as the last two entries of dm
00168     float sx, sy;
00169 
00170     sx = dm(7);
00171     sy = dm(8);
00172 
00173     int xcent = origin[0];
00174     int ycent = origin[1];
00175     int zcent = origin[2];
00176 
00177     int nx = volsize[0];
00178     int ny = volsize[1];
00179 
00180     dm1 = dm(1);
00181     dm4 = dm(4);
00182  
00183     if ( nx > 2*ri ) {
00184         for (i = 1; i <= nraysloc; i++) {
00185 
00186             zc = cord(1,myptrstart+i)-zcent;
00187             yc = cord(2,myptrstart+i)-ycent;
00188             xc = cord(3,myptrstart+i)-xcent;
00189             xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
00190             yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;
00191 
00192             for (j = ptrs(myptrstart+i); j< ptrs(myptrstart+i+1); j++) {
00193                jj = j-ptrs[myptrstart]+1;
00194                iqx = ifix(xb);
00195                iqy = ifix(yb);
00196 
00197                 // jj is the index on the local volume
00198                ct   = x(jj);
00199 
00200                // dipx =  xb - (float)(iqx);
00201                // dipy = (yb - (float)(iqy)) * ct;
00202                    dipx =  xb - iqx;
00203                    dipy = (yb - iqy) * ct;
00204 
00205                dipy1m = ct - dipy;
00206                dipx1m = 1.0 - dipx;
00207 
00208                         if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
00209                // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
00210                y(iqx  ,iqy)   +=  dipx1m*dipy1m;
00211                         if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
00212                // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
00213                y(iqx+1,iqy)   +=  dipx*dipy1m; 
00214                         if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
00215                // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
00216                y(iqx+1,iqy+1) +=  dipx*dipy;         
00217                         if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
00218                // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
00219                y(iqx  ,iqy+1) +=  dipx1m*dipy;
00220                xb += dm1;
00221                yb += dm4;
00222            }
00223         }
00224     }
00225     else {
00226         fprintf(stderr, " nx must be greater than 2*ri\n");
00227         exit(1);
00228     }
00229     return status;
00230 }
#undef x

int getcb2sph Vec3i  volsize,
int  ri,
Vec3i  origin,
int  nnz0,
int *  ptrs,
int *  cord
 

Definition at line 15 of file project3d_Cart.cpp.

References cord, nnz, nrays, nx, ny, ptrs, status, and EMAN::Vec3i.

00018              :
00019          volsize - dimensions of 3D cube
00020               ri - radius of sphere
00021           origin - origin of sphere
00022             nnz0 - number of non-zeros
00023 
00024         Output:
00025               ptrs - vector containing pointers to the beginning of each ray
00026               cord - vector containing the x,y,z coordinates of the beginning of each ray
00027 */
00028 
00029 {
00030     int    xs, ys, zs, xx, yy, zz, rs, r2;
00031     int    ix, iy, iz, jnz, nnz, nrays;
00032     int    ftm = 0, status = 0;  
00033 
00034     int xcent = (int)origin[0];
00035     int ycent = (int)origin[1];
00036     int zcent = (int)origin[2];
00037 
00038     int nx = (int)volsize[0];
00039     int ny = (int)volsize[1];
00040     int nz = (int)volsize[2];
00041 
00042     r2      = ri*ri;
00043     nnz     = 0;
00044     nrays    = 0;
00045     ptrs(1) = 1;
00046 
00047     for (ix = 1; ix <= nx; ix++) {
00048        xs  = ix-xcent;
00049        xx  = xs*xs;
00050        for ( iy = 1; iy <= ny; iy++ ) {
00051            ys = iy-ycent;
00052            yy = ys*ys;
00053            jnz = 0;
00054 
00055            ftm = 1;
00056            // not the most efficient implementation
00057            for (iz = 1; iz <= nz; iz++) {
00058                zs = iz-zcent;
00059                zz = zs*zs;
00060                rs = xx + yy + zz;
00061                if (rs <= r2) {
00062                   jnz++;
00063                   nnz++;
00064      //             sphere(nnz) = cube(iz, iy, ix); 
00065 
00066                   //  record the coordinates of the first nonzero ===
00067                   if (ftm) {
00068                      nrays++;
00069                      cord(1,nrays) = iz; 
00070                      cord(2,nrays) = iy; 
00071                      cord(3,nrays) = ix;
00072                      ftm = 0;
00073                   }
00074                }
00075             } // end for (iz..)
00076             if (jnz > 0) {
00077                 ptrs(nrays+1) = ptrs(nrays) + jnz;
00078             }  // endif (jnz)
00079        } // end for iy
00080     } // end for ix
00081     if (nnz != nnz0) status = -1;
00082     return status;
00083 }

int ReadAngTrandDist_Cart MPI_Comm  comm_2d,
MPI_Comm  comm_row,
int *  dim,
float *  angleshift,
char *  angfname,
int  nloc
 

Definition at line 187 of file utilcomm_Cart.cpp.

References COL, EMDeleteArray(), ierr, and ROW.

Referenced by main().

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 }

int ReadStackandDist_Cart MPI_Comm  comm,
EMData ***  images2D,
char *  stackfname,
int *  nloc
 

Definition at line 11 of file utilcomm_Cart.cpp.

References COL, EMDeleteArray(), EMDeletePtr(), EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::EMUtil::get_image_count(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), ierr, img_ptr, imgdata, nx, ny, EMAN::EMData::read_image(), ROW, and setpart_gc1().

Referenced by main().

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 }

int setpart_gc1 MPI_Comm  comm_2d,
int  nangs,
int *  psize,
int *  nbase
 

Definition at line 323 of file utilcomm_Cart.cpp.

00326              :
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 }

int setpart_gr1 MPI_Comm  comm_2d,
int  nnz,
int *  nnzpart,
int *  nnzbase
 

Definition at line 364 of file utilcomm_Cart.cpp.

00367              :
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 }

int sphpart MPI_Comm  comm_2d,
int  nrays,
int *  ptrs,
int *  nnzbase,
int *  ptrstart
 

Definition at line 85 of file project3d_Cart.cpp.

00088              :
00089             comm_2d - Cartesian communicator
00090               nrays - total number of rays
00091                ptrs - vector containing pointers
00092             nnzbase - ideal volume partition of nnz
00093         Output:
00094            ptrstart - vector containing all the starting pointers for each column processor group
00095            nraysloc - actual number of local rays
00096 */
00097 {
00098         int ROW = 0, COL = 1;
00099         int dims[2], periods[2], mycoords[2];
00100         int nraysloc = 0;
00101         
00102         // Get information associated with comm_2d
00103         MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00104 
00105         int count = 1;
00106         if (mycoords[COL] == 0){ // First column starts out with the first ray
00107           nraysloc++;
00108         }
00109         ptrstart[0] = 0;
00110 
00111         for(int i=1; i<nrays ; i++){
00112           if (ptrs[i]-1 <= nnzbase[count] && ptrs[i+1]-1 >= nnzbase[count]){ 
00113                 //nnzbase is between or equal to ptrs
00114 
00115             if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
00116               if(mycoords[COL] == count-1){  // ray belongs to count-1
00117                 nraysloc++;
00118               }
00119               ptrstart[count] = i+1;
00120               count++;
00121 
00122             } else { //nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]
00123                 if(mycoords[COL] == count){// ray belongs to count and it's a new starting ray
00124                   nraysloc++;
00125                 }
00126                 ptrstart[count] = i;
00127                 count++;
00128             }
00129 
00130         }
00131         else {  //ptrs[i]-1 > nnzbase[count] so ray belongs to count-1
00132           if(mycoords[COL] == count-1){
00133             nraysloc++;
00134           }
00135 
00136         }
00137         } // end for
00138         ptrstart[dims[COL]] = nrays;
00139         return nraysloc;
00140 
00141 }


Generated on Tue Jun 11 13:47:32 2013 for EMAN2 by  doxygen 1.3.9.1