/[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 2260 by jfenwick, Tue Feb 10 04:50:10 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      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // 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        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
597        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
598        for (i=olookup.begin();i!=olookup.end();i++)
599        {
600               tptr->addTag(i->first);
601            }
602            for (i=mlookup.begin();i!=mlookup.end();i++) {
603               tptr->addTag(i->first);
604            }
605        // now we know that *this has all the required tags but they aren't guaranteed to be in
606        // the same order
607    
608        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
609        if ((selfrank==otherrank) && (otherrank==maskrank))
610        {
611            for (i=tlookup.begin();i!=tlookup.end();i++)
612            {
613                // get the target offset
614                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
615                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
616                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
617                for (int j=0;j<getDataPointSize();++j)
618                {
619                    if (mvec[j+moff]>0)
620                    {
621                        self[j+toff]=ovec[j+ooff];
622                    }
623                }
624                }
625            // now for the default value
626            for (int j=0;j<getDataPointSize();++j)
627            {
628                if (mvec[j+mptr->getDefaultOffset()]>0)
629                {
630                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
631                }
632            }
633        }
634        else    // other is a scalar
635        {
636            for (i=tlookup.begin();i!=tlookup.end();i++)
637            {
638                // get the target offset
639                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
640                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
641                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
642                for (int j=0;j<getDataPointSize();++j)
643                {
644                    if (mvec[j+moff]>0)
645                    {
646                        self[j+toff]=ovec[ooff];
647                    }
648                }
649                }
650            // now for the default value
651            for (int j=0;j<getDataPointSize();++j)
652            {
653                if (mvec[j+mptr->getDefaultOffset()]>0)
654                {
655                    self[j+tptr->getDefaultOffset()]=ovec[0];
656                }
657            }
658        }
659    
660        return;         // ugly
661      }
662      // mixed scalar and non-scalar operation
663      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
664      {
665            size_t num_points=self.size();
666        // OPENMP 3.0 allows unsigned loop vars.
667        #if defined(_OPENMP) && (_OPENMP < 200805)
668        long i;
669        #else
670        size_t i;
671        #endif  
672        size_t psize=getDataPointSize();    
673        #pragma omp parallel for private(i) schedule(static)
674        for (i=0;i<num_points;++i)
675        {
676            if (mvec[i]>0)
677            {
678                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
679            }                   // dest point
680        }
681        return;         // ugly!
682      }
683      // tagged data is already taken care of so we only need to worry about shapes
684      // special cases with scalars are already dealt with so all we need to worry about is shape
685      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
686      {
687        ostringstream oss;
688        oss <<"Error - size mismatch in arguments to copyWithMask.";
689        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
690        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
691        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
692        throw DataException(oss.str());
693    }    }
694    size_t num_points=self.size();    size_t num_points=self.size();
695    
# Line 691  Data::resolve() Line 836  Data::resolve()
836  Data  Data
837  Data::oneOver() const  Data::oneOver() const
838  {  {
839    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
840    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
841  }  }
842    
843  Data  Data
844  Data::wherePositive() const  Data::wherePositive() const
845  {  {
846    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
847    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
848  }  }
849    
850  Data  Data
851  Data::whereNegative() const  Data::whereNegative() const
852  {  {
853    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
854    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
855  }  }
856    
857  Data  Data
858  Data::whereNonNegative() const  Data::whereNonNegative() const
859  {  {
860    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
861    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
862  }  }
863    
864  Data  Data
865  Data::whereNonPositive() const  Data::whereNonPositive() const
866  {  {
867    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
868    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
869  }  }
870    
871  Data  Data
872  Data::whereZero(double tol) const  Data::whereZero(double tol) const
873  {  {
874    Data dataAbs=abs();  //   Data dataAbs=abs();
875    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
876       MAKELAZYOPOFF(EZ,tol)
877       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
878    
879  }  }
880    
881  Data  Data
882  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
883  {  {
884    Data dataAbs=abs();  //   Data dataAbs=abs();
885    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
886      MAKELAZYOPOFF(NEZ,tol)
887      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
888    
889  }  }
890    
891  Data  Data
# Line 767  bool Line 898  bool
898  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
899  {  {
900    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());  
 //     }  
 //   }  
901  }  }
902    
903  Data  Data
# Line 911  Data:: getValueOfDataPoint(int dataPoint Line 1032  Data:: getValueOfDataPoint(int dataPoint
1032    //    //
1033    // return the array    // return the array
1034    return numArray;    return numArray;
   
