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];
00021
00022
00023 MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
00024 MPI_Comm_rank(comm_2d, &my2dpid);
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
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];
00051
00052 EMData *img_ptr;
00053 int img_index;
00054 float *imgdata;
00055
00056
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
00073 imgdata = img_ptr->get_data();
00074
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
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 {
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
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
00124 }
00125
00126
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
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
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
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);
00199 MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords);
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) {
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) {
00222 for (int iproc = 0; iproc < dims[ROW]; iproc++) {
00223
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
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) {
00256
00257
00258 MPI_Send(&nloc, 1, MPI_INT, 0, my2dpid, comm_2d);
00259
00260 MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
00261 }
00262 }
00263
00264
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
00284 double aba, abaloc;
00285 aba = 0.0;
00286 abaloc = 0.0;
00287 int klp, klploc;
00288 klp = 0;
00289 klploc = 0;
00290
00291
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
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
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
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 {
00335 int ROW = 0, COL = 1;
00336 int dims[2], periods[2], mycoords[2];
00337 int nangsloc, nrem;
00338
00339
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
00359
00360 return nangsloc;
00361 }
00362
00363
00364 int setpart_gr1(MPI_Comm comm_2d, int nnz, int *nnzpart, int *nnzbase)
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 {
00376 int ROW = 0, COL = 1;
00377 int dims[2], periods[2], mycoords[2];
00378 int nnzloc, nrem;
00379
00380
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