/[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 2049 by phornby, Mon Nov 17 08:54:33 2008 UTC revision 2146 by jfenwick, Wed Dec 10 02:59:46 2008 UTC
# Line 26  Line 26 
26  #include "EscriptParams.h"  #include "EscriptParams.h"
27    
28  extern "C" {  extern "C" {
29  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
30  }  }
31    
32  #include <fstream>  #include <fstream>
33  #include <algorithm>  #include <algorithm>
34  #include <vector>  #include <vector>
35  #include <functional>  #include <functional>
36    #include <sstream>  // so we can throw messages about ranks
37    
38  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
39  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
# Line 46  using namespace escript; Line 47  using namespace escript;
47  // ensure the current object is not a DataLazy  // ensure the current object is not a DataLazy
48  // The idea was that we could add an optional warning whenever a resolve is forced  // The idea was that we could add an optional warning whenever a resolve is forced
49  #define FORCERESOLVE if (isLazy()) {resolve();}  #define FORCERESOLVE if (isLazy()) {resolve();}
50    #define AUTOLAZYON escriptParams.getAUTOLAZY()
51    #define MAKELAZYOP(X)   if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
52      {\
53        DataLazy* c=new DataLazy(borrowDataPtr(),X);\
54        return Data(c);\
55      }
56    #define MAKELAZYOPOFF(X,Y) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
57      {\
58        DataLazy* c=new DataLazy(borrowDataPtr(),X,Y);\
59        return Data(c);\
60      }
61    
62    #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
63      {\
64        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
65            m_data=c->getPtr();\
66        return (*this);\
67      }
68    
69    // like the above but returns a new data rather than *this
70    #define MAKELAZYBIN(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
71      {\
72        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
73        return Data(c);\
74      }
75    
76    #define MAKELAZYBIN2(L,R,X)   if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
77      {\
78        DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
79        return Data(c);\
80      }
81    
82  Data::Data()  Data::Data()
83  {  {
# Line 130  Data::Data(const Data& inData, Line 162  Data::Data(const Data& inData,
162      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space
163        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
164        {           // Even though this is constant, we still need to check whether interpolation is allowed        {           // Even though this is constant, we still need to check whether interpolation is allowed
165      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
166        }        }
167        // if the data is not lazy, this will just be a cast to DataReady        // if the data is not lazy, this will just be a cast to DataReady
168        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
# Line 176  Data::Data(const numeric::array& value, Line 208  Data::Data(const numeric::array& value,
208    initialise(value,what,expanded);    initialise(value,what,expanded);
209    m_protected=false;    m_protected=false;
210  }  }
 /*  
 Data::Data(const DataArrayView& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }*/  
211    
212  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
213           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
# Line 212  Data::Data(const object& value, Line 236  Data::Data(const object& value,
236    
237    // extract the shape of the numarray    // extract the shape of the numarray
238    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);
239  // /*  for (int i=0; i < asNumArray.getrank(); i++) {    if (DataTypes::getRank(tempShape)==0) {
 //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
 //   }*/  
 //   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 // /*  DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(asNumArray);*/  
 //   temp_data.copyFromNumArray(asNumArray);  
   
   //  
   // Create DataConstant using the given value and all other parameters  
   // copied from other. If value is a rank 0 object this Data  
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
240    
241    
242      // get the space for the data vector      // get the space for the data vector
243      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
244      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
245      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromNumArray(asNumArray,1);
246    
247      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
248    
249      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);
     //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());  
 //     initialise(temp2_dataView, other.getFunctionSpace(), false);  
   
250      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
 //     boost::shared_ptr<DataAbstract> sp(t);  
 //     m_data=sp;  
