00001 #include "mpi.h"
00002 #include "stdlib.h"
00003 #include "emdata.h"
00004
00005 #include "ali3d_unified_mpi.h"
00006 #include "utilcomm.h"
00007
00008 using namespace EMAN;
00009
00010 int ReadVandBcast(MPI_Comm comm, EMData *volume, char *volfname)
00011 {
00012 int mypid, ierr, ndim, nx, ny, nz, mpierr=0;
00013 MPI_Comm_rank(comm,&mypid);
00014 FILE *fp=NULL;
00015
00016 ierr = 0;
00017 if (mypid == 0) {
00018
00019 fp = fopen(volfname,"r");
00020 if (!fp) {
00021 ierr = 1;
00022 printf("failed to open %s\n", volfname);
00023 }
00024 else {
00025 fclose(fp);
00026 volume->read_image(volfname);
00027 ndim = volume->get_ndim();
00028 nx = volume->get_xsize();
00029 ny = volume->get_ysize();
00030 nz = volume->get_zsize();
00031 printf("ndim = %d, nx = %d, ny = %d, nz = %d\n", ndim, nx, ny, nz);
00032 }
00033 }
00034 mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00035 if (ierr==0) {
00036 mpierr = MPI_Bcast(&nx, 1, MPI_INT, 0, comm);
00037 mpierr = MPI_Bcast(&ny, 1, MPI_INT, 0, comm);
00038 mpierr = MPI_Bcast(&nz, 1, MPI_INT, 0, comm);
00039 if (mypid !=0) volume->set_size(nx,ny,nz);
00040
00041 float * voldata = volume->get_data();
00042 mpierr = MPI_Bcast(voldata, nx*ny*nz, MPI_FLOAT, 0, comm);
00043 ierr = mpierr;
00044 }
00045 return ierr;
00046 }
00047
00048 int ReadStackandDist(MPI_Comm comm, EMData ***images2D, char *stackfname, int *nloc)
00049 {
00050 int ncpus, mypid, ierr, mpierr=0;
00051 MPI_Status mpistatus;
00052 EMUtil *my_util;
00053 int nima;
00054 FILE *fp=NULL;
00055
00056 MPI_Comm_size(comm,&ncpus);
00057 MPI_Comm_rank(comm,&mypid);
00058
00059 ierr = 0;
00060 if (mypid == 0) {
00061 fp = fopen(stackfname,"r");
00062 if (!fp) {
00063 ierr = 1;
00064 printf("failed to open %s\n", stackfname);
00065 }
00066 else {
00067 fclose(fp);
00068 nima = my_util->get_image_count(stackfname);
00069 }
00070 }
00071 mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00072 if (ierr == 0) {
00073 int *psize = new int[ncpus];
00074 int *nbase = new int[ncpus];
00075
00076 mpierr = MPI_Bcast(&nima, 1, MPI_INT, 0, comm);
00077
00078 *nloc = setpart(comm, nima, psize, nbase);
00079 *images2D = new EMData*[*nloc];
00080
00081 EMData *img_ptr;
00082 int img_index;
00083 float *imgdata;
00084
00085
00086 img_ptr = new EMData();
00087 img_ptr->read_image(stackfname, 0);
00088 int nx = img_ptr->get_xsize();
00089 int ny = img_ptr->get_ysize();
00090
00091 float s2x, s2y;
00092
00093 if (mypid == 0) {
00094 printf("Master node reading and distributing %d images...\n", nima);
00095 for ( int ip = 0 ; ip < ncpus ; ++ip ) {
00096 for ( int i = 0 ; i < psize[ip] ; ++i ) {
00097 img_index = nbase[ip] + i;
00098 if (ip != 0) {
00099 img_ptr->read_image(stackfname, img_index);
00100
00101 imgdata = img_ptr->get_data();
00102
00103 try {
00104 s2x = (*images2D)[i]->get_attr("s2x");
00105 } catch ( std::exception& e ) {
00106 s2x = 0.0;
00107 }
00108 try {
00109 s2y = (*images2D)[i]->get_attr("s2y");
00110 } catch ( std::exception& e ) {
00111 s2y = 0.0;
00112 }
00113
00114 MPI_Send(imgdata, nx*ny, MPI_FLOAT, ip, ip, comm);
00115 MPI_Send(&s2x, 1, MPI_FLOAT, ip, ip, comm);
00116 MPI_Send(&s2y, 1, MPI_FLOAT, ip, ip, comm);
00117 } else {
00118 (*images2D)[i] = new EMData();
00119 (*images2D)[i]->read_image(stackfname, img_index);
00120 try {
00121 s2x = (*images2D)[i]->get_attr("s2x");
00122 } catch ( std::exception& e ) {
00123 (*images2D)[i]->set_attr("s2x",0.0);
00124 }
00125 try {
00126 s2y = (*images2D)[i]->get_attr("s2y");
00127 } catch ( std::exception& e ) {
00128 (*images2D)[i]->set_attr("s2y",0.0);
00129 }
00130 }
00131 }
00132 printf("finished reading data for processor %d\n", ip);
00133 }
00134 } else {
00135 for ( int i = 0 ; i < psize[mypid] ; ++i ) {
00136 (*images2D)[i] = new EMData();
00137 (*images2D)[i]->set_size(nx, ny, 1);
00138 imgdata = (*images2D)[i]->get_data();
00139 MPI_Recv(imgdata, nx*ny, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00140 MPI_Recv(&s2x, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00141 MPI_Recv(&s2y, 1, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00142 (*images2D)[i]->set_attr("s2x",s2x);
00143 (*images2D)[i]->set_attr("s2y",s2y);
00144 }
00145 printf("received %d images for processor %d\n", *nloc, mypid);
00146 }
00147 if (mypid == 0) printf("finished reading and distributing data\n");
00148 ierr = mpierr;
00149
00150 EMDeletePtr(img_ptr);
00151 EMDeleteArray(psize);
00152 EMDeleteArray(nbase);
00153 }
00154 return ierr;
00155 }
00156
00157 int ReadAngTrandDist(MPI_Comm comm, float *angleshift, char *paramfname, int nloc)
00158 {
00159
00160
00161
00162 int mypid, ierr=0, mpierr=0, nimgs=0, ncpus;
00163 double t0;
00164 MPI_Status mpistatus;
00165 FILE *fp=NULL;
00166
00167 float *iobuffer = new float[5*nloc];
00168 if (!iobuffer) {
00169 fprintf(stderr,"failed to allocate buffer to read angles shifts\n");
00170 ierr = -1;
00171 goto EXIT;
00172 }
00173
00174 MPI_Comm_rank(comm,&mypid);
00175 MPI_Comm_size(comm,&ncpus);
00176
00177 t0 = MPI_Wtime();
00178
00179 if (mypid ==0) {
00180 fp = fopen(paramfname,"r");
00181 if (!fp) ierr = 1;
00182 }
00183 MPI_Bcast(&ierr, 1, MPI_INT, 0, comm);
00184
00185 if ( ierr ) {
00186 if (mypid ==0) fprintf(stderr,"failed to open %s\n", paramfname);
00187 }
00188 else {
00189 if (mypid == 0) {
00190 for (int iproc = 0; iproc < ncpus; iproc++) {
00191
00192 if (iproc > 0) {
00193 MPI_Recv(&nimgs, 1, MPI_INT, iproc, iproc, comm, &mpistatus);
00194
00195 for (int i = 0; i < nimgs; i++) {
00196 fscanf(fp,"%f %f %f %f %f",
00197 &iobuffer[5*i+0],
00198 &iobuffer[5*i+1],
00199 &iobuffer[5*i+2],
00200 &iobuffer[5*i+3],
00201 &iobuffer[5*i+4]);
00202 }
00203 MPI_Send(iobuffer,5*nimgs,MPI_FLOAT,iproc,iproc,comm);
00204 }
00205 else {
00206 for (int i = 0; i < nloc; i++) {
00207 fscanf(fp,"%f %f %f %f %f",
00208 &angleshift[5*i+0],
00209 &angleshift[5*i+1],
00210 &angleshift[5*i+2],
00211 &angleshift[5*i+3],
00212 &angleshift[5*i+4]);
00213 }
00214 }
00215 }
00216 fclose(fp);
00217 }
00218 else {
00219
00220 MPI_Send(&nloc, 1, MPI_INT, 0, mypid, comm);
00221
00222 MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, mypid, comm, &mpistatus);
00223 }
00224 }
00225 EMDeleteArray(iobuffer);
00226 if (mypid == 0)
00227 printf("I/O time for reading angles & shifts = %11.3e\n",
00228 MPI_Wtime() - t0);
00229 EXIT:
00230 return ierr;
00231 }
00232
00233 int CleanStack(MPI_Comm comm, EMData ** image_stack, int nloc, int ri, Vec3i volsize, Vec3i origin)
00234 {
00235 int nx = volsize[0];
00236 int ny = volsize[1];
00237
00238 float * rhs = new float[nx*ny];
00239 float * imgdata;
00240
00241
00242 double aba, abaloc;
00243 aba = 0.0;
00244 abaloc = 0.0;
00245 int klp, klploc;
00246 klp = 0;
00247 klploc = 0;
00248
00249
00250 for ( int i = 0 ; i < nloc ; ++i ) {
00251 imgdata = image_stack[i]->get_data();
00252 asta2(imgdata, nx, ny, ri, &abaloc, &klploc);
00253 }
00254
00255 MPI_Allreduce(&abaloc, &aba, 1, MPI_DOUBLE, MPI_SUM, comm);
00256 MPI_Allreduce(&klploc, &klp, 1, MPI_INT, MPI_SUM, comm);
00257
00258 aba /= klp;
00259
00260
00261 int x_summand, y_summand;
00262 int r_squared = ri * ri;
00263 for ( int i = 0 ; i < nloc ; ++i ) {
00264 imgdata = image_stack[i]->get_data();
00265 for ( int j = 0 ; j < nx ; ++j) {
00266 x_summand = (j - origin[0]) * (j - origin[0]);
00267 for ( int k = 0 ; k < ny ; ++k ) {
00268 y_summand = (k - origin[1]) * (k - origin[1]);
00269 if ( x_summand + y_summand <= r_squared) {
00270 imgdata[j*ny + k] -= aba;
00271 }
00272 }
00273 }
00274 }
00275 EMDeleteArray(rhs);
00276 return 0;
00277 }
00278
00279 int setpart(MPI_Comm comm, int nang, int *psize, int *nbase)
00280 {
00281 int ncpus, mypid, nangloc, nrem;
00282
00283 MPI_Comm_size(comm,&ncpus);
00284 MPI_Comm_rank(comm,&mypid);
00285
00286 nangloc = nang/ncpus;
00287 nrem = nang - ncpus*nangloc;
00288 if (mypid < nrem) nangloc++;
00289
00290 for (int i = 0; i < ncpus; i++) {
00291 psize[i] = nang/ncpus;
00292 if (i < nrem) psize[i]++;
00293 }
00294
00295 nbase[0] = 0;
00296 for (int i = 1; i < ncpus; i++) {
00297 nbase[i] = nbase[i-1] + psize[i-1];
00298 }
00299
00300 return nangloc;
00301 }
00302