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

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

Parent Directory | Revision Log | View Patch Patch

-revision 2720 by jfenwick,
Tue Sep 22 00:48:00 2009 UTC
+revision 2721 by jfenwick,
Fri Oct 16 05:40:12 2009 UTC
 Line 87 
 using namespace escript;
      return Data(c);\
    }
+ #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
  namespace
  {
-Line 552 
 Data::copy(const Data& other)
+Line 554 
 Data::copy(const Data& other)
  Data
  Data::delay()
  {
-   DataLazy* dl=new DataLazy(m_data);
+   if (!isLazy())
-   return Data(dl);
+   {
+       DataLazy* dl=new DataLazy(m_data);
+       return Data(dl);
+   }
+   return *this;
  }
  void
-Line 1199 
 Data::setValueOfDataPointToArray(int dat
+Line 1205 
 Data::setValueOfDataPointToArray(int dat
    if (isProtected()) {
          throw DataException("Error - attempt to update protected Data object.");
    }
-   forceResolve();
    WrappedArray w(obj);
    //
-Line 1213 
 Data::setValueOfDataPointToArray(int dat
+Line 1218 
 Data::setValueOfDataPointToArray(int dat
      if (w.getShape()[i]!=getDataPointShape()[i])
         throw DataException("Shape of array does not match Data object rank");
    }
+   exclusiveWrite();
    //
    // make sure data is expanded:
    //
-Line 1236 
 Data::setValueOfDataPoint(int dataPointN
+Line 1244 
 Data::setValueOfDataPoint(int dataPointN
    }
    //
    // make sure data is expanded:
-   forceResolve();
+   exclusiveWrite();
    if (!isExpanded()) {
      expand();
    }
-Line 1318 
 Data::integrateToTuple()
+Line 1326 
 Data::integrateToTuple()
  {
    if (isLazy())
    {
-     expand();
+     expand();   // Can't do a non-resolving version of this without changing the domain object
-   }
+   }         // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
    return integrateWorker();
  }
-Line 1546 
 Data::Lsup()
+Line 1554 
 Data::Lsup()
  {
     if (isLazy())
     {
-     resolve();
+     if (!actsExpanded() || CHECK_DO_CRES)
+     {
+         resolve();
+     }
+     else
+     {
+ #ifdef PASO_MPI
+         return lazyAlgWorker<AbsMax>(0,MPI_MAX);
+ #else
+         return lazyAlgWorker<AbsMax>(0,0);
+ #endif
+     }
     }
     return LsupWorker();
  }
-Line 1566 
 Data::sup()
+Line 1585 
 Data::sup()
  {
     if (isLazy())
     {
-     resolve();
+     if (!actsExpanded() || CHECK_DO_CRES)
+     {
+         resolve();
+     }
+     else
+     {
+ #ifdef PASO_MPI
+         return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
+ #else
+         return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, 0);
+ #endif
+     }
     }
     return supWorker();
  }
-Line 1586 
 Data::inf()
+Line 1616 
 Data::inf()
  {
     if (isLazy())
     {
-     resolve();
+     if (!actsExpanded() || CHECK_DO_CRES)
+     {
+         resolve();
+     }
+     else
+     {
+ #ifdef PASO_MPI
+         return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
+ #else
+         return lazyAlgWorker<FMin>(numeric_limits<double>::max(), 0);
+ #endif
+     }
     }
     return infWorker();
  }
+ template <class BinaryOp>
+ double
+ Data::lazyAlgWorker(double init, int mpiop_type)
+ {
+    if (!isLazy() || !m_data->actsExpanded())
+    {
+     throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
+    }
+    DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
+    EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
+    BufferGroup* bg=allocSampleBuffer();
+    double val=init;
+    int i=0;
+    const size_t numsamples=getNumSamples();
+    const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
+    BinaryOp operation;
+    #pragma omp parallel private(i)
+    {
+     double localtot=init;
+     #pragma omp for schedule(static) private(i)
+     for (i=0;i<numsamples;++i)
+     {
+         size_t roffset=0;
+         const DataTypes::ValueType* v=dl->resolveSample(*bg, i, roffset);
+         // Now we have the sample, run operation on all points
+         for (size_t j=0;j<samplesize;++j)
+         {
+         localtot=operation(localtot,(*v)[j+roffset]);
+         }
+     }
+     #pragma omp critical
+     val=operation(val,localtot);
+    }
+    freeSampleBuffer(bg);
+ #ifdef PASO_MPI
+    double globalValue;
+    MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
+    return globalValue;
+ #else
+    return val;
+ #endif
+ }
  double
  Data::LsupWorker() const
  {
-Line 1648 
 Data::infWorker() const
+Line 1732 
 Data::infWorker() const
  Data
  Data::maxval() const
  {
-   if (isLazy())
+    MAKELAZYOP(MAXVAL)
-   {
-     Data temp(*this);   // to get around the fact that you can't resolve a const Data
-     temp.resolve();
-     return temp.maxval();
-   }
    //
    // set the initial maximum value to min possible double
    FMax fmax_func;
-Line 1663 
 Data::maxval() const
+Line 1742 
 Data::maxval() const
  Data
  Data::minval() const
  {
-   if (isLazy())
+   MAKELAZYOP(MINVAL)
-   {
-     Data temp(*this);   // to get around the fact that you can't resolve a const Data
-     temp.resolve();
-     return temp.minval();
-   }
    //
    // set the initial minimum value to max possible double
    FMin fmin_func;
-Line 1854 
 Data::transpose(int axis_offset) const
+Line 1928 
 Data::transpose(int axis_offset) const
          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
       }
       for (int i=0; i<rank; i++) {
         int index = (axis_offset+i)%rank;
         ev_shape.push_back(s[index]); // Append to new shape
       }

 Legend:



Removed from v.2720
 


changed lines


 
Added in v.2721
 Legend:



Removed from v.2720
 


changed lines


 
Added in v.2721
-Removed from v.2720
+Added in v.2721

	ViewVC Help
Powered by ViewVC 1.1.26