/[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 2677 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; Line 87  using namespace escript;
87      return Data(c);\      return Data(c);\
88    }    }
89    
90    #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
91    
92  namespace  namespace
93  {  {
94    
# Line 552  Data::copy(const Data& other) Line 554  Data::copy(const Data& other)
554  Data  Data
555  Data::delay()  Data::delay()
556  {  {
557    DataLazy* dl=new DataLazy(m_data);    if (!isLazy())
558    return Data(dl);    {
559          DataLazy* dl=new DataLazy(m_data);
560          return Data(dl);
561      }
562      return *this;
563  }  }
564    
565  void  void
# Line 1199  Data::setValueOfDataPointToArray(int dat Line 1205  Data::setValueOfDataPointToArray(int dat
1205    if (isProtected()) {    if (isProtected()) {
1206          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1207    }    }
   forceResolve();  
1208    
1209    WrappedArray w(obj);    WrappedArray w(obj);
1210    //    //
# Line 1213  Data::setValueOfDataPointToArray(int dat Line 1218  Data::setValueOfDataPointToArray(int dat
1218      if (w.getShape()[i]!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1219         throw DataException("Shape of array does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1220    }    }
1221    
1222      exclusiveWrite();
1223    
1224    //    //
1225    // make sure data is expanded:    // make sure data is expanded:
1226    //    //
# Line 1236  Data::setValueOfDataPoint(int dataPointN Line 1244  Data::setValueOfDataPoint(int dataPointN
1244    }    }
1245    //    //
1246    // make sure data is expanded:    // make sure data is expanded:
1247    forceResolve();    exclusiveWrite();
1248    if (!isExpanded()) {    if (!isExpanded()) {
1249      expand();      expand();
1250    }    }
# Line 1318  Data::integrateToTuple() Line 1326  Data::integrateToTuple()
1326  {  {
1327    if (isLazy())    if (isLazy())
1328    {    {
1329      expand();      expand();   // Can't do a non-resolving version of this without changing the domain object
1330    }    }         // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1331    return integrateWorker();    return integrateWorker();
1332    
1333  }  }
# Line 1546  Data::Lsup() Line 1554  Data::Lsup()
1554  {  {
1555     if (isLazy())     if (isLazy())
1556     {     {
1557      resolve();      if (!actsExpanded() || CHECK_DO_CRES)
1558        {
1559            resolve();
1560        }
1561        else
1562        {
1563    #ifdef PASO_MPI
1564            return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1565    #else
1566            return lazyAlgWorker<AbsMax>(0,0);
1567    #endif
1568        }
1569     }     }
1570     return LsupWorker();     return LsupWorker();
1571  }  }
# Line 1566  Data::sup() Line 1585  Data::sup()
1585  {  {
1586     if (isLazy())     if (isLazy())
1587     {     {
1588      resolve();      if (!actsExpanded() || CHECK_DO_CRES)
1589        {
1590            resolve();
1591        }
1592        else
1593        {
1594    #ifdef PASO_MPI
1595            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1596    #else
1597            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, 0);
1598    #endif
1599        }
1600     }     }
1601     return supWorker();     return supWorker();
1602  }  }
# Line 1586  Data::inf() Line 1616  Data::inf()
1616  {  {
1617     if (isLazy())     if (isLazy())
1618     {     {
1619      resolve();      if (!actsExpanded() || CHECK_DO_CRES)
1620        {
1621            resolve();
1622        }
1623        else
1624        {
1625    #ifdef PASO_MPI
1626            return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
1627    #else
1628            return lazyAlgWorker<FMin>(numeric_limits<double>::max(), 0);
1629    #endif
1630        }
1631     }     }
1632     return infWorker();     return infWorker();
1633  }  }
1634    
1635    template <class BinaryOp>
1636    double
1637    Data::lazyAlgWorker(double init, int mpiop_type)
1638    {
1639       if (!isLazy() || !m_data->actsExpanded())
1640       {
1641        throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1642       }
1643       DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1644       EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1645       BufferGroup* bg=allocSampleBuffer();
1646       double val=init;
1647       int i=0;
1648       const size_t numsamples=getNumSamples();
1649       const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1650       BinaryOp operation;
1651       #pragma omp parallel private(i)
1652       {
1653        double localtot=init;
1654        #pragma omp for schedule(static) private(i)
1655        for (i=0;i<numsamples;++i)
1656        {
1657            size_t roffset=0;
1658            const DataTypes::ValueType* v=dl->resolveSample(*bg, i, roffset);
1659            // Now we have the sample, run operation on all points
1660            for (size_t j=0;j<samplesize;++j)
1661            {
1662            localtot=operation(localtot,(*v)[j+roffset]);
1663            }
1664        }
1665        #pragma omp critical
1666        val=operation(val,localtot);
1667       }
1668       freeSampleBuffer(bg);
1669    #ifdef PASO_MPI
1670       double globalValue;
1671       MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1672       return globalValue;
1673    #else
1674       return val;
1675    #endif
1676    }
1677    
1678  double  double
1679  Data::LsupWorker() const  Data::LsupWorker() const
1680  {  {
# Line 1648  Data::infWorker() const Line 1732  Data::infWorker() const
1732  Data  Data
1733  Data::maxval() const  Data::maxval() const
1734  {  {
1735    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();  
   }  
1736    //    //
1737    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1738    FMax fmax_func;    FMax fmax_func;
# Line 1663  Data::maxval() const Line 1742  Data::maxval() const
1742  Data  Data
1743  Data::minval() const  Data::minval() const
1744  {  {
1745    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();  
   }  
1746    //    //
1747    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1748    FMin fmin_func;    FMin fmin_func;
# Line 1854  Data::transpose(int axis_offset) const Line 1928  Data::transpose(int axis_offset) const
1928          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1929       }       }
1930       for (int i=0; i<rank; i++) {       for (int i=0; i<rank; i++) {
1931    
1932         int index = (axis_offset+i)%rank;         int index = (axis_offset+i)%rank;
1933         ev_shape.push_back(s[index]); // Append to new shape         ev_shape.push_back(s[index]); // Append to new shape
1934       }       }

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

  ViewVC Help
Powered by ViewVC 1.1.26