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

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].
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.
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 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
4058
4059        // now we check to see if the radius is acceptable
4060        // That is, would not cross multiple ranks in MPI
4061
4063        {
4064            throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");
4065        }
4066
4067
4068
4069
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;
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;

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++));
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      {          {