/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

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

revision 2673 by jfenwick, Fri Sep 18 05:33:10 2009 UTC revision 2677 by jfenwick, Tue Sep 22 00:48:00 2009 UTC
# Line 3046  Data Line 3046  Data
3046  Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3047                   double undef, bool check_boundaries)                   double undef, bool check_boundaries)
3048  {  {
3049        table.convertArray();       // critical!  Calling getElt on an unconverted array is not thread safe
3050      int error=0;      int error=0;
3051      if ((getDataPointRank()!=0))      if ((getDataPointRank()!=0))
3052      {      {
# Line 3064  Data::interpolateFromTable1D(const Wrapp Line 3065  Data::interpolateFromTable1D(const Wrapp
3065          expand();          expand();
3066      }      }
3067      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3068      do                                   // to make breaks useful      try
3069      {      {
3070          try          int numpts=getNumDataPoints();
3071            const DataVector& adat=getReady()->getVectorRO();
3072            DataVector& rdat=res.getReady()->getVectorRW();
3073            int twidth=table.getShape()[0]-1;
3074            bool haserror=false;
3075            int l=0;
3076            #pragma omp parallel for private(l) schedule(static)
3077            for (l=0;l<numpts; ++l)
3078          {          {
3079              int numpts=getNumDataPoints();            #pragma omp flush(haserror)       // In case haserror was in register
3080              const DataVector& adat=getReady()->getVectorRO();            if (!haserror)        
3081              DataVector& rdat=res.getReady()->getVectorRW();            {
3082              int twidth=table.getShape()[0]-1;             int lerror=0;
3083              for (int l=0; l<numpts; ++l)             try
3084               {
3085                 do         // so we can use break
3086                 {
3087                double a=adat[l];
3088                int x=static_cast<int>(((a-Amin)/Astep));
3089                if (check_boundaries) {
3090                    if ((a<Amin) || (x<0))
3091                    {
3092                        lerror=1;
3093                        break;  
3094                    }
3095                    if (a>Amin+Astep*twidth)
3096                    {
3097                            lerror=4;
3098                        break;
3099                    }
3100                }
3101                if (x<0) x=0;
3102                if (x>twidth) x=twidth;
3103    
3104                if (x==twidth)          // value is on the far end of the table
3105              {              {
3106                  double a=adat[l];                  double e=table.getElt(x);
3107                  int x=static_cast<int>(((a-Amin)/Astep));                  if (e>undef)
3108                                  if (check_boundaries) {                  {
3109                      if ((a<Amin) || (x<0))                      lerror=2;
3110                      {                      break;
                         error=1;  
                         break;    
                     }  
                     if (a>Amin+Astep*twidth)  
                     {  
                         error=4;  
                         break;  
                     }  
                                 }  
                                 if (x<0) x=0;  
                                 if (x>twidth) x=twidth;  
   
                 if (x==twidth)          // value is on the far end of the table  
                 {  
                     double e=table.getElt(x);  
                     if (e>undef)  
                     {  
                         error=2;  
                         break;  
                     }  
                     rdat[l]=e;  
                 }  
                 else        // x and y are in bounds  
                 {  
                     double e=table.getElt(x);  
                     double w=table.getElt(x+1);  
                     if ((e>undef) || (w>undef))  
                     {  
                         error=2;  
                         break;  
                     }  
             // map x*Astep <= a << (x+1)*Astep to [-1,1]  
                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;  
                     rdat[l]=((1-la)*e + (1+la)*w)/2;  
3111                  }                  }
3112                    rdat[l]=e;
3113              }              }
3114          } catch (DataException d)              else        // x and y are in bounds
3115          {              {
3116              error=3;                  double e=table.getElt(x);
3117              break;                  double w=table.getElt(x+1);
3118          }                  if ((e>undef) || (w>undef))
3119      } while (false);                  {
3120                        lerror=2;
3121                        break;
3122                    }
3123            // map x*Astep <= a << (x+1)*Astep to [-1,1]
3124                    double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3125                    rdat[l]=((1-la)*e + (1+la)*w)/2;
3126                }  
3127                  } while (false);
3128                } catch (DataException d)
3129                {
3130                    lerror=3;
3131                }
3132                if (lerror!=0)
3133                {
3134                #pragma omp critical    // Doco says there is a flush associated with critical
3135                {
3136                    haserror=true;  // We only care that one error is recorded. We don't care which
3137                    error=lerror;   // one
3138                }
3139                }
3140              } // if (!error)
3141            }   // parallelised for
3142        } catch (DataException d)
3143        {
3144            error=3;        // this is outside the parallel region so assign directly
3145        }
3146  #ifdef PASO_MPI  #ifdef PASO_MPI
3147      int rerror=0;      int rerror=0;
3148      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
# Line 3145  Data Line 3167  Data
3167  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3168                         double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)                         double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
3169  {  {
3170        table.convertArray();       // critical!   Calling getElt on an unconverted array is not thread safe
3171      int error=0;      int error=0;
3172      if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))      if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3173      {      {
# Line 3173  Data::interpolateFromTable2D(const Wrapp Line 3196  Data::interpolateFromTable2D(const Wrapp
3196      B.expand();      B.expand();
3197      }      }
3198      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3199      do                                   // to make breaks useful      try
3200      {      {
3201      try      int numpts=getNumDataPoints();
3202        const DataVector& adat=getReady()->getVectorRO();
3203        const DataVector& bdat=B.getReady()->getVectorRO();
3204        DataVector& rdat=res.getReady()->getVectorRW();
3205        const DataTypes::ShapeType& ts=table.getShape();
3206        int twx=ts[0]-1;    // table width x
3207        int twy=ts[1]-1;    // table width y
3208        bool haserror=false;
3209        int l=0;
3210        #pragma omp parallel for private(l) schedule(static)
3211        for (l=0; l<numpts; ++l)
3212      {      {
3213          int numpts=getNumDataPoints();         #pragma omp flush(haserror)      // In case haserror was in register
3214          const DataVector& adat=getReady()->getVectorRO();         if (!haserror)      
3215          const DataVector& bdat=B.getReady()->getVectorRO();         {
3216          DataVector& rdat=res.getReady()->getVectorRW();           int lerror=0;
3217          const DataTypes::ShapeType& ts=table.getShape();           try
3218          int twx=ts[0]-1;    // table width x           {
3219          int twy=ts[1]-1;    // table width y             do
3220          for (int l=0; l<numpts; ++l)             {
         {  
3221          double a=adat[l];          double a=adat[l];
3222          double b=bdat[l];          double b=bdat[l];
3223          int x=static_cast<int>(((a-Amin)/Astep));          int x=static_cast<int>(((a-Amin)/Astep));
# Line 3193  Data::interpolateFromTable2D(const Wrapp Line 3225  Data::interpolateFromTable2D(const Wrapp
3225                  if (check_boundaries) {                  if (check_boundaries) {
3226                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )
3227                  {                  {
3228                      error=1;                  #pragma omp critical
3229                      break;                    {
3230                        lerror=1;
3231                    }
3232                    break;  
3233                  }                  }
3234                  if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )                  if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )
3235                  {                  {
3236                      error=4;                  #pragma omp critical
3237                      break;                  {
3238                        lerror=4;
3239                    }
3240                    break;
3241                  }                  }
3242                  }                  }
3243                  if (x<0) x=0;                  if (x<0) x=0;
# Line 3212  Data::interpolateFromTable2D(const Wrapp Line 3250  Data::interpolateFromTable2D(const Wrapp
3250                   double sw=table.getElt(x,y);                   double sw=table.getElt(x,y);
3251                   if ((sw>undef))                   if ((sw>undef))
3252                   {                   {
3253                   error=2;                   #pragma omp critical
3254                     {
3255                         lerror=2;
3256                     }
3257                   break;                   break;
3258                   }                   }
3259                   rdat[l]=sw;                   rdat[l]=sw;
# Line 3222  Data::interpolateFromTable2D(const Wrapp Line 3263  Data::interpolateFromTable2D(const Wrapp
3263                   double nw=table.getElt(x,y+1);                   double nw=table.getElt(x,y+1);
3264                   if ((sw>undef) || (nw>undef))                   if ((sw>undef) || (nw>undef))
3265                   {                   {
3266                   error=2;                   #pragma omp critical
3267                     {
3268                        lerror=2;
3269                     }
3270                   break;                   break;
3271                   }                   }
3272                   double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;                   double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
# Line 3235  Data::interpolateFromTable2D(const Wrapp Line 3279  Data::interpolateFromTable2D(const Wrapp
3279                   double se=table.getElt(x+1,y);                   double se=table.getElt(x+1,y);
3280                   if ((sw>undef) || (se>undef) )                   if ((sw>undef) || (se>undef) )
3281                   {                   {
3282                   error=2;                   #pragma omp critical
3283                     {
3284                        lerror=2;
3285                     }
3286                   break;                   break;
3287                   }                   }
3288                   double la = 2.0*(a-Amin-(x*Astep))/Astep-1;                   double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
# Line 3248  Data::interpolateFromTable2D(const Wrapp Line 3295  Data::interpolateFromTable2D(const Wrapp
3295                   double ne=table.getElt(x+1,y+1);                   double ne=table.getElt(x+1,y+1);
3296                   if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))                   if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3297                   {                   {
3298                   error=2;                  #pragma omp critical
3299                   break;                  {
3300                         lerror=2;
3301                    }
3302                    break;
3303                   }                   }
3304                   // map x*Astep <= a << (x+1)*Astep to [-1,1]                   // map x*Astep <= a << (x+1)*Astep to [-1,1]
3305                   // same with b                   // same with b
# Line 3259  Data::interpolateFromTable2D(const Wrapp Line 3309  Data::interpolateFromTable2D(const Wrapp
3309                        (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;                        (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3310                       }                       }
3311                  }                  }
3312          }             } while (false);
3313      } catch (DataException d)            } catch (DataException d)
3314      {            {
3315            lerror=3;
3316              }
3317              if (lerror!=0)
3318              {
3319                #pragma omp critical    // Doco says there is a flush associated with critical
3320                {
3321                    haserror=true;  // We only care that one error is recorded. We don't care which
3322                    error=lerror;   // one
3323                }      
3324              }
3325            }  // if (!error)
3326        }   // parallel for
3327        } catch (DataException d)
3328        {
3329          error=3;          error=3;
3330          break;      }
     }  
     } while (false);  
3331  #ifdef PASO_MPI  #ifdef PASO_MPI
3332      int rerror=0;      int rerror=0;
3333      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );

Legend:
Removed from v.2673  
changed lines
  Added in v.2677

  ViewVC Help
Powered by ViewVC 1.1.26