/[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 2037 by jfenwick, Thu Nov 13 06:17:12 2008 UTC revision 2246 by jfenwick, Thu Feb 5 06:04:55 2009 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 538  Data::copyWithMask(const Data& other, Line 541  Data::copyWithMask(const Data& other,
541    {    {
542      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
543    }    }
544      unsigned int selfrank=getDataPointRank();
545      unsigned int otherrank=other2.getDataPointRank();
546      unsigned int maskrank=mask2.getDataPointRank();
547      if ((selfrank==0) && (otherrank>0 || maskrank>0))
548      {
549        // to get here we must be copying from a large object into a scalar
550        // I am not allowing this.
551        // If you are calling copyWithMask then you are considering keeping some existing values
552        // and so I'm going to assume that you don't want your data objects getting a new shape.
553        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
554      }
555    // Now we iterate over the elements    // Now we iterate over the elements
556    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReadyPtr()->getVector();
557    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVector();
558    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVector();
559    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
560      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
561      {
562        // Not allowing this combination.
563        // it is not clear what the rank of the target should be.
564        // Should it be filled with the scalar (rank stays the same);
565        // or should the target object be reshaped to be a scalar as well.
566        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
567      }
568      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
569      {
570        if (mvec[0]>0)      // copy whole object if scalar is >0
571        {
572            copy(other);
573        }
574        return;
575      }
576      if (isTagged())       // so all objects involved will also be tagged
577      {
578        // note the !
579        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
580            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
581        {
582            throw DataException("copyWithMask, shape mismatch.");
583        }
584    
585        // We need to consider the possibility that tags are missing or in the wrong order
586        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
587        // all tags which are not explicitly defined
588    
589        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
590        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
591        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
592    
593        // first, add any tags missing from other or mask
594        const DataTagged::DataMapType& olookup=optr->getTagLookup();
595            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
596        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
597        for (i=olookup.begin();i!=olookup.end();i++)
598        {
599               tptr->addTag(i->first);
600            }
601            for (i=mlookup.begin();i!=mlookup.end();i++) {
602               tptr->addTag(i->first);
603            }
604        // now we know that *this has all the required tags but they aren't guaranteed to be in
605        // the same order
606    
607        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
608        if ((selfrank==otherrank) && (otherrank==maskrank))
609        {
610            for (i=olookup.begin();i!=olookup.end();i++)
611            {
612                // get the target offset
613                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
614                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
615                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
616                for (int i=0;i<getDataPointSize();++i)
617                {
618                    if (mvec[i+moff]>0)
619                    {
620                        self[i+toff]=ovec[i+ooff];
621                    }
622                }
623                }
624        }
625        else    // other is a scalar
626        {
627            for (i=mlookup.begin();i!=mlookup.end();i++)
628            {
629                // get the target offset
630                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
631                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
632                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
633                for (int i=0;i<getDataPointSize();++i)
634                {
635                    if (mvec[i+moff]>0)
636                    {
637                        self[i+toff]=ovec[ooff];
638                    }
639                }
640                }
641        }
642    
643        return;         // ugly
644      }
645      // mixed scalar and non-scalar operation
646      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
647    {    {
648      throw DataException("Error - size mismatch in arguments to copyWithMask.");          size_t num_points=self.size();
649        // OPENMP 3.0 allows unsigned loop vars.
650        #if defined(_OPENMP) && (_OPENMP < 200805)
651        long i;
652        #else
653        size_t i;
654        #endif      
655        #pragma omp parallel for private(i) schedule(static)
656        for (i=0;i<num_points;++i)
657        {
658            if (mvec[i]>0)
659            {
660                self[i]=ovec[0];
661            }
662        }
663        return;         // ugly!
664      }
665      // tagged data is already taken care of so we only need to worry about shapes
666      // special cases with scalars are already dealt with so all we need to worry about is shape
667      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
668      {
669        ostringstream oss;
670        oss <<"Error - size mismatch in arguments to copyWithMask.";
671        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
672        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
673        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
674        throw DataException(oss.str());
675    }    }
676    size_t num_points=self.size();    size_t num_points=self.size();
677    
# Line 691  Data::resolve() Line 818  Data::resolve()
818  Data  Data
819  Data::oneOver() const  Data::oneOver() const
820  {  {
821    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
822    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
823  }  }
824    
825  Data  Data
826  Data::wherePositive() const  Data::wherePositive() const
827  {  {
828    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
829    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
830  }  }
831    
832  Data  Data
833  Data::whereNegative() const  Data::whereNegative() const
834  {  {
835    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
836    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
837  }  }
838    
839  Data  Data
840  Data::whereNonNegative() const  Data::whereNonNegative() const
841  {  {
842    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
843    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
844  }  }
845    
846  Data  Data
847  Data::whereNonPositive() const  Data::whereNonPositive() const
848  {  {
849    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
850    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
851  }  }
852    
853  Data  Data
854  Data::whereZero(double tol) const  Data::whereZero(double tol) const
855  {  {
856    Data dataAbs=abs();  //   Data dataAbs=abs();
857    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
858       MAKELAZYOPOFF(EZ,tol)
859       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
860    
861  }  }
862    
863  Data  Data
864  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
865  {  {
866    Data dataAbs=abs();  //   Data dataAbs=abs();
867    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
868      MAKELAZYOPOFF(NEZ,tol)
869      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
870    
871  }  }
872    
873  Data  Data
# Line 767  bool Line 880  bool
880  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
881  {  {
882    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());  
 //     }  
 //   }  
883  }  }
884    
885  Data  Data
# Line 911  Data:: getValueOfDataPoint(int dataPoint Line 1014  Data:: getValueOfDataPoint(int dataPoint
1014    //    //
1015    // return the array    // return the array
1016    return numArray;    return numArray;
   
1017  }  }
1018    
1019  void  void
# Line 1151  Data::integrateWorker() const Line 1253  Data::integrateWorker() const
1253    // calculate the integral values    // calculate the integral values
1254    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1255    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1256      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1257      if (dom==0)
1258      {            
1259        throw DataException("Can not integrate over non-continuous domains.");
1260      }
1261  #ifdef PASO_MPI  #ifdef PASO_MPI
1262    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1263    // 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
1264    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1265    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
# Line 1162  Data::integrateWorker() const Line 1269  Data::integrateWorker() const
1269    delete[] tmp;    delete[] tmp;
1270    delete[] tmp_local;    delete[] tmp_local;
1271  #else  #else
1272    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1273  #endif  #endif
1274    
1275    //    //
# Line 1223  Data::integrateWorker() const Line 1330  Data::integrateWorker() const
1330  Data  Data
1331  Data::sin() const  Data::sin() const
1332  {  {
1333    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1334    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1335  }  }
1336    
1337  Data  Data
1338  Data::cos() const  Data::cos() const
1339  {  {
1340    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1341    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1342  }  }
1343    
1344  Data  Data
1345  Data::tan() const  Data::tan() const
1346  {  {
1347    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1348    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1349  }  }
1350    
1351  Data  Data
1352  Data::asin() const  Data::asin() const
1353  {  {
1354    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1355    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1356  }  }
1357    
1358  Data  Data
1359  Data::acos() const  Data::acos() const
1360  {  {
1361    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1362    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1363  }  }
1364    
# Line 1279  Data::acos() const Line 1366  Data::acos() const
1366  Data  Data
1367  Data::atan() const  Data::atan() const
1368  {  {
1369    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1370    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1371  }  }
1372    
1373  Data  Data
1374  Data::sinh() const  Data::sinh() const
1375  {  {
1376    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1377    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1378  }  }
1379    
1380  Data  Data
1381  Data::cosh() const  Data::cosh() const
1382  {  {
1383    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1384    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1385  }  }
1386    
1387  Data  Data
1388  Data::tanh() const  Data::tanh() const
1389  {  {
1390    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1391    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1392  }  }
1393    
# Line 1324  Data::tanh() const Line 1395  Data::tanh() const
1395  Data  Data
1396  Data::erf() const  Data::erf() const
1397  {  {
1398  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1399    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1400  #else  #else
1401    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1402    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1403  #endif  #endif
1404  }  }
# Line 1339  Data::erf() const Line 1406  Data::erf() const
1406  Data  Data
1407  Data::asinh() const  Data::asinh() const
1408  {  {
1409    if (isLazy())    MAKELAZYOP(ASINH)
1410    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1411    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1412  #else  #else
1413    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1354  Data::asinh() const Line 1417  Data::asinh() const
1417  Data  Data
1418  Data::acosh() const  Data::acosh() const
1419  {  {
1420    if (isLazy())    MAKELAZYOP(ACOSH)
1421    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1422    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1423  #else  #else
1424    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1369  Data::acosh() const Line 1428  Data::acosh() const
1428  Data  Data
1429  Data::atanh() const  Data::atanh() const
1430  {  {
1431    if (isLazy())    MAKELAZYOP(ATANH)
1432    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1433    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1434  #else  #else
1435    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1383  Data::atanh() const Line 1438  Data::atanh() const
1438    
1439  Data  Data
1440  Data::log10() const  Data::log10() const
1441  {  if (isLazy())  {
1442    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1443    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1444  }  }
1445    
1446  Data  Data
1447  Data::log() const  Data::log() const
1448  {  {
1449    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1450    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1451  }  }
1452    
1453  Data  Data
1454  Data::sign() const  Data::sign() const
1455  {  {
1456    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1457    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1458  }  }
1459    
1460  Data  Data
1461  Data::abs() const  Data::abs() const
1462  {  {
1463    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1464    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1465  }  }
1466    
1467  Data  Data
1468  Data::neg() const  Data::neg() const
1469  {  {
1470    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1471    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1472  }  }
1473    
# Line 1448  Data::pos() const Line 1484  Data::pos() const
1484    
1485  Data  Data
1486  Data::exp() const  Data::exp() const
1487  {    {
1488    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1489    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1490  }  }
1491    
1492  Data  Data
1493  Data::sqrt() const  Data::sqrt() const
1494  {  {
1495    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1496    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1497  }  }
1498    
# Line 1483  Data::Lsup() Line 1511  Data::Lsup()
1511  {  {
1512     if (isLazy())     if (isLazy())
1513     {     {
1514      expand();      resolve();
1515     }     }
1516     return LsupWorker();     return LsupWorker();
1517  }  }
# Line 1503  Data::sup() Line 1531  Data::sup()
1531  {  {
1532     if (isLazy())     if (isLazy())
1533     {     {
1534      expand();      resolve();
1535     }     }
1536     return supWorker();     return supWorker();
1537  }  }
# Line 1523  Data::inf() Line 1551  Data::inf()
1551  {  {
1552     if (isLazy())     if (isLazy())
1553     {     {
1554      expand();      resolve();
1555     }     }
1556     return infWorker();     return infWorker();
1557  }  }
# Line 1615  Data::minval() const Line 1643  Data::minval() const
1643  Data  Data
1644  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1645  {  {
1646         if (isLazy())
1647         {
1648        Data temp(*this);
1649        temp.resolve();
1650        return temp.swapaxes(axis0,axis1);
1651         }
1652       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1653       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1654       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1706  Data::symmetric() const
1706       else {       else {
1707          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.");
1708       }       }
1709       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1710       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1711       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1712       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1716  Data::symmetric() const
1716  Data  Data
1717  Data::nonsymmetric() const  Data::nonsymmetric() const
1718  {  {
1719       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1720       // check input       // check input
1721       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1722       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1748  Data::nonsymmetric() const
1748       }       }
1749  }  }
1750    
   
 // 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).  
1751  Data  Data
1752  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1753  {  {    
1754       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1755         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1756       {       {
1757      Data temp(*this);   // to get around the fact that you can't resolve a const Data      throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
     temp.resolve();  
     return temp.trace(axis_offset);  
1758       }       }
1759       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1760       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1806  Data::trace(int axis_offset) const
1806  Data  Data
1807  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1808  {      {    
1809       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);  
      }  
1810       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1811       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1812       // 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 1908  Data::calc_minGlobalDataPoint(int& ProcN
1908    double next,local_min;    double next,local_min;
1909    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1910    
1911    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1912    {    {
1913      local_min=min;      local_min=min;
1914      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
# Line 1921  Data::calc_minGlobalDataPoint(int& ProcN Line 1935  Data::calc_minGlobalDataPoint(int& ProcN
1935      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPoint(lowi,lowj);
1936      int lowProc = 0;      int lowProc = 0;
1937      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1938      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1939        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1940    
1941      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1942          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1987  Data::operator+=(const Data& right) Line 2002  Data::operator+=(const Data& right)
2002    if (isProtected()) {    if (isProtected()) {
2003          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2004    }    }
2005    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2006    {    binaryOp(right,plus<double>());
2007      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);  
   }  
2008  }  }
2009    
2010  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2014  Data::operator+=(const boost::python::ob
2014          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2015    }    }
2016    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2017    if (isLazy())    MAKELAZYBINSELF(tmp,ADD)
2018    {    binaryOp(tmp,plus<double>());
2019      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);  
   }  
2020  }  }
2021    
2022  // 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 2033  Data::operator-=(const Data& right)
2033    if (isProtected()) {    if (isProtected()) {
2034          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2035    }    }
2036    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
2037    {    binaryOp(right,minus<double>());
2038      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);  
   }  
2039  }  }
2040    
2041  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2045  Data::operator-=(const boost::python::ob
2045          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2046    }    }
2047    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2048    if (isLazy())    MAKELAZYBINSELF(tmp,SUB)
2049    {    binaryOp(tmp,minus<double>());
2050      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);  
   }  