251      m_data=DataAbstract_ptr(t);      m_data=DataAbstract_ptr(t);
252    
253    } else {    } else {
254      //      //
255      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
 //     initialise(temp_dataView, other.getFunctionSpace(), false);  
256      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());
257  //     boost::shared_ptr<DataAbstract> sp(t);  //     boost::shared_ptr<DataAbstract> sp(t);
258  //     m_data=sp;  //     m_data=sp;
# Line 435  Data::getShapeTuple() const Line 438  Data::getShapeTuple() const
438  // It can't work out what type the function is based soley on its name.  // It can't work out what type the function is based soley on its name.
439  // There are ways to fix this involving creating function pointer variables for each form  // There are ways to fix this involving creating function pointer variables for each form
440  // but there doesn't seem to be a need given that the methods have the same name from the python point of view  // but there doesn't seem to be a need given that the methods have the same name from the python point of view
441  Data*  Data
442  Data::copySelf()  Data::copySelf()
443  {  {
444     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
445     return new Data(temp);     return Data(temp);
446  }  }
447    
448  void  void
# Line 691  Data::resolve() Line 694  Data::resolve()
694  Data  Data
695  Data::oneOver() const  Data::oneOver() const
696  {  {
697    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
698    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
699  }  }
700    
701  Data  Data
702  Data::wherePositive() const  Data::wherePositive() const
703  {  {
704    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
705    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
706  }  }
707    
708  Data  Data
709  Data::whereNegative() const  Data::whereNegative() const
710  {  {
711    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
712    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
713  }  }
714    
715  Data  Data
716  Data::whereNonNegative() const  Data::whereNonNegative() const
717  {  {
718    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
719    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
720  }  }
721    
722  Data  Data
723  Data::whereNonPositive() const  Data::whereNonPositive() const
724  {  {
725    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
726    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
727  }  }
728    
# Line 767  bool Line 750  bool
750  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
751  {  {
752    return getFunctionSpace().probeInterpolation(functionspace);    return getFunctionSpace().probeInterpolation(functionspace);
 //   if (getFunctionSpace()==functionspace) {  
 //     return true;  
 //   } else {  
 //     const_Domain_ptr domain=getDomain();  
 //     if  (*domain==*functionspace.getDomain()) {  
 //       return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());  
 //     } else {  
 //       return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());  
 //     }  
 //   }  
753  }  }
754    
755  Data  Data
# Line 911  Data:: getValueOfDataPoint(int dataPoint Line 884  Data:: getValueOfDataPoint(int dataPoint
884    //    //
885    // return the array    // return the array
886    return numArray;    return numArray;
   
