/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4575 by jfenwick, Mon Dec 9 22:59:31 2013 UTC revision 4581 by jfenwick, Tue Dec 10 09:32:50 2013 UTC
# Line 3966  namespace Line 3966  namespace
3966    
3967    
3968    
3969    /* This routine produces a Data object filled with smoothed random data.
3970    The dimensions of the rectangle being filled are internal[0] x internal[1] points.
3971    A parameter radius  dives the size of the stencil used for the smoothing.
3972    A point on the left hand edge for example, will still require `radius` extra points to the left
3973    in order to complete the stencil.
3974    
3975    All local calculation is done on an array called `src`, with
3976    dimensions = ext[0] * ext[1].
3977    Where ext[i]= internal[i]+2*radius.
3978    
3979    Now for MPI there is overlap to deal with. We need to share both the overlapping
3980    values themselves but also the external region.
3981    
3982    In a hypothetical 1-D case:
3983    
3984    
3985    1234567
3986    would be split into two ranks thus:
3987    123(4)  (4)567     [4 being a shared element]
3988    
3989    If the radius is 2.   There will be padding elements on the outside:
3990    
3991    pp123(4)  (4)567pp
3992    
3993    To ensure that 4 can be correctly computed on both ranks, values from the other rank
3994    need to be known.
3995    
3996    pp123(4)56   23(4)567pp
3997    
3998    Now in our case, we wout set all the values 23456 on the left rank and send them to the
3999    right hand rank.
4000    
4001    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
4002    inset=2*radius+1    
4003    This is to ensure that values at distance `radius` from the shared/overlapped element
4004    that ripley has.
4005    
4006    
4007    */
4008  escript::Data Rectangle::randomFill(long seed, const boost::python::tuple& filter) const  escript::Data Rectangle::randomFill(long seed, const boost::python::tuple& filter) const
4009  {  {
4010  //     if (m_mpiInfo->size!=1)  //     if (m_mpiInfo->size!=1)
# Line 3999  escript::Data Rectangle::randomFill(long Line 4038  escript::Data Rectangle::randomFill(long
4038            
4039      // In fact it needs to be stricter than this, if a rank has neighbours on both sides, the borders can't overlap.      // In fact it needs to be stricter than this, if a rank has neighbours on both sides, the borders can't overlap.
4040            
4041        // basically, 2*inset<ext[i]
4042        
4043        
4044  #endif      #endif    
4045      double sigma=0.5;      double sigma=0.5;
4046      boost::python::extract<double> ex2(filter[2]);      boost::python::extract<double> ex2(filter[2]);
# Line 4006  escript::Data Rectangle::randomFill(long Line 4048  escript::Data Rectangle::randomFill(long
4048      {      {
4049          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Sigma must be a postive floating point number.");
4050      }          }    
     size_t numpoints[2];  
     numpoints[0]=m_ownNE[0]+1;  
     numpoints[1]=m_ownNE[1]+1;  
     size_t padding=max((unsigned)max((m_NE[0]-m_ownNE[0])/2, (m_NE[1]-m_ownNE[1])/2), radius);  
     size_t width=(numpoints[0]+2*padding);      // width of one row in points  
     size_t height=(numpoints[1]+2*padding); // height of one row in points  
     size_t dsize=width * height; // size of padded source grid  
4051            
4052      double* src=new double[dsize];      size_t internal[2];
4053      esysUtils::randomFillArray(seed, src, dsize);        internal[0]=m_NE[0]+1;  // number of points in the internal region
4054        internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of
4055        size_t ext[2];
4056        ext[0]=internal[0]+2*radius;    // includes points we need as input
4057        ext[1]=internal[1]+2*radius;    // for smoothing
4058        
4059        // now we check to see if the radius is acceptable
4060        // That is, would not cross multiple ranks in MPI
4061    
4062        if ((2*radius>=internal[0]) || (2*radius>=internal[1]))
4063        {
4064            throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");
4065        }
4066        
4067        
4068        
4069        
4070        size_t inset=2*radius+1;
4071        size_t Eheight=ext[1]-2*inset;  // how high the E (shared) region is
4072        size_t Swidth=ext[0]-2*inset;
4073    
4074        
4075    //     size_t numpoints[2];
4076    //     numpoints[0]=m_ownNE[0]+1;
4077    //     numpoints[1]=m_ownNE[1]+1;
4078    //     size_t padding=max((unsigned)max((m_NE[0]-m_ownNE[0])/2, (m_NE[1]-m_ownNE[1])/2), radius);
4079    //     size_t width=(numpoints[0]+2*padding);   // width of one row in points
4080    //     size_t height=(numpoints[1]+2*padding);  // height of one row in points
4081    //     size_t dsize=width * height; // size of padded source grid
4082        
4083        double* src=new double[ext[0]*ext[1]];
4084        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]);  
4085                
4086      // Now we need to copy the regions owned by other ranks over here        // Now we need to copy the regions owned by other ranks over here  
4087  #ifdef ESYS_MPI      //#ifdef ESYS_MPI    
4088            
4089    
4090    
4091        double* SWin=new double[inset*inset];  memset(SWin, 0, inset*inset*sizeof(double));
4092        double* SEin=new double[inset*inset];  memset(SEin, 0, inset*inset*sizeof(double));
4093        double* NWin=new double[inset*inset];  memset(NWin, 0, inset*inset*sizeof(double));
4094        double* Sin=new double[inset*Swidth];  memset(Sin, 0, inset*Swidth*sizeof(double));
4095        double* Win=new double[inset*Eheight];  memset(Win, 0, inset*Eheight*sizeof(double));
4096    
4097        double* NEout=new double[inset*inset];  memset(NEout, 0, inset*inset*sizeof(double));
4098        uint base=ext[0]-inset+(ext[1]-inset)*ext[0];
4099        for (uint i=0;i<inset;++i)
4100        {
4101        memcpy(NEout+inset*i, src+base, inset*sizeof(double));
4102        base+=ext[0];
4103        }
4104        double* NWout=new double[inset*inset];  memset(NWout, 0, inset*inset*sizeof(double));
4105        base=(ext[1]-inset)*ext[0];
4106        for (uint i=0;i<inset;++i)
4107        {
4108        memcpy(NWout+inset*i, src+base, inset*sizeof(double));
4109        base+=ext[0];
4110        }
4111        double* SEout=new double[inset*inset];  memset(SEout, 0, inset*inset*sizeof(double));
4112        base=ext[0]-inset;
4113        for (int i=0;i<inset;++i)
4114        {
4115        memcpy(SEout+inset*i, src+base, inset*sizeof(double));
4116        base+=ext[0];
4117        }
4118        double* Nout=new double[inset*Swidth];  memset(Nout, 0, inset*Swidth*sizeof(double));
4119        base=inset+(ext[1]-inset)*ext[0];
4120        for (uint i=0;i<inset;++i)
4121        {
4122        memcpy(Nout, src+base, Swidth*sizeof(double));
4123        base+=ext[0];
4124        }
4125        double* Eout=new double[inset*Eheight];  memset(Eout, 0, inset*Eheight*sizeof(double));
4126        base=ext[0]-inset+inset*ext[0];
4127        for (uint i=0;i<inset;++i)
4128        {
4129        memcpy(Eout, src+base, inset);
4130        base+=ext[0];
4131        }
4132            
4133      dim_t X=m_mpiInfo->rank%m_NX[0];  #ifdef ESYS_MPI    
     dim_t Y=m_mpiInfo->rank/m_NX[0];  
     dim_t row=m_NX[0];  