2051  }  }
2052    
2053  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 2056  Data::operator*=(const Data& right)
2056    if (isProtected()) {    if (isProtected()) {
2057          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2058    }    }
2059    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2060    {    binaryOp(right,multiplies<double>());
2061      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);  
   }  
2062  }  }
2063    
2064  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2068  Data::operator*=(const boost::python::ob
2068          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2069    }    }
2070    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2071    if (isLazy())    MAKELAZYBINSELF(tmp,MUL)
2072    {    binaryOp(tmp,multiplies<double>());
2073      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);  
   }  
2074  }  }
2075    
2076  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2079  Data::operator/=(const Data& right)
2079    if (isProtected()) {    if (isProtected()) {
2080          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2081    }    }
2082    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2083    {    binaryOp(right,divides<double>());
2084      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);  
   }  
2085  }  }
2086    
2087  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2091  Data::operator/=(const boost::python::ob
2091          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2092    }    }
2093    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2094    if (isLazy())    MAKELAZYBINSELF(tmp,DIV)
2095    {    binaryOp(tmp,divides<double>());
2096      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);  
   }  
2097  }  }
2098    
2099  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2113  Data::powO(const boost::python::object&
2113  Data  Data
2114  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2115  {  {
2116    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2117    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2118  }  }
2119    
# Line 2175  Data::powD(const Data& right) const Line 2122  Data::powD(const Data& right) const
2122  Data  Data
2123  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2124  {  {
2125    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2126    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2127  }  }
2128    
# Line 2188  escript::operator+(const Data& left, con Line 2131  escript::operator+(const Data& left, con
2131  Data  Data
2132  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2133  {  {
2134    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2135    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2136  }  }
2137    
# Line 2201  escript::operator-(const Data& left, con Line 2140  escript::operator-(const Data& left, con
2140  Data  Data
2141  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2142  {  {
2143    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2144    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2145  }  }
2146    
# Line 2214  escript::operator*(const Data& left, con Line 2149  escript::operator*(const Data& left, con
2149  Data  Data
2150  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2151  {  {
2152    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2153    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2154  }  }
2155    
# Line 2227  escript::operator/(const Data& left, con Line 2158  escript::operator/(const Data& left, con
2158  Data  Data
2159  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2160  {  {
2161    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2162    {    MAKELAZYBIN2(left,tmp,ADD)
2163      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);  
2164  }  }
2165    
2166  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2168  escript::operator+(const Data& left, con
2168  Data  Data
2169  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2170  {  {
2171    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2172    {    MAKELAZYBIN2(left,tmp,SUB)
2173      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);  
2174  }  }
2175    
2176  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2178  escript::operator-(const Data& left, con
2178  Data  Data
2179  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2180  {  {
2181    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2182    {    MAKELAZYBIN2(left,tmp,MUL)
2183      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);  
2184  }  }
2185    
2186  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2188  escript::operator*(const Data& left, con
2188  Data  Data
2189  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2190  {  {
2191    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2192    {    MAKELAZYBIN2(left,tmp,DIV)
2193      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);  
2194  }  }
2195    
2196  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2198  escript::operator/(const Data& left, con
2198  Data  Data
2199  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2200  {  {
2201    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2202    {    MAKELAZYBIN2(tmp,right,ADD)
2203      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;  
2204  }  }
2205    
2206  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2208  escript::operator+(const boost::python::
2208  Data  Data
2209  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2210  {  {
2211    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2212    {    MAKELAZYBIN2(tmp,right,SUB)
2213      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;  
2214  }  }
2215    
2216  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2218  escript::operator-(const boost::python::
2218  Data  Data
2219  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2220  {  {
2221    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2222    {    MAKELAZYBIN2(tmp,right,MUL)
2223      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;  
2224  }  }
2225    
2226  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2228  escript::operator*(const boost::python::
2228  Data  Data
2229  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2230  {  {
2231    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2232    {    MAKELAZYBIN2(tmp,right,DIV)
2233      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;  
2234  }  }
2235    
2236    
# Line 2461  Data::setTaggedValue(int tagKey, Line 2368  Data::setTaggedValue(int tagKey,
2368    }    }
2369    
2370    DataVector temp_data2;    DataVector temp_data2;
2371    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromNumArray(asNumArray,1);
2372    
2373    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2374  }  }
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2419  escript::C_GeneralTensorProduct(Data& ar
2419    // 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().
2420    
2421    // deal with any lazy data    // deal with any lazy data
2422    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2423    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2424      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2425      {
2426        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2427        return Data(c);
2428      }
2429    
2430    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2431    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2502  escript::C_GeneralTensorProduct(Data& ar
2502       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
2503    }    }
2504    
2505      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2506      {
2507         ostringstream os;
2508         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2509         throw DataException(os.str());
2510      }
2511    
2512    // Declare output Data object    // Declare output Data object
2513    Data res;    Data res;
2514    
# Line 2923  Data::borrowReadyPtr() const Line 2842  Data::borrowReadyPtr() const
2842  std::string  std::string
2843  Data::toString() const  Data::toString() const
2844  {  {
2845      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2846      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2847        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2848      {      {
2849      stringstream temp;      stringstream temp;
2850      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
# Line 3020  Data::dump(const std::string fileName) c Line 2940  Data::dump(const std::string fileName) c
2940            return m_data->dump(fileName);            return m_data->dump(fileName);
2941      }      }
2942    }    }
2943    catch (exception& e)    catch (std::exception& e)
2944    {    {
2945          cout << e.what() << endl;          cout << e.what() << endl;
2946    }    }

Legend:
Removed from v.2037  
changed lines
  Added in v.2246

  ViewVC Help
Powered by ViewVC 1.1.26