887  }  }
888    
889  void  void
# Line 1151  Data::integrateWorker() const Line 1123  Data::integrateWorker() const
1123    // calculate the integral values    // calculate the integral values
1124    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1125    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1126      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1127      if (dom==0)
1128      {            
1129        throw DataException("Can not integrate over non-continuous domains.");
1130      }
1131  #ifdef PASO_MPI  #ifdef PASO_MPI
1132    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1133    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1134    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1135    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
# Line 1162  Data::integrateWorker() const Line 1139  Data::integrateWorker() const
1139    delete[] tmp;    delete[] tmp;
1140    delete[] tmp_local;    delete[] tmp_local;
1141  #else  #else
1142    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1143  #endif  #endif
1144    
1145    //    //
# Line 1223  Data::integrateWorker() const Line 1200  Data::integrateWorker() const
1200  Data  Data
1201  Data::sin() const  Data::sin() const
1202  {  {
1203    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1204    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1205  }  }
1206    
1207  Data  Data
1208  Data::cos() const  Data::cos() const
1209  {  {
1210    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1211    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1212  }  }
1213    
1214  Data  Data
1215  Data::tan() const  Data::tan() const
1216  {  {
1217    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1218    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1219  }  }
1220    
1221  Data  Data
1222  Data::asin() const  Data::asin() const
1223  {  {
1224    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1225    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1226  }  }
1227    
1228  Data  Data
1229  Data::acos() const  Data::acos() const
1230  {  {
1231    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1232    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1233  }  }
1234    
# Line 1279  Data::acos() const Line 1236  Data::acos() const
1236  Data  Data
1237  Data::atan() const  Data::atan() const
1238  {  {
1239    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1240    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1241  }  }
1242    
1243  Data  Data
1244  Data::sinh() const  Data::sinh() const
1245  {  {
1246    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1247    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1248  }  }
1249    
1250  Data  Data
1251  Data::cosh() const  Data::cosh() const
1252  {  {
1253    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1254    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1255  }  }
1256    
1257  Data  Data
1258  Data::tanh() const  Data::tanh() const
1259  {  {
1260    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1261    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1262  }  }
1263    
# Line 1327  Data::erf() const Line 1268  Data::erf() const
1268  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1269    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1270  #else  #else
1271    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1272    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1273  #endif  #endif
1274  }  }
# Line 1339  Data::erf() const Line 1276  Data::erf() const
1276  Data  Data
1277  Data::asinh() const  Data::asinh() const
1278  {  {
1279    if (isLazy())    MAKELAZYOP(ASINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
1280  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1281    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1282  #else  #else
# Line 1354  Data::asinh() const Line 1287  Data::asinh() const
1287  Data  Data
1288  Data::acosh() const  Data::acosh() const
1289  {  {
1290    if (isLazy())    MAKELAZYOP(ACOSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
1291  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1292    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1293  #else  #else
# Line 1369  Data::acosh() const Line 1298  Data::acosh() const
1298  Data  Data
1299  Data::atanh() const  Data::atanh() const
1300  {  {
1301    if (isLazy())    MAKELAZYOP(ATANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
1302  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1303    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1304  #else  #else
# Line 1383  Data::atanh() const Line 1308  Data::atanh() const
1308    
1309  Data  Data
1310  Data::log10() const  Data::log10() const
1311  {  if (isLazy())  {
1312    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1313    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1314  }  }
1315    
1316  Data  Data
1317  Data::log() const  Data::log() const
1318  {  {
1319    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1320    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1321  }  }
1322    
1323  Data  Data
1324  Data::sign() const  Data::sign() const
1325  {  {
1326    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1327    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1328  }  }
1329    
1330  Data  Data
1331  Data::abs() const  Data::abs() const
1332  {  {
1333    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1334    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1335  }  }
1336    
1337  Data  Data
1338  Data::neg() const  Data::neg() const
1339  {  {
1340    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1341    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1342  }  }
1343    
# Line 1448  Data::pos() const Line 1354  Data::pos() const
1354    
1355  Data  Data
1356  Data::exp() const  Data::exp() const
1357  {    {
1358    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1359    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1360  }  }
1361    
1362  Data  Data
1363  Data::sqrt() const  Data::sqrt() const
1364  {  {
1365    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1366    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1367  }  }
1368    
# Line 1483  Data::Lsup() Line 1381  Data::Lsup()
1381  {  {
1382     if (isLazy())     if (isLazy())
1383     {     {
1384      expand();      resolve();
1385     }     }
1386     return LsupWorker();     return LsupWorker();
1387  }  }
# Line 1503  Data::sup() Line 1401  Data::sup()
1401  {  {
1402     if (isLazy())     if (isLazy())
1403     {     {
1404      expand();      resolve();
1405     }     }
1406     return supWorker();     return supWorker();
1407  }  }
# Line 1523  Data::inf() Line 1421  Data::inf()
1421  {  {
1422     if (isLazy())     if (isLazy())
1423     {     {
1424      expand();      resolve();
1425     }     }
1426     return infWorker();     return infWorker();
1427  }  }
# Line 1615  Data::minval() const Line 1513  Data::minval() const
1513  Data  Data
1514  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1515  {  {
1516         if (isLazy())
1517         {
1518        Data temp(*this);
1519        temp.resolve();
1520        return temp.swapaxes(axis0,axis1);
1521         }
1522       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1523       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1524       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1576  Data::symmetric() const
1576       else {       else {
1577          throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");          throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1578       }       }
1579       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1580       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1581       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1582       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1586  Data::symmetric() const
1586  Data  Data
1587  Data::nonsymmetric() const  Data::nonsymmetric() const
1588  {  {
1589       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1590       // check input       // check input
1591       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1592       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1618  Data::nonsymmetric() const
1618       }       }
1619  }  }
1620    
   
 // Doing a lazy version of this would require some thought.  
 // First it needs a parameter (which DataLazy doesn't support at the moment).  
 // (secondly although it does not apply to trace) we can't handle operations which return  
 // multiple results (like eigenvectors_values) or return values of different shapes to their input  
 // (like eigenvalues).  
1621  Data  Data
1622  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1623  {  {    
1624       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
      {  
     Data temp(*this);   // to get around the fact that you can't resolve a const Data  
     temp.resolve();  
     return temp.trace(axis_offset);  
      }  
1625       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1626       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1627          DataTypes::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
# Line 1787  Data::trace(int axis_offset) const Line 1672  Data::trace(int axis_offset) const
1672  Data  Data
1673  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1674  {      {    
1675       if (isLazy())       MAKELAZYOPOFF(TRANS,axis_offset)
      {  
     Data temp(*this);   // to get around the fact that you can't resolve a const Data  
     temp.resolve();  
     return temp.transpose(axis_offset);  
      }  
1676       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1677       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1678       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
# Line 1894  Data::calc_minGlobalDataPoint(int& ProcN Line 1774  Data::calc_minGlobalDataPoint(int& ProcN
1774    double next,local_min;    double next,local_min;
1775    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1776    
1777    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1778    {    {
1779      local_min=min;      local_min=min;
1780      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
# Line 1921  Data::calc_minGlobalDataPoint(int& ProcN Line 1801  Data::calc_minGlobalDataPoint(int& ProcN
1801      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPoint(lowi,lowj);
1802      int lowProc = 0;      int lowProc = 0;
1803      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1804      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1805        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1806    
1807      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1808          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1987  Data::operator+=(const Data& right) Line 1868  Data::operator+=(const Data& right)
1868    if (isProtected()) {    if (isProtected()) {
1869          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1870    }    }
1871    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
1872    {    binaryOp(right,plus<double>());
1873      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
1874  }  }
1875    
1876  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 1880  Data::operator+=(const boost::python::ob
1880          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1881    }    }
1882    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1883    if (isLazy())    MAKELAZYBINSELF(tmp,ADD)
1884    {    binaryOp(tmp,plus<double>());
1885      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD);   // for lazy + is equivalent to +=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,plus<double>());  
     return (*this);  
   }  
1886  }  }
1887    
1888  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
# Line 2034  Data::operator-=(const Data& right) Line 1899  Data::operator-=(const Data& right)
1899    if (isProtected()) {    if (isProtected()) {
1900          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1901    }    }
1902    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
1903    {    binaryOp(right,minus<double>());
1904      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
1905  }  }
1906    
1907  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 1911  Data::operator-=(const boost::python::ob
1911          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1912    }    }
1913    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1914    if (isLazy())    MAKELAZYBINSELF(tmp,SUB)
1915    {    binaryOp(tmp,minus<double>());
1916      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB);   // for lazy - is equivalent to -=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,minus<double>());  
     return (*this);  
   }  
1917  }  }
1918    
1919  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 1922  Data::operator*=(const Data& right)
1922    if (isProtected()) {    if (isProtected()) {
1923          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1924    }    }
1925    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
1926    {    binaryOp(right,multiplies<double>());
1927      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
1928  }  }
1929    
1930  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 1934  Data::operator*=(const boost::python::ob
1934          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1935    }    }
1936    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1937    if (isLazy())    MAKELAZYBINSELF(tmp,MUL)
1938    {    binaryOp(tmp,multiplies<double>());
1939      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL);   // for lazy * is equivalent to *=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,multiplies<double>());  
     return (*this);  
   }  
1940  }  }
1941    
1942  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 1945  Data::operator/=(const Data& right)
1945    if (isProtected()) {    if (isProtected()) {
1946          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1947    }    }
1948    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
1949    {    binaryOp(right,divides<double>());
1950      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
1951  }  }
1952    
1953  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 1957  Data::operator/=(const boost::python::ob
1957          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1958    }    }
1959    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1960    if (isLazy())    MAKELAZYBINSELF(tmp,DIV)
1961    {    binaryOp(tmp,divides<double>());
1962      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV);   // for lazy / is equivalent to /=    return (*this);
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,divides<double>());  
     return (*this);  
   }  