1035  }  }
1036    
1037  void  void
# Line 1151  Data::integrateWorker() const Line 1271  Data::integrateWorker() const
1271    // calculate the integral values    // calculate the integral values
1272    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1273    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1274      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1275      if (dom==0)
1276      {            
1277        throw DataException("Can not integrate over non-continuous domains.");
1278      }
1279  #ifdef PASO_MPI  #ifdef PASO_MPI
1280    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1281    // 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
1282    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1283    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
# Line 1162  Data::integrateWorker() const Line 1287  Data::integrateWorker() const
1287    delete[] tmp;    delete[] tmp;
1288    delete[] tmp_local;    delete[] tmp_local;
1289  #else  #else
1290    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1291  #endif  #endif
1292    
1293    //    //
# Line 1223  Data::integrateWorker() const Line 1348  Data::integrateWorker() const
1348  Data  Data
1349  Data::sin() const  Data::sin() const
1350  {  {
1351    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1352    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1353  }  }
1354    
1355  Data  Data
1356  Data::cos() const  Data::cos() const
1357  {  {
1358    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1359    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1360  }  }
1361    
1362  Data  Data
1363  Data::tan() const  Data::tan() const
1364  {  {
1365    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1366    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1367  }  }
1368    
1369  Data  Data
1370  Data::asin() const  Data::asin() const
1371  {  {
1372    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1373    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1374  }  }
1375    
1376  Data  Data
1377  Data::acos() const  Data::acos() const
1378  {  {
1379    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1380    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1381  }  }
1382    
# Line 1279  Data::acos() const Line 1384  Data::acos() const
1384  Data  Data
1385  Data::atan() const  Data::atan() const
1386  {  {
1387    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1388    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1389  }  }
1390    
1391  Data  Data
1392  Data::sinh() const  Data::sinh() const
1393  {  {
1394    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1395    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1396  }  }
1397    
1398  Data  Data
1399  Data::cosh() const  Data::cosh() const
1400  {  {
1401    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1402    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1403  }  }
1404    
1405  Data  Data
1406  Data::tanh() const  Data::tanh() const
1407  {  {
1408    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1409    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1410  }  }
1411    
# Line 1327  Data::erf() const Line 1416  Data::erf() const
1416  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1417    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1418  #else  #else
1419    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1420    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1421  #endif  #endif
1422  }  }
# Line 1339  Data::erf() const Line 1424  Data::erf() const
1424  Data  Data
1425  Data::asinh() const  Data::asinh() const
1426  {  {
1427    if (isLazy())    MAKELAZYOP(ASINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
1428  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1429    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1430  #else  #else
# Line 1354  Data::asinh() const Line 1435  Data::asinh() const
1435  Data  Data
1436  Data::acosh() const  Data::acosh() const
1437  {  {
1438    if (isLazy())    MAKELAZYOP(ACOSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
1439  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1440    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1441  #else  #else
# Line 1369  Data::acosh() const Line 1446  Data::acosh() const
1446  Data  Data
1447  Data::atanh() const  Data::atanh() const
1448  {  {
1449    if (isLazy())    MAKELAZYOP(ATANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
1450  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1451    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1452  #else  #else
# Line 1383  Data::atanh() const Line 1456  Data::atanh() const
1456    
1457  Data  Data
1458  Data::log10() const  Data::log10() const
1459  {  if (isLazy())  {
1460    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1461    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1462  }  }
1463    
1464  Data  Data
1465  Data::log() const  Data::log() const
1466  {  {
1467    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1468    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1469  }  }
1470    
1471  Data  Data
1472  Data::sign() const  Data::sign() const
1473  {  {
1474    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1475    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1476  }  }
1477    
1478  Data  Data
1479  Data::abs() const  Data::abs() const
1480  {  {
1481    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1482    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1483  }  }
1484    
1485  Data  Data
1486  Data::neg() const  Data::neg() const
1487  {  {
1488    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1489    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1490  }  }
1491    
# Line 1448  Data::pos() const Line 1502  Data::pos() const
1502    
1503  Data  Data
1504  Data::exp() const  Data::exp() const
1505  {    {
1506    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1507    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1508  }  }
1509    
1510  Data  Data
1511  Data::sqrt() const  Data::sqrt() const
1512  {  {
1513    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1514    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1515  }  }
1516    
# Line 1483  Data::Lsup() Line 1529  Data::Lsup()
1529  {  {
1530     if (isLazy())     if (isLazy())
1531     {     {
1532      expand();      resolve();
1533     }     }
1534     return LsupWorker();     return LsupWorker();
1535  }  }
# Line 1503  Data::sup() Line 1549  Data::sup()
1549  {  {
1550     if (isLazy())     if (isLazy())
1551     {     {
1552      expand();      resolve();
1553     }     }
1554     return supWorker();     return supWorker();
1555  }  }
# Line 1523  Data::inf() Line 1569  Data::inf()
1569  {  {
1570     if (isLazy())     if (isLazy())
1571     {     {
1572      expand();      resolve();
1573     }     }
1574     return infWorker();     return infWorker();
1575  }  }
# Line 1615  Data::minval() const Line 1661  Data::minval() const
1661  Data  Data
1662  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1663  {  {
1664         if (isLazy())
1665         {
1666        Data temp(*this);
1667        temp.resolve();
1668        return temp.swapaxes(axis0,axis1);
1669         }
1670       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1671       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1672       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1724  Data::symmetric() const
1724       else {       else {
1725          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.");
1726       }       }
1727       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1728       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1729       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1730       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1734  Data::symmetric() const
1734  Data  Data
1735  Data::nonsymmetric() const  Data::nonsymmetric() const
1736  {  {
1737       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1738       // check input       // check input
1739       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1740       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1766  Data::nonsymmetric() const
1766       }       }
1767  }  }
1768    
   
 // 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).  
1769  Data  Data
1770  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1771  {  {    
1772       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1773         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1774       {       {
1775      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);  
1776       }       }
1777       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1778       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1824  Data::trace(int axis_offset) const
1824  Data  Data
1825  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1826  {      {    
1827       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);  
      }  
1828       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1829       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1830       // 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 1926  Data::calc_minGlobalDataPoint(int& ProcN
1926    double next,local_min;    double next,local_min;
1927    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1928    
1929    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1930    {    {
1931      local_min=min;      local_min=min;
1932      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
# Line 1921  Data::calc_minGlobalDataPoint(int& ProcN Line 1953  Data::calc_minGlobalDataPoint(int& ProcN
1953      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPoint(lowi,lowj);
1954      int lowProc = 0;      int lowProc = 0;
1955      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1956      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1957        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1958    
1959      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1960          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1987  Data::operator+=(const Data& right) Line 2020  Data::operator+=(const Data& right)
2020    if (isProtected()) {    if (isProtected()) {
2021          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2022    }    }
2023    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2024    {    binaryOp(right,plus<double>());
2025      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);  
   }  
2026  }  }
2027    
2028  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2032  Data::operator+=(const boost::python::ob
2032          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2033    }    }
2034    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2035    if (isLazy())    MAKELAZYBINSELF(tmp,ADD)
2036    {    binaryOp(tmp,plus<double>());
2037      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);  
   }  
2038  }  }
2039    
2040  // 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 2051  Data::operator-=(const Data& right)
2051    if (isProtected()) {    if (isProtected()) {
2052          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2053    }    }
2054    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
2055    {    binaryOp(right,minus<double>());
2056      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);  
   }  
2057  }  }
2058    
2059  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2063  Data::operator-=(const boost::python::ob
2063          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2064    }    }
2065    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2066    if (isLazy())    MAKELAZYBINSELF(tmp,SUB)
2067    {    binaryOp(tmp,minus<double>());
2068      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);  
   }  