4134    
4135      MPI_Request reqs[10];      MPI_Request reqs[10];
4136      MPI_Status stats[10];      MPI_Status stats[10];
4137      short rused=0;      short rused=0;
     double* SWin=new double[radius*radius];  memset(SWin, 0, radius*radius*sizeof(double));  
     double* SEin=new double[radius*radius];  memset(SEin, 0, radius*radius*sizeof(double));  
     double* NWin=new double[radius*radius];  memset(NWin, 0, radius*radius*sizeof(double));  
     double* Sin=new double[radius*numpoints[0]];  memset(Sin, 0, radius*numpoints[0]*sizeof(double));  
     double* Win=new double[radius*numpoints[1]];  memset(Win, 0, radius*numpoints[1]*sizeof(double));  
   
     double* NEout=new double[radius*radius];  memset(NEout, 0, radius*radius*sizeof(double));  
     double* NWout=new double[radius*radius];  memset(NWout, 0, radius*radius*sizeof(double));  
     double* SEout=new double[radius*radius];  memset(SEout, 0, radius*radius*sizeof(double));  
     double* Nout=new double[radius*numpoints[0]];  memset(Nout, 0, radius*numpoints[0]*sizeof(double));  
     double* Eout=new double[radius*numpoints[1]];  memset(Eout, 0, radius*numpoints[1]*sizeof(double));  
