EMAN2
Functions
utilcomm_Cart.h File Reference
#include "mpi.h"
#include "stdlib.h"
#include "emdata.h"
#include "alignoptions.h"
Include dependency graph for utilcomm_Cart.h:
This graph shows which files directly or indirectly include this file:

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, nraysloc, nx, ny, ptrs, status, x, and y.

Referenced by LBD_Cart(), main(), recons3d_CGLS_mpi_Cart(), and recons3d_sirt_mpi_Cart().

{
    int       i, j, iqx,iqy, xc, yc, zc, jj;
    float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
    int       status = 0; 

    int xcent = origin[0];
    int ycent = origin[1];
    int zcent = origin[2];

    int nx = volsize[0];
    int ny = volsize[1];

    // Phi: adding the shift parameters that get passed in as the last two entries of dm
    float sx, sy;

    sx = dm(7);
    sy = dm(8);

    if ( nx > 2*ri) {
        for (i = 1; i <= nraysloc; i++) {
            zc = cord(1,myptrstart+i) - zcent;
            yc = cord(2,myptrstart+i) - ycent;
            xc = cord(3,myptrstart+i) - xcent;

            xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent + sx;
            yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent + sy;

            for (j = ptrs(myptrstart+i); j <ptrs(myptrstart+i+1); j++) {
                jj = j-ptrs[myptrstart]+1;

                // jj is the index on the local volume
                iqx = ifix(xb);
                iqy = ifix(yb);

                dx = xb - iqx;
                dy = yb - iqy;
                dx1m = 1.0 - dx;
                dy1m = 1.0 - dy;
                dxdy = dx*dy;
/*
c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
c     &                     + dx1m*dy  *x(iqx  , iqy+1)
c     &                     + dx  *dy1m*x(iqx+1, iqy)
c     &                     + dx  *dy  *x(iqx+1, iqy+1)  
c
c              --- faster version of the above commented out
c                  code (derived by summing the following table 
c                  of coefficients along  the colunms) ---
c
c                        1         dx        dy      dxdy
c                     ------   --------  --------  -------
c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)  
c                                        x(i,j+1) -x(i,j+1)
c                              x(i+1,j)           -x(i+1,j)
c                                                x(i+1,j+1) 
c
*/
                // Phi: add index checking, now that shifts are being used
                if ( iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1 ) {
                    y(jj) += x(iqx,iqy);
                    if ( iqx + 1 <= nx && iqx + 1 >= 1 ) {
                        y(jj) += dx*(-x(iqx,iqy)+x(iqx+1,iqy));
                    }
                    if ( iqy + 1 <= ny && iqy + 1 >= 1 ) {
                        y(jj) += dy*(-x(iqx,iqy)+x(iqx,iqy+1));
                    }
                    if ( iqx + 1 <= nx && iqy + 1 <= ny && iqx + 1 >= 1 && iqy + 1 >= 1 ) {
                        y(jj) += dxdy*( x(iqx,iqy) - x(iqx,iqy+1) -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
                    }
                }

//                y(j) += x(iqx,iqy)
//                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
//                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
//                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1) 
//                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );

               xb += dm(1);
               yb += dm(4);
            } // end for j
        } // end for i
     }
    else {
        fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
    }

    return status;
}
int CleanStack_Cart ( MPI_Comm  comm,
EMData **  image_stack,
int  nloc,
int  ri,
Vec3i  volsize,
Vec3i  origin 
)
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, nraysloc, nx, ny, ptrs, status, x, and y.

Referenced by LBD_Cart(), main(), recons3d_CGLS_mpi_Cart(), and recons3d_sirt_mpi_Cart().

{
    int    iqx, iqy, i, j, xc, yc, zc, jj;
    float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
    int    status = 0;
    
    // Phi: adding the shift parameters that get passed in as the last two entries of dm
    float sx, sy;

    sx = dm(7);
    sy = dm(8);

    int xcent = origin[0];
    int ycent = origin[1];
    int zcent = origin[2];

    int nx = volsize[0];
    int ny = volsize[1];

    dm1 = dm(1);
    dm4 = dm(4);
 
    if ( nx > 2*ri ) {
        for (i = 1; i <= nraysloc; i++) {

            zc = cord(1,myptrstart+i)-zcent;
            yc = cord(2,myptrstart+i)-ycent;
            xc = cord(3,myptrstart+i)-xcent;
            xb = zc* dm(1) +yc* dm(2) +xc* dm(3) + xcent + sx;
            yb = zc* dm(4) +yc* dm(5) +xc* dm(6) + ycent + sy;

            for (j = ptrs(myptrstart+i); j< ptrs(myptrstart+i+1); j++) {
               jj = j-ptrs[myptrstart]+1;
               iqx = ifix(xb);
               iqy = ifix(yb);

                // jj is the index on the local volume
               ct   = x(jj);

               // dipx =  xb - (float)(iqx);
               // dipy = (yb - (float)(iqy)) * ct;
                   dipx =  xb - iqx;
                   dipy = (yb - iqy) * ct;

               dipy1m = ct - dipy;
               dipx1m = 1.0 - dipx;

                        if (iqx <= nx && iqy <= ny && iqx >= 1 && iqy >= 1) 
               // y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
               y(iqx  ,iqy)   +=  dipx1m*dipy1m;
                        if (iqx + 1 <= nx && iqy <= ny && iqx >= 0 && iqy >= 1) 
               // y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m; 
               y(iqx+1,iqy)   +=  dipx*dipy1m; 
                        if (iqx + 1 <= nx && iqy + 1 <= ny && iqx >= 0 && iqy >= 0) 
               // y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;         
               y(iqx+1,iqy+1) +=  dipx*dipy;         
                        if (iqx <= nx && iqy + 1 <= ny && iqx >= 1 && iqy >= 0) 
               // y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
               y(iqx  ,iqy+1) +=  dipx1m*dipy;
               xb += dm1;
               yb += dm4;
           }
        }
    }
    else {
        fprintf(stderr, " nx must be greater than 2*ri\n");
        exit(1);
    }
    return status;
}
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, and status.

Referenced by main(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), and recons3d_sirt_mpi_Cart().

{
    int    xs, ys, zs, xx, yy, zz, rs, r2;
    int    ix, iy, iz, jnz, nnz, nrays;
    int    ftm = 0, status = 0;  

    int xcent = (int)origin[0];
    int ycent = (int)origin[1];
    int zcent = (int)origin[2];

    int nx = (int)volsize[0];
    int ny = (int)volsize[1];
    int nz = (int)volsize[2];

    r2      = ri*ri;
    nnz     = 0;
    nrays    = 0;
    ptrs(1) = 1;

    for (ix = 1; ix <= nx; ix++) {
       xs  = ix-xcent;
       xx  = xs*xs;
       for ( iy = 1; iy <= ny; iy++ ) {
           ys = iy-ycent;
           yy = ys*ys;
           jnz = 0;

           ftm = 1;
           // not the most efficient implementation
           for (iz = 1; iz <= nz; iz++) {
               zs = iz-zcent;
               zz = zs*zs;
               rs = xx + yy + zz;
               if (rs <= r2) {
                  jnz++;
                  nnz++;
     //             sphere(nnz) = cube(iz, iy, ix); 

                  //  record the coordinates of the first nonzero ===
                  if (ftm) {
                     nrays++;
                     cord(1,nrays) = iz; 
                     cord(2,nrays) = iy; 
                     cord(3,nrays) = ix;
                     ftm = 0;
                  }
               }
            } // end for (iz..)
            if (jnz > 0) {
                ptrs(nrays+1) = ptrs(nrays) + jnz;
            }  // endif (jnz)
       } // end for iy
    } // end for ix
    if (nnz != nnz0) status = -1;
    return status;
}
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().

{
   int ierr = 0;
   int ROW = 0, COL = 1;
   int my2dpid, mycoords[2];
   int srpid, srcoords[2], keep_dims[2];
   MPI_Status mpistatus;
   FILE *fp = NULL;

   int nimgs=0;

   MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
   MPI_Cart_coords(comm_2d, my2dpid, 2, mycoords); // Get my coordinates

   float * iobuffer   = new float[5*nloc];
   if (!iobuffer) {
      fprintf(stderr,"failed to allocate buffer to read angles shifts\n");
      ierr = -1;
      goto EXIT;
   }

   if (mycoords[COL] == 0 && mycoords[ROW] == 0) { //I am Proc (0,0)
      fp = fopen(angfname,"r");
      if (!fp)  ierr = 1;
   }
   MPI_Bcast(&ierr, 1, MPI_INT, 0, comm_2d);

   if ( ierr ) {
      if (mycoords[COL] == 0 && mycoords[ROW] == 0) 
          fprintf(stderr,"failed to open %s\n", angfname);
      ierr = MPI_Finalize();
      goto EXIT;
   }
   else {
       if (mycoords[COL] == 0 && mycoords[ROW] == 0) { //I am Proc (0,0)
          for (int iproc = 0; iproc < dims[ROW]; iproc++) {
             // figure out the number of images assigned to processor (iproc,0)
             if (iproc > 0) {
                srcoords[COL] = 0;
                srcoords[ROW] = iproc;
                MPI_Cart_rank(comm_2d, srcoords, &srpid);

                MPI_Recv(&nimgs, 1, MPI_INT, srpid, srpid, comm_2d, &mpistatus);

                // Read the next nimgs set of angles and shifts
                for (int i = 0; i < nimgs; i++) {
                   fscanf(fp,"%f %f %f %f %f", 
                          &iobuffer[5*i+0],
                          &iobuffer[5*i+1],
                          &iobuffer[5*i+2],
                          &iobuffer[5*i+3],
                          &iobuffer[5*i+4]);
                }
                MPI_Send(iobuffer,5*nimgs,MPI_FLOAT,srpid,srpid,comm_2d);
             }
             else {
                for (int i = 0; i < nloc; i++) {
                   fscanf(fp,"%f %f %f %f %f", 
                          &angleshift[5*i+0],
                          &angleshift[5*i+1],
                          &angleshift[5*i+2],
                          &angleshift[5*i+3],
                          &angleshift[5*i+4]);
                }
             }
          }
          fclose(fp);
       }
       else if (mycoords[COL] == 0 && mycoords[ROW] != 0) { //I am in the first column
          // send image count to the master processor (mypid = 0)

          MPI_Send(&nloc, 1, MPI_INT, 0, my2dpid, comm_2d);
          // Receive angleshifts
          MPI_Recv(angleshift, 5*nloc, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
       }
  }

  // Now have all the processors in group g_c_0 broadcast the angles along the row communicator
  srcoords[ROW] = 0;
  MPI_Cart_rank(comm_row, srcoords, &srpid);
  MPI_Bcast(angleshift, 5*nloc, MPI_FLOAT, srpid, comm_row);

  EMDeleteArray(iobuffer);

EXIT:
   return ierr;
}
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_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().

{
    int ncpus, ierr, mpierr=0;
    MPI_Status mpistatus;
    EMUtil *my_util;
    int nima; 
    FILE *fp=NULL;
    int ROW = 0, COL = 1;
    int my2dpid, mycoords[2], dims[2], periods[2];
    int srpid, srcoords[2]; // Send/receive info

    // Get dims and my coordinates
    MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
    MPI_Comm_rank(comm_2d, &my2dpid); //Get my pid in the new 2D topology
    ncpus = dims[ROW]*dims[COL];

    ierr = 0;
    if (mycoords[ROW] == 0 && mycoords[COL] == 0) {
        fp = fopen(stackfname,"r");
        if (!fp) {
            ierr = 1;
            printf("failed to open %s\n", stackfname);
        }
        else {
            fclose(fp);
            nima = my_util->get_image_count(stackfname);
/**********************************************/
//    nima = 84;
/****************************/
        }
    }
    mpierr = MPI_Bcast(&ierr, 1, MPI_INT, 0, comm_2d);
    if (ierr == 0) {
        int *psize = new int[dims[ROW]];
        int *nbase = new int[dims[ROW]];
    
        mpierr = MPI_Bcast(&nima, 1, MPI_INT, 0, comm_2d);
    
        *nloc = setpart_gc1(comm_2d, nima, psize, nbase);
        *images2D = new EMData*[*nloc]; // NB!: whoever calls ReadStackandDist must delete this!

        EMData *img_ptr;
        int img_index;
        float *imgdata;

        // read the first image to get size
        img_ptr = new EMData();
        img_ptr->read_image(stackfname, 0);
        int nx = img_ptr->get_xsize();
        int ny = img_ptr->get_ysize();

        float s2x, s2y;

        if (mycoords[COL] == 0 && mycoords[ROW] == 0) {
            printf("Master node reading and distributing %d images along the first column of processors\n", nima);
            for ( int ip = 0 ; ip < dims[ROW] ; ++ip ) {
                for ( int i = 0 ; i < psize[ip] ; ++i ) {
                    img_index = nbase[ip] + i;

                    if (ip != 0) {
                        img_ptr->read_image(stackfname, img_index);
                        // get a pointer to the image's data
                        imgdata = img_ptr->get_data();
                        // find the x/y shift values if it has them, otherwise set them to 0.0
                        try {
                            s2x = (*images2D)[i]->get_attr("s2x");
                        } catch ( std::exception& e ) {
                            s2x = 0.0;
                        }
                        try {
                            s2y = (*images2D)[i]->get_attr("s2y");
                        } catch ( std::exception& e ) {
                            s2y = 0.0;
                        }
                        // send these to processor (ip,0)
                        srcoords[COL] = 0;
                        srcoords[ROW] = ip;
                        MPI_Cart_rank(comm_2d, srcoords, &srpid);

                        MPI_Send(imgdata, nx*ny, MPI_FLOAT, srpid, srpid, comm_2d);
                        MPI_Send(&s2x, 1, MPI_FLOAT, srpid, srpid, comm_2d);
                        MPI_Send(&s2y, 1, MPI_FLOAT, srpid, srpid, comm_2d);
                    } else { // ip == 0                             
                        (*images2D)[i] = new EMData();
                        (*images2D)[i]->read_image(stackfname, img_index);
                        try {
                            s2x = (*images2D)[i]->get_attr("s2x");
                        } catch ( std::exception& e ) {
                            (*images2D)[i]->set_attr("s2x",0.0);
                        }
                        try {
                            s2y = (*images2D)[i]->get_attr("s2y");
                        } catch ( std::exception& e ) {
                            (*images2D)[i]->set_attr("s2y",0.0);
                        }
                    }
                }
                //printf("finished reading data for processor %d\n", ip);
            }

        } else if (mycoords[COL] == 0 && mycoords[ROW] != 0) {
            for ( int i = 0 ; i < psize[mycoords[ROW]] ; ++i ) {
                (*images2D)[i] = new EMData();
                (*images2D)[i]->set_size(nx, ny, 1);
                imgdata = (*images2D)[i]->get_data();
                
                MPI_Recv(imgdata, nx*ny, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
                MPI_Recv(&s2x, 1, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
                MPI_Recv(&s2y, 1, MPI_FLOAT, 0, my2dpid, comm_2d, &mpistatus);
                (*images2D)[i]->set_attr("s2x",s2x);
                (*images2D)[i]->set_attr("s2y",s2y);
            }
            //printf("received %d images for processor (%d,%d)\n", *nloc, mycoords[ROW], mycoords[COL]);
        }

// Now have all the processors in group g_c_0 broadcast the images and angles along the row communicator
    if (mycoords[COL] == 0 && mycoords[ROW] == 0)
      printf("First column of processors distributing images to other columns..\n");

    if (mycoords[COL] == 0) {
      for (int img_index = 0; img_index < *nloc; img_index++){
        imgdata = (*images2D)[img_index]->get_data();
        
        // find the x/y shift values if it has them, otherwise set them to 0.0
        try {
            s2x = (*images2D)[img_index]->get_attr("s2x");
        } catch ( std::exception& e ) {
            s2x = 0.0;
        }
        try {
            s2y = (*images2D)[img_index]->get_attr("s2y");
        } catch ( std::exception& e ) {
            s2y = 0.0;
        }

        for(int ip=1; ip< dims[COL]; ip++){
// printf("Proc (%d, %d) sending image %d of %d to Proc (%d,%d) \n", mycoords[ROW], mycoords[COL], img_index, *nloc, mycoords[ROW], ip); 
          srcoords[ROW] = mycoords[ROW];
          srcoords[COL] = ip;
          MPI_Cart_rank(comm_2d, srcoords, &srpid);

          MPI_Send(imgdata, nx*ny, MPI_FLOAT, srpid, img_index, comm_2d);
          MPI_Send(&s2x, 1, MPI_FLOAT, srpid, img_index, comm_2d);
          MPI_Send(&s2y, 1, MPI_FLOAT, srpid, img_index, comm_2d);
        }
      }
    }
    if (mycoords[COL] != 0) {
  //      printf("Proc (%d, %d) receiving...", mycoords[ROW], mycoords[COL]);
        for(int img_index = 0; img_index < *nloc; img_index++){
          (*images2D)[img_index] = new EMData();
          (*images2D)[img_index]->set_size(nx, ny, 1);
          imgdata = (*images2D)[img_index]->get_data();
          
          srcoords[ROW] = mycoords[ROW];
          srcoords[COL] = 0;
          MPI_Cart_rank(comm_2d, srcoords, &srpid);

          MPI_Recv(imgdata, nx*ny, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);
          MPI_Recv(&s2x, 1, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);
          MPI_Recv(&s2y, 1, MPI_FLOAT, srpid, img_index, comm_2d, &mpistatus);

          (*images2D)[img_index]->set_attr("s2x",s2x);
          (*images2D)[img_index]->set_attr("s2y",s2y);
         }
    }

        ierr = mpierr; 

        EMDeletePtr(img_ptr);
        EMDeleteArray(psize);
        EMDeleteArray(nbase);
    }
    return ierr;
}
int setpart_gc1 ( MPI_Comm  comm_2d,
int  nangs,
int *  psize,
int *  nbase 
)

Definition at line 323 of file utilcomm_Cart.cpp.

{
        int ROW = 0, COL = 1;
        int dims[2], periods[2], mycoords[2];
        int nangsloc, nrem;

        // Get information associated with comm_2d
        MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
 
        nangsloc = nangs/dims[ROW];
        nrem = nangs - dims[ROW]*nangsloc;
        if (mycoords[ROW] < nrem) nangsloc++;

        for (int i = 0; i < dims[ROW]; i++) {
                psize[i] = nangs/dims[ROW];
                if (i < nrem) psize[i]++;
        }
 
        for (int i = 0; i < dims[ROW]; i++) {
          if(i==0)
            nbase[i] = 0;
          else
            nbase[i] = nbase[i-1] + psize[i-1];
        }
   
//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]);

        return nangsloc;
}
int setpart_gr1 ( MPI_Comm  comm_2d,
int  nnz,
int *  nnzpart,
int *  nnzbase 
)

Definition at line 364 of file utilcomm_Cart.cpp.

{
        int ROW = 0, COL = 1;
        int dims[2], periods[2], mycoords[2];
        int nnzloc, nrem;
        
        // Get information associated with comm_2d
        MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);
 
        nnzloc = nnz/dims[COL];
        nrem = nnz - dims[COL]*nnzloc;
        if (mycoords[COL] < nrem) nnzloc++;

        for (int i = 0; i < dims[COL]; i++) {
                nnzpart[i] = nnz/dims[COL];
                if (i < nrem) nnzpart[i]++;
        }
 
        nnzbase[0] = 0; 
        for (int i = 1; i < dims[COL]; i++) {
                nnzbase[i] = nnzbase[i-1] + nnzpart[i-1];
        }
        nnzbase[dims[COL]] = nnz;
        return nnzloc;
}
int sphpart ( MPI_Comm  comm_2d,
int  nrays,
int *  ptrs,
int *  nnzbase,
int *  ptrstart 
)

Definition at line 85 of file project3d_Cart.cpp.

{
        int ROW = 0, COL = 1;
        int dims[2], periods[2], mycoords[2];
        int nraysloc = 0;
        
        // Get information associated with comm_2d
        MPI_Cart_get(comm_2d, 2, dims, periods, mycoords);

        int count = 1;
        if (mycoords[COL] == 0){ // First column starts out with the first ray
          nraysloc++;
        }
        ptrstart[0] = 0;

        for(int i=1; i<nrays ; i++){
          if (ptrs[i]-1 <= nnzbase[count] && ptrs[i+1]-1 >= nnzbase[count]){ 
                //nnzbase is between or equal to ptrs

            if (nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]){
              if(mycoords[COL] == count-1){  // ray belongs to count-1
                nraysloc++;
              }
              ptrstart[count] = i+1;
              count++;

            } else { //nnzbase[count]-(ptrs[i]-1)>= ptrs[i+1]-1-nnzbase[count]
                if(mycoords[COL] == count){// ray belongs to count and it's a new starting ray
                  nraysloc++;
                }
                ptrstart[count] = i;
                count++;
            }

        }
        else {  //ptrs[i]-1 > nnzbase[count] so ray belongs to count-1
          if(mycoords[COL] == count-1){
            nraysloc++;
          }

        }
        } // end for
        ptrstart[dims[COL]] = nrays;
        return nraysloc;

}