2069  }  }
2070    
2071  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 2074  Data::operator*=(const Data& right)
2074    if (isProtected()) {    if (isProtected()) {
2075          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2076    }    }
2077    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2078    {    binaryOp(right,multiplies<double>());
2079      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);  
   }  
2080  }  }
2081    
2082  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2086  Data::operator*=(const boost::python::ob
2086          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2087    }    }
2088    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2089    if (isLazy())    MAKELAZYBINSELF(tmp,MUL)
2090    {    binaryOp(tmp,multiplies<double>());
2091      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);  
   }  
2092  }  }
2093    
2094  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2097  Data::operator/=(const Data& right)
2097    if (isProtected()) {    if (isProtected()) {
2098          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2099    }    }
2100    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2101    {    binaryOp(right,divides<double>());
2102      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);  
   }  
2103  }  }
2104    
2105  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2109  Data::operator/=(const boost::python::ob
2109          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2110    }    }
2111    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2112    if (isLazy())    MAKELAZYBINSELF(tmp,DIV)
2113    {    binaryOp(tmp,divides<double>());
2114      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);  
   }  
2115  }  }
2116    
2117  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2131  Data::powO(const boost::python::object&
2131  Data  Data
2132  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2133  {  {
2134    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2135    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2136  }  }
2137    
# Line 2175  Data::powD(const Data& right) const Line 2140  Data::powD(const Data& right) const
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,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2144    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2145  }  }
2146    
# Line 2188  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,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2153    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2154  }  }
2155    
# Line 2201  escript::operator-(const Data& left, con Line 2158  escript::operator-(const Data& left, con
2158  Data  Data
2159  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2160  {  {
2161    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2162    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2163  }  }
2164    
# Line 2214  escript::operator*(const Data& left, con Line 2167  escript::operator*(const Data& left, con
2167  Data  Data
2168  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2169  {  {
2170    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2171    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2172  }  }
2173    
# Line 2227  escript::operator/(const Data& left, con Line 2176  escript::operator/(const Data& left, con
2176  Data  Data
2177  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2178  {  {
2179    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2180    {    MAKELAZYBIN2(left,tmp,ADD)
2181      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);  
2182  }  }
2183    
2184  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2186  escript::operator+(const Data& left, con
2186  Data  Data
2187  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2188  {  {
2189    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2190    {    MAKELAZYBIN2(left,tmp,SUB)
2191      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);  
2192  }  }
2193    
2194  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2196  escript::operator-(const Data& left, con
2196  Data  Data
2197  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2198  {  {
2199    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2200    {    MAKELAZYBIN2(left,tmp,MUL)
2201      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);  
2202  }  }
2203    
2204  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2206  escript::operator*(const Data& left, con
2206  Data  Data
2207  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2208  {  {
2209    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2210    {    MAKELAZYBIN2(left,tmp,DIV)
2211      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);  
2212  }  }
2213    
2214  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2216  escript::operator/(const Data& left, con
2216  Data  Data
2217  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2218  {  {
2219    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2220    {    MAKELAZYBIN2(tmp,right,ADD)
2221      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;  
2222  }  }
2223    
2224  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2226  escript::operator+(const boost::python::
2226  Data  Data
2227  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2228  {  {
2229    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2230    {    MAKELAZYBIN2(tmp,right,SUB)
2231      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;  
2232  }  }
2233    
2234  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2236  escript::operator-(const boost::python::
2236  Data  Data
2237  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2238  {  {
2239    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2240    {    MAKELAZYBIN2(tmp,right,MUL)
2241      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;  
2242  }  }
2243    
2244  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2246  escript::operator*(const boost::python::
2246  Data  Data
2247  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2248  {  {
2249    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2250    {    MAKELAZYBIN2(tmp,right,DIV)
2251      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;  
2252  }  }
2253    
2254    
# Line 2461  Data::setTaggedValue(int tagKey, Line 2386  Data::setTaggedValue(int tagKey,
2386    }    }
2387    
2388    DataVector temp_data2;    DataVector temp_data2;
2389    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromNumArray(asNumArray,1);
2390    
2391    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2392  }  }
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2437  escript::C_GeneralTensorProduct(Data& ar
2437    // 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().
2438    
2439    // deal with any lazy data    // deal with any lazy data
2440    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2441    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2442      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2443      {
2444        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2445        return Data(c);
2446      }
2447    
2448    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2449    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2520  escript::C_GeneralTensorProduct(Data& ar
2520       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
2521    }    }
2522    
2523      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2524      {
2525         ostringstream os;
2526         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2527         throw DataException(os.str());
2528      }
2529    
2530    // Declare output Data object    // Declare output Data object
2531    Data res;    Data res;
2532    
# Line 2923  Data::borrowReadyPtr() const Line 2860  Data::borrowReadyPtr() const
2860  std::string  std::string
2861  Data::toString() const  Data::toString() const
2862  {  {
2863      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2864      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2865        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2866      {      {
2867      stringstream temp;      stringstream temp;
2868      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 2958  Data::dump(const std::string fileName) c
2958            return m_data->dump(fileName);            return m_data->dump(fileName);
2959      }      }
2960    }    }
2961    catch (exception& e)    catch (std::exception& e)
2962    {    {
2963          cout << e.what() << endl;          cout << e.what() << endl;
2964    }    }

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

  ViewVC Help
Powered by ViewVC 1.1.26