1963  }  }
1964    
1965  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 1979  Data::powO(const boost::python::object&
1979  Data  Data
1980  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1981  {  {
1982    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
1983    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1984  }  }
1985    
# Line 2175  Data::powD(const Data& right) const Line 1988  Data::powD(const Data& right) const
1988  Data  Data
1989  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1990  {  {
1991    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
1992    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
1993  }  }
1994    
# Line 2188  escript::operator+(const Data& left, con Line 1997  escript::operator+(const Data& left, con
1997  Data  Data
1998  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1999  {  {
2000    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2001    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2002  }  }
2003    
# Line 2201  escript::operator-(const Data& left, con Line 2006  escript::operator-(const Data& left, con
2006  Data  Data
2007  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2008  {  {
2009    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2010    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2011  }  }
2012    
# Line 2214  escript::operator*(const Data& left, con Line 2015  escript::operator*(const Data& left, con
2015  Data  Data
2016  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2017  {  {
2018    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2019    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2020  }  }
2021    
# Line 2227  escript::operator/(const Data& left, con Line 2024  escript::operator/(const Data& left, con
2024  Data  Data
2025  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2026  {  {
2027    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2028    {    MAKELAZYBIN2(left,tmp,ADD)
2029      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);    return left+tmp;
     return Data(c);  
   }  
   return left+Data(right,left.getFunctionSpace(),false);  
2030  }  }
2031    
2032  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2034  escript::operator+(const Data& left, con
2034  Data  Data
2035  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2036  {  {
2037    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2038    {    MAKELAZYBIN2(left,tmp,SUB)
2039      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);    return left-tmp;
     return Data(c);  
   }  
   return left-Data(right,left.getFunctionSpace(),false);  
2040  }  }
2041    
2042  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2044  escript::operator-(const Data& left, con
2044  Data  Data
2045  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2046  {  {
2047    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2048    {    MAKELAZYBIN2(left,tmp,MUL)
2049      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);    return left*tmp;
     return Data(c);  
   }  
   return left*Data(right,left.getFunctionSpace(),false);  
2050  }  }
2051    
2052  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2054  escript::operator*(const Data& left, con
2054  Data  Data
2055  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2056  {  {
2057    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2058    {    MAKELAZYBIN2(left,tmp,DIV)
2059      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);    return left/tmp;
     return Data(c);  
   }  
   return left/Data(right,left.getFunctionSpace(),false);  
2060  }  }
2061    
2062  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2064  escript::operator/(const Data& left, con
2064  Data  Data
2065  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2066  {  {
2067    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2068    {    MAKELAZYBIN2(tmp,right,ADD)
2069      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);    return tmp+right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)+right;  
2070  }  }
2071    
2072  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2074  escript::operator+(const boost::python::
2074  Data  Data
2075  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2076  {  {
2077    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2078    {    MAKELAZYBIN2(tmp,right,SUB)
2079      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);    return tmp-right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)-right;  
2080  }  }
2081    
2082  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2084  escript::operator-(const boost::python::
2084  Data  Data
2085  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2086  {  {
2087    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2088    {    MAKELAZYBIN2(tmp,right,MUL)
2089      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);    return tmp*right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)*right;  
2090  }  }
2091    
2092  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2094  escript::operator*(const boost::python::
2094  Data  Data
2095  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2096  {  {
2097    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2098    {    MAKELAZYBIN2(tmp,right,DIV)
2099      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);    return tmp/right;
     return Data(c);  
   }  
   return Data(left,right.getFunctionSpace(),false)/right;  
2100  }  }
2101    
2102    
# Line 2461  Data::setTaggedValue(int tagKey, Line 2234  Data::setTaggedValue(int tagKey,
2234    }    }
2235    
2236    DataVector temp_data2;    DataVector temp_data2;
2237    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromNumArray(asNumArray,1);
2238    
2239    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2240  }  }
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2285  escript::C_GeneralTensorProduct(Data& ar
2285    // SM is the product of the last axis_offset entries in arg_0.getShape().    // SM is the product of the last axis_offset entries in arg_0.getShape().
2286    
2287    // deal with any lazy data    // deal with any lazy data
2288    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2289    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2290      if (arg_0.isLazy() || arg_1.isLazy())
2291      {
2292        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2293        return Data(c);
2294      }
2295    
2296    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2297    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2368  escript::C_GeneralTensorProduct(Data& ar
2368       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2369    }    }
2370    
2371      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2372      {
2373         ostringstream os;
2374         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2375         throw DataException(os.str());
2376      }
2377    
2378    // Declare output Data object    // Declare output Data object
2379    Data res;    Data res;
2380    
# Line 2923  Data::borrowReadyPtr() const Line 2708  Data::borrowReadyPtr() const
2708  std::string  std::string
2709  Data::toString() const  Data::toString() const
2710  {  {
2711      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2712      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2713        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2714      {      {
2715      stringstream temp;      stringstream temp;
2716      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();

Legend:
Removed from v.2049  
changed lines
  Added in v.2146

  ViewVC Help
Powered by ViewVC 1.1.26