/[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 4529 by caltinay, Fri Oct 25 01:23:27 2013 UTC revision 4572 by jfenwick, Sun Dec 8 00:54:37 2013 UTC
# Line 3990  escript::Data Rectangle::randomFill(long Line 3990  escript::Data Rectangle::randomFill(long
3990          throw RipleyException("Radius of gaussian filter must be a positive integer.");          throw RipleyException("Radius of gaussian filter must be a positive integer.");
3991      }      }
3992      unsigned int radius=ex1();      unsigned int radius=ex1();
3993    #ifdef ESYS_MPI    
3994    
3995        // Need to check to see that radius would not cause the overlap to cover
3996        // more than one cell (eg each rank holds 4 columns and the radius is 5).
3997        // Also need to take special care with narrow cells
3998        
3999        
4000        // In fact it needs to be stricter than this, if a rank has neighbours on both sides, the borders can't overlap.
4001        
4002    #endif    
4003      double sigma=0.5;      double sigma=0.5;
4004      boost::python::extract<double> ex2(filter[2]);      boost::python::extract<double> ex2(filter[2]);
4005      if (!ex2.check() || (sigma=ex2())<=0)      if (!ex2.check() || (sigma=ex2())<=0)
# Line 4008  escript::Data Rectangle::randomFill(long Line 4018  escript::Data Rectangle::randomFill(long
4018      esysUtils::randomFillArray(seed, src, dsize);        esysUtils::randomFillArray(seed, src, dsize);  
4019                
4020      // 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  
4021      #ifdef ESYS_MPI    
4022        
4023        
4024        dim_t X=m_mpiInfo->rank%m_NX[0];
4025        dim_t Y=m_mpiInfo->rank/m_NX[0];
4026    
4027        MPI_Request reqs[10];
4028        MPI_Status stats[10];
4029        short rused=0;
4030        double* SWin=new double[radius*radius];
4031        double* SEin=new double[radius*radius];
4032        double* NWin=new double[radius*radius];
4033        double* Sin=new double[radius*m_NX[0]];
4034        double* Win=new double[radius*m_NX[0]];
4035    
4036        double* NEout=new double[radius*radius];
4037        double* NWout=new double[radius*radius];
4038        double* SEout=new double[radius*radius];
4039        double* Nout=new double[radius*m_NX[0]];
4040        double* Eout=new double[radius*m_NX[0]];
4041        
4042        int comserr=0;
4043        if (Y!=0)   // not on bottom row,
4044        {
4045        if (X!=0)   // not on the left hand edge
4046        {
4047            // recv bottomleft from SW
4048            comserr|=MPI_Irecv(SWin, radius*radius, MPI_DOUBLE, (X-1)+(Y-1)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4049            comserr|=MPI_Irecv(Win, numpoints[1]*radius, MPI_DOUBLE, X-1+Y*m_NX[0], 10, m_mpiInfo->comm, reqs+(rused++));    
4050        }
4051        else
4052        {
4053            comserr|=MPI_Irecv(SWin, radius*radius, MPI_DOUBLE, (Y-1)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4054        }
4055        comserr|=MPI_Irecv(Sin, numpoints[0]*radius, MPI_DOUBLE, X+(Y-1)*m_NX[0], 8, m_mpiInfo->comm, reqs+(rused++));
4056        comserr|=MPI_Irecv(SEin, radius*radius, MPI_DOUBLE, X+(Y-1)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4057          
4058        }
4059        else        // on the bottom row
4060        {
4061        if (X!=0)
4062        {
4063            comserr|=MPI_Irecv(Win, numpoints[1]*radius, MPI_DOUBLE, X-1+Y*m_NX[0], 10, m_mpiInfo->comm, reqs+(rused++));
4064            comserr|=MPI_Irecv(NWin, radius*radius, MPI_DOUBLE, X-1+Y*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4065        }
4066        if (X!=(m_NX[0]-1))
4067        {
4068            comserr|=MPI_Isend(SEout, radius*radius, MPI_DOUBLE, X+1+(Y)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));    
4069        }
4070        }
4071        if (Y!=(m_NX[1]-1)) // not on the top row
4072        {
4073        comserr|=MPI_Isend(Nout, radius*numpoints[0], MPI_DOUBLE, X+(Y+1)*m_NX[0], 8, m_mpiInfo->comm, reqs+(rused++));
4074        comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+(Y+1)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4075        if (X!=(m_NX[0]-1)) // not on right hand edge
4076        {
4077            comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+1+(Y+1)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4078        }
4079        if (X==0)   // left hand edge
4080        {
4081            comserr|=MPI_Isend(NWout, radius*radius, MPI_DOUBLE, (Y+1)*m_NX[0],7, m_mpiInfo->comm, reqs+(rused++));
4082        }  
4083        }
4084        if (X!=(m_NX[0]-1)) // not on right hand edge
4085        {
4086        comserr|=MPI_Isend(NEout, radius*radius, MPI_DOUBLE, X+1+(Y)*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4087        comserr|=MPI_Isend(Eout, numpoints[1]*radius, MPI_DOUBLE, X+1+(Y)*m_NX[0], 10, m_mpiInfo->comm, reqs+(rused++));
4088        comserr|=MPI_Irecv(NWin, radius*radius, MPI_DOUBLE, (X-1)+Y*m_NX[0], 7, m_mpiInfo->comm, reqs+(rused++));
4089        }
4090    
4091        if (!comserr)
4092        {
4093            comserr=MPI_Waitall(rused, reqs, stats);
4094        }
4095        
4096        if (comserr)
4097        {
4098        // Yes this is throwing an exception as a result of an MPI error.
4099        // and no we don't inform the other ranks that we are doing this.
4100        // however, we have no reason to believe coms work at this point anyway
4101            throw RipleyException("Error in coms for randomFill");      
4102        }
4103        
4104        
4105        delete[] SWin;
4106        delete[] SEin;
4107        delete[] NWin;
4108        delete[] Sin;
4109        delete[] Win;
4110    
4111        delete[] NEout;
4112        delete[] NWout;
4113        delete[] SEout;
4114        delete[] Nout;
4115        delete[] Eout;
4116        
4117        
4118        
4119        
4120    #endif    
4121      // Lets call that done for now      // Lets call that done for now
4122      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
4123      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);

Legend:
Removed from v.4529  
changed lines
  Added in v.4572

  ViewVC Help
Powered by ViewVC 1.1.26