4138            
4139        dim_t X=m_mpiInfo->rank%m_NX[0];
4140        dim_t Y=m_mpiInfo->rank/m_NX[0];
4141        dim_t row=m_NX[0];
4142        bool swused=false;      // These vars will be true if data needs to be copied out of them
4143        bool seused=false;
4144        bool nwused=false;
4145        bool sused=false;
4146        bool wused=false;    
4147        
4148        
4149        
4150    
4151  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << SWin << " " << SEin << " " << NWin << " "<< Sin << " "<< Win << " "<< NEout << " "<< NWout << " "  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << SWin << " " << SEin << " " << NWin << " "<< Sin << " "<< Win << " "<< NEout << " "<< NWout << " "
4152  //<< SEout << " "  //<< SEout << " "
4153  //<< Nout << " "  //<< Nout << " "
# Line 4053  escript::Data Rectangle::randomFill(long Line 4161  escript::Data Rectangle::randomFill(long
4161      if (X!=0)   // not on the left hand edge      if (X!=0)   // not on the left hand edge
4162      {      {
4163          // recv bottomleft from SW          // recv bottomleft from SW
4164  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SW (7) from " << (X-1)+(Y-1)*row << endl;  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SW (7) from " << (X-1)+(Y-1)*row << endl;
4165          comserr|=MPI_Irecv(SWin, radius*radius, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
4166  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv W (10) from " <<  X-1+Y*row << endl;                  swused=true;
4167          comserr|=MPI_Irecv(Win, numpoints[1]*radius, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv W (10) from " <<  X-1+Y*row << endl;      
4168                comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));
4169            wused=true;
4170      }      }
4171      else      else
4172      {      {
4173  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SW (7) from " << (Y-1)*row  << endl;                cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SW (7) from " << (Y-1)*row  << endl;            
4174          comserr|=MPI_Irecv(SWin, radius*radius, MPI_DOUBLE, (Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
4175            swused=true;
4176      }      }
4177  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv S (8) from " << X+(Y-1)*row  << endl;          cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv S (8) from " << X+(Y-1)*row  << endl;          
4178      comserr|=MPI_Irecv(Sin, numpoints[0]*radius, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));      comserr|=MPI_Irecv(Sin, Swidth*inset, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));
4179  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SE (7) from " << X+(Y-1)*row << endl;          sused=true;
4180      comserr|=MPI_Irecv(SEin, radius*radius, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv SE (7) from " << X+(Y-1)*row << endl;      
4181        comserr|=MPI_Irecv(SEin, inset*inset, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
4182        seused=true;
4183    
4184                
4185      }      }
# Line 4075  escript::Data Rectangle::randomFill(long Line 4187  escript::Data Rectangle::randomFill(long
4187      {      {
4188      if (X!=0)      if (X!=0)
4189      {      {
4190  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv W (10) from " << X-1+Y*row << endl;      cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv W (10) from " << X-1+Y*row << endl;      
4191          comserr|=MPI_Irecv(Win, numpoints[1]*radius, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));
4192  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv NW (7) from " << X-1+Y*row << endl;                wused=true;
4193          comserr|=MPI_Irecv(NWin, radius*radius, MPI_DOUBLE, X-1+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv NW (7) from " << X-1+Y*row << endl;        
4194            comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, X-1+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));
4195            nwused=true;
4196      }      }
4197      if (X!=(row-1))      if (X!=(row-1))
4198      {      {
4199  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send SE (7) to " <<  X+1+(Y)*row << endl;            cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send SE (7) to " <<  X+1+(Y)*row << endl;        
4200          comserr|=MPI_Isend(SEout, radius*radius, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));            comserr|=MPI_Isend(SEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));    
4201            
4202      }      }
4203      }      }
4204            
4205      if (Y!=(m_NX[1]-1)) // not on the top row      if (Y!=(m_NX[1]-1)) // not on the top row
4206      {      {
4207  //cerr << "Y=" << Y << "  (numpoints[1]-1)=" << (numpoints[1]-1) << endl;  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send N (8) to " << X+(Y+1)*row  << endl;
4208  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send N (8) to " << X+(Y+1)*row  << endl;      comserr|=MPI_Isend(Nout, inset*Swidth, MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));
4209      comserr|=MPI_Isend(Nout, radius*numpoints[0], MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));  cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " <<  X+(Y+1)*row << endl;
4210  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " <<  X+(Y+1)*row << endl;      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
     comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
4211      if (X!=(row-1)) // not on right hand edge      if (X!=(row-1)) // not on right hand edge
4212      {      {
4213  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " <<  X+1+(Y+1)*row << endl;              cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " <<  X+1+(Y+1)*row << endl;              
4214          comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));
4215      }      }
4216      if (X==0)   // left hand edge      if (X==0)   // left hand edge
4217      {      {
4218  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NW (7) to " << (Y+1)*row << endl;        cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NW (7) to " << (Y+1)*row << endl;    
4219          comserr|=MPI_Isend(NWout, radius*radius, MPI_DOUBLE, (Y+1)*row,7, m_mpiInfo->comm, reqs+(rused++));              comserr|=MPI_Isend(NWout, inset*inset, MPI_DOUBLE, (Y+1)*row,7, m_mpiInfo->comm, reqs+(rused++));      
4220      }        }  
4221      }      }
4222      if (X!=(row-1)) // not on right hand edge      if (X!=(row-1)) // not on right hand edge
4223      {      {
4224  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " << X+1+(Y)*row << endl;        cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send NE (7) to " << X+1+(Y)*row << endl;      
4225      comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));
4226  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send E (10) to " << X+1+(Y)*row << endl;    cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Send E (10) to " << X+1+(Y)*row << endl;    
4227      comserr|=MPI_Isend(Eout, numpoints[1]*radius, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));      comserr|=MPI_Isend(Eout, Eheight*inset, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));
4228      }      }
4229      if (X!=0)      if (X!=0)
4230      {      {
4231  //cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv NW (7) from " << (X-1)+Y*row << endl;      cerr << m_mpiInfo->rank << "[" << __LINE__ << "]: " << "Recv NW (7) from " << (X-1)+Y*row << endl;      
4232                
4233      comserr|=MPI_Irecv(NWin, radius*radius, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));      comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));
4234        nwused=true;
4235                
4236                
4237      }      }
# Line 4136  escript::Data Rectangle::randomFill(long Line 4250  escript::Data Rectangle::randomFill(long
4250          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
4251      }      }
4252            
4253    
4254        // Need to copy the values back from the buffers
4255        // Copy from SW
4256            
4257        if (swused)
4258        {
4259        base=0;
4260        for (uint i=0;i<inset;++i)
4261        {
4262            memcpy(src+base, SWin+i*inset, inset*sizeof(double));
4263            base+=ext[0];
4264        }
4265        }
4266        if (seused)
4267        {
4268            base=ext[0]-inset;
4269        for (uint i=0;i<inset;++i)
4270        {
4271            memcpy(src+base, SEin+i*inset, inset*sizeof(double));
4272            base+=ext[0];
4273        }
4274        }
4275        if (nwused)
4276        {
4277            base=(ext[1]-inset)*ext[0];
4278        for (uint i=0;i<inset;++i)
4279        {
4280            memcpy(src+base, NWin+i*inset, inset*sizeof(double));
4281            base+=ext[0];
4282        }
4283        }
4284        if (sused)
4285        {
4286           base=inset;
4287           for (uint i=0;i<inset;++i)
4288           {
4289           memcpy(src+base, Sin+i*Swidth, Swidth*sizeof(double));
4290           base+=ext[0];
4291           }
4292        }
4293        if (wused)
4294        {
4295        base=inset*ext[0];
4296        for (uint i=0;i<Eheight;++i)
4297        {
4298            memcpy(src+base, Win+i*inset, inset*sizeof(double));
4299            base+=ext[0];
4300        }
4301          
4302        }
4303      delete[] SWin;      delete[] SWin;
4304      delete[] SEin;      delete[] SEin;
4305      delete[] NWin;      delete[] NWin;
# Line 4163  escript::Data Rectangle::randomFill(long Line 4326  escript::Data Rectangle::randomFill(long
4326      {      {
4327          for (size_t x=0;x<(m_ownNE[0]+1);++x)          for (size_t x=0;x<(m_ownNE[0]+1);++x)
4328      {          {    
4329          dv[x+y*(m_ownNE[0]+1)]=Convolve2D(convolution, src, x+radius, y+radius, radius, width);          dv[x+y*(m_ownNE[0]+1)]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);
4330      }      }
4331      }      }
4332      delete[] convolution;      delete[] convolution;

Legend:
Removed from v.4575  
changed lines
  Added in v.4581

  ViewVC Help
Powered by ViewVC 1.1.26