/[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 2459 by jfenwick, Thu Jun 4 06:17:54 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>
40  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
41    #include "WrappedArray.h"
42    
43  using namespace std;  using namespace std;
44  using namespace boost::python;  using namespace boost::python;
# Line 45  using namespace escript; Line 47  using namespace escript;
47    
48  // ensure the current object is not a DataLazy  // ensure the current object is not a DataLazy
49  // 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
50  #define FORCERESOLVE if (isLazy()) {resolve();}  // #define forceResolve() if (isLazy()) {#resolve();}
51    
52    #define AUTOLAZYON escriptParams.getAUTOLAZY()
53    #define MAKELAZYOP(X)   if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
54      {\
55        DataLazy* c=new DataLazy(borrowDataPtr(),X);\
56        return Data(c);\
57      }
58    #define MAKELAZYOPOFF(X,Y) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
59      {\
60        DataLazy* c=new DataLazy(borrowDataPtr(),X,Y);\
61        return Data(c);\
62      }
63    
64    #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
65      {\
66        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
67    /*         m_data=c->getPtr();*/     set_m_data(c->getPtr());\
68        return (*this);\
69      }
70    
71    // like the above but returns a new data rather than *this
72    #define MAKELAZYBIN(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
73      {\
74        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
75        return Data(c);\
76      }
77    
78    #define MAKELAZYBIN2(L,R,X)   if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
79      {\
80        DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
81        return Data(c);\
82      }
83    
84    // Do not use the following unless you want to make copies on assignment rather than
85    // share data.  CopyOnWrite should make this unnescessary.
86    // #define ASSIGNMENT_MEANS_DEEPCOPY
87    
88    namespace
89    {
90    
91    template <class ARR>
92    inline
93    boost::python::tuple
94    pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
95    {
96        using namespace boost::python;
97        using boost::python::tuple;
98        using boost::python::list;
99    
100        list l;
101        unsigned int dim=shape[0];
102        for (size_t i=0;i<dim;++i)
103        {
104        l.append(v[i+offset]);
105        }
106        return tuple(l);
107    }
108    
109    template <class ARR>
110    inline
111    boost::python::tuple
112    pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
113    {
114        using namespace boost::python;
115        using boost::python::tuple;
116        using boost::python::list;
117    
118        unsigned int shape0=shape[0];
119        unsigned int shape1=shape[1];
120        list lj;
121        for (size_t j=0;j<shape0;++j)
122        {
123            list li;
124        for (size_t i=0;i<shape1;++i)
125        {
126            li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
127        }
128        lj.append(tuple(li));
129        }
130        return tuple(lj);
131    }
132    
133    template <class ARR>
134    inline
135    boost::python::tuple
136    pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
137    {
138        using namespace boost::python;
139        using boost::python::tuple;
140        using boost::python::list;
141    
142        unsigned int shape0=shape[0];
143        unsigned int shape1=shape[1];
144        unsigned int shape2=shape[2];
145    
146        list lk;
147        for (size_t k=0;k<shape0;++k)
148        {
149            list lj;
150        for (size_t j=0;j<shape1;++j)
151        {
152            list li;
153            for (size_t i=0;i<shape2;++i)
154            {
155                    li.append(v[DataTypes::getRelIndex(shape,k,j,i)]);
156                }
157            lj.append(tuple(li));
158            }
159            lk.append(tuple(lj));
160        }
161        return tuple(lk);
162    }
163    
164    template <class ARR>
165    inline
166    boost::python::tuple
167    pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
168    {
169        using namespace boost::python;
170        using boost::python::tuple;
171        using boost::python::list;
172    
173        unsigned int shape0=shape[0];
174        unsigned int shape1=shape[1];
175        unsigned int shape2=shape[2];
176        unsigned int shape3=shape[3];
177    
178        list ll;
179        for (size_t l=0;l<shape0;++l)
180        {
181            list lk;
182        for (size_t k=0;k<shape1;++k)
183        {
184                list lj;
185                for (size_t j=0;j<shape2;++j)
186                {
187                    list li;
188                    for (size_t i=0;i<shape3;++i)
189                    {
190                        li.append(v[DataTypes::getRelIndex(shape,l,k,j,i)]);
191                    }
192                    lj.append(tuple(li));
193                }
194                lk.append(tuple(lj));
195        }
196            ll.append(tuple(lk));
197        }
198        return tuple(ll);
199    }
200    
201    
202    // This should be safer once the DataC RO changes have been brought in
203    template <class ARR>
204    boost::python::tuple
205    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
206    {
207       using namespace boost::python;
208       using boost::python::list;
209       int rank=shape.size();
210       if (rank==0)
211       {
212        return make_tuple(v[0]);
213       }
214       else if (rank==1)
215       {
216            return pointToTuple1(shape,v,0);
217       }
218       else if (rank==2)
219       {
220        return pointToTuple2(shape,v,0);
221       }
222       else if (rank==3)
223       {
224        return pointToTuple3(shape,v,0);
225       }
226       else if (rank==4)
227       {
228        return pointToTuple4(shape,v,0);
229       }
230       else
231         throw DataException("Unknown rank in pointToTuple.");
232    }
233    
234    }  // anonymous namespace
235    
236  Data::Data()  Data::Data()
237        : m_shared(false), m_lazy(false)
238  {  {
239    //    //
240    // Default data is type DataEmpty    // Default data is type DataEmpty
241    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
242    m_data=temp->getPtr();  //   m_data=temp->getPtr();
243      set_m_data(temp->getPtr());
244    m_protected=false;    m_protected=false;
245  }  }
246    
# Line 60  Data::Data(double value, Line 248  Data::Data(double value,
248             const tuple& shape,             const tuple& shape,
249             const FunctionSpace& what,             const FunctionSpace& what,
250             bool expanded)             bool expanded)
251        : m_shared(false), m_lazy(false)
252  {  {
253    DataTypes::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
254    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
# Line 76  Data::Data(double value, Line 265  Data::Data(double value,
265         const DataTypes::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
266         const FunctionSpace& what,         const FunctionSpace& what,
267             bool expanded)             bool expanded)
268        : m_shared(false), m_lazy(false)
269  {  {
270    int len = DataTypes::noValues(dataPointShape);    int len = DataTypes::noValues(dataPointShape);
   
271    DataVector temp_data(len,value,len);    DataVector temp_data(len,value,len);
 //   DataArrayView temp_dataView(temp_data, dataPointShape);  
   
 //   initialise(temp_dataView, what, expanded);  
272    initialise(temp_data, dataPointShape, what, expanded);    initialise(temp_data, dataPointShape, what, expanded);
   
273    m_protected=false;    m_protected=false;
274  }  }
275    
276  Data::Data(const Data& inData)  Data::Data(const Data& inData)
277        : m_shared(false), m_lazy(false)
278  {  {
279    m_data=inData.m_data;    set_m_data(inData.m_data);
280    m_protected=inData.isProtected();    m_protected=inData.isProtected();
281  }  }
282    
283    
284  Data::Data(const Data& inData,  Data::Data(const Data& inData,
285             const DataTypes::RegionType& region)             const DataTypes::RegionType& region)
286        : m_shared(false), m_lazy(false)
287  {  {
288    DataAbstract_ptr dat=inData.m_data;    DataAbstract_ptr dat=inData.m_data;
289    if (inData.isLazy())    if (inData.isLazy())
# Line 110  Data::Data(const Data& inData, Line 297  Data::Data(const Data& inData,
297    //    //
298    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
299    DataAbstract* tmp = dat->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
300    m_data=DataAbstract_ptr(tmp);    set_m_data(DataAbstract_ptr(tmp));
301    m_protected=false;    m_protected=false;
302    
303  }  }
304    
305  Data::Data(const Data& inData,  Data::Data(const Data& inData,
306             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
307        : m_shared(false), m_lazy(false)
308  {  {
309    if (inData.isEmpty())    if (inData.isEmpty())
310    {    {
311      throw DataException("Error - will not interpolate for instances of DataEmpty.");      throw DataException("Error - will not interpolate for instances of DataEmpty.");
312    }    }
313    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
314      m_data=inData.m_data;      set_m_data(inData.m_data);
315    }    }
316    else    else
317    {    {
# Line 130  Data::Data(const Data& inData, Line 319  Data::Data(const Data& inData,
319      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
320        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
321        {           // 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
322      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
323        }        }
324        // 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
325        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
326        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());          DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
327        m_data=DataAbstract_ptr(dc);  //       m_data=DataAbstract_ptr(dc);
328          set_m_data(DataAbstract_ptr(dc));
329      } else {      } else {
330        Data tmp(0,inData.getDataPointShape(),functionspace,true);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
331        // Note: Must use a reference or pointer to a derived object        // Note: Must use a reference or pointer to a derived object
# Line 149  Data::Data(const Data& inData, Line 339  Data::Data(const Data& inData,
339        } else {        } else {
340          inDataDomain->interpolateACross(tmp,inData);          inDataDomain->interpolateACross(tmp,inData);
341        }        }
342        m_data=tmp.m_data;  //       m_data=tmp.m_data;
343          set_m_data(tmp.m_data);
344      }      }
345    }    }
346    m_protected=false;    m_protected=false;
347  }  }
348    
349  Data::Data(DataAbstract* underlyingdata)  Data::Data(DataAbstract* underlyingdata)
350        : m_shared(false), m_lazy(false)
351  {  {
352  //  m_data=shared_ptr<DataAbstract>(underlyingdata);      set_m_data(underlyingdata->getPtr());
     m_data=underlyingdata->getPtr();  
353      m_protected=false;      m_protected=false;
354  }  }
355    
356  Data::Data(DataAbstract_ptr underlyingdata)  Data::Data(DataAbstract_ptr underlyingdata)
357        : m_shared(false), m_lazy(false)
358  {  {
359      m_data=underlyingdata;      set_m_data(underlyingdata);
360      m_protected=false;      m_protected=false;
361  }  }
362    
   
 Data::Data(const numeric::array& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }  
 /*  
 Data::Data(const DataArrayView& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }*/  
   
363  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
364           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
365                   const FunctionSpace& what,                   const FunctionSpace& what,
366                   bool expanded)                   bool expanded)
367        : m_shared(false), m_lazy(false)
368  {  {
369     initialise(value,shape,what,expanded);     initialise(value,shape,what,expanded);
370     m_protected=false;     m_protected=false;
# Line 198  Data::Data(const DataTypes::ValueType& v Line 374  Data::Data(const DataTypes::ValueType& v
374  Data::Data(const object& value,  Data::Data(const object& value,
375         const FunctionSpace& what,         const FunctionSpace& what,
376             bool expanded)             bool expanded)
377        : m_shared(false), m_lazy(false)
378  {  {
379    numeric::array asNumArray(value);    WrappedArray w(value);
380    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
381    m_protected=false;    m_protected=false;
382  }  }
383    
384    
385  Data::Data(const object& value,  Data::Data(const object& value,
386             const Data& other)             const Data& other)
387        : m_shared(false), m_lazy(false)
388  {  {
389    numeric::array asNumArray(value);    WrappedArray w(value);
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);  
 // /*  for (int i=0; i < asNumArray.getrank(); i++) {  
 //     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.  
390    
391    if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {    // extract the shape of the array
392      const DataTypes::ShapeType& tempShape=w.getShape();
393      if (w.getRank()==0) {
394    
395    
396      // get the space for the data vector      // get the space for the data vector
397      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
398      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
399      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
400    
401      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
402    
403      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);      DataVector temp2_data(len, temp_data[0], len);
     //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());  
 //     initialise(temp2_dataView, other.getFunctionSpace(), false);  
   
404      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
405  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
406  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
407    
408    } else {    } else {
409      //      //
410      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
411  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
412      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
413  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
414    }    }
415    m_protected=false;    m_protected=false;
416  }  }
417    
418  Data::~Data()  Data::~Data()
419  {  {
420      set_m_data(DataAbstract_ptr());
421  }  }
422    
423    
424    // only call in thread safe contexts.
425    // This method should be atomic
426    void Data::set_m_data(DataAbstract_ptr p)
427    {
428      if (m_data.get()!=0)  // release old ownership
429      {
430        m_data->removeOwner(this);
431      }
432      if (p.get()!=0)
433      {
434        m_data=p;
435        m_data->addOwner(this);
436        m_shared=m_data->isShared();
437        m_lazy=m_data->isLazy();
438      }
439    }
440    
441  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
442                   const FunctionSpace& what,                   const FunctionSpace& what,
443                   bool expanded)                   bool expanded)
444  {  {
# Line 277  Data::initialise(const boost::python::nu Line 449  Data::initialise(const boost::python::nu
449    // within the shared_ptr constructor.    // within the shared_ptr constructor.
450    if (expanded) {    if (expanded) {
451      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
452  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
453  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
454    } else {    } else {
455      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
456  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
457  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
458    }    }
459  }  }
460    
# Line 302  Data::initialise(const DataTypes::ValueT Line 472  Data::initialise(const DataTypes::ValueT
472    // within the shared_ptr constructor.    // within the shared_ptr constructor.
473    if (expanded) {    if (expanded) {
474      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
475  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
476  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
477    } else {    } else {
478      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
479  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
480  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
481    }    }
482  }  }
483    
484    
 // void  
 // Data::CompareDebug(const Data& rd)  
 // {  
 //  using namespace std;  
 //  bool mismatch=false;  
 //  std::cout << "Comparing left and right" << endl;  
 //  const DataTagged* left=dynamic_cast<DataTagged*>(m_data.get());  
 //  const DataTagged* right=dynamic_cast<DataTagged*>(rd.m_data.get());  
 //    
 //  if (left==0)  
 //  {  
 //      cout << "left arg is not a DataTagged\n";  
 //      return;  
 //  }  
 //    
 //  if (right==0)  
 //  {  
 //      cout << "right arg is not a DataTagged\n";  
 //      return;  
 //  }  
 //  cout << "Num elements=" << left->getVector().size() << ":" << right->getVector().size() << std::endl;  
 //  cout << "Shapes ";  
 //  if (left->getShape()==right->getShape())  
 //  {  
 //      cout << "ok\n";  
 //  }  
 //  else  
 //  {  
 //      cout << "Problem: shapes do not match\n";  
 //      mismatch=true;  
 //  }  
 //  int lim=left->getVector().size();  
 //  if (right->getVector().size()) lim=right->getVector().size();  
 //  for (int i=0;i<lim;++i)  
 //  {  
 //      if (left->getVector()[i]!=right->getVector()[i])  
 //      {  
 //          cout << "[" << i << "] value mismatch " << left->getVector()[i] << ":" << right->getVector()[i] << endl;  
 //          mismatch=true;  
 //      }  
 //  }  
 //  
 //  // still need to check the tag map  
 //  // also need to watch what is happening to function spaces, are they copied or what?  
 //  
 //  const DataTagged::DataMapType& mapleft=left->getTagLookup();  
 //  const DataTagged::DataMapType& mapright=right->getTagLookup();  
 //  
 //  if (mapleft.size()!=mapright.size())  
 //  {  
 //      cout << "Maps are different sizes " << mapleft.size() << ":" << mapright.size() << endl;  
 //      mismatch=true;  
 //      cout << "Left map\n";  
 //      DataTagged::DataMapType::const_iterator i,j;  
 //      for (i=mapleft.begin();i!=mapleft.end();++i) {  
 //          cout << "(" << i->first << "=>" << i->second << ")\n";  
 //      }  
 //      cout << "Right map\n";  
 //      for (i=mapright.begin();i!=mapright.end();++i) {  
 //          cout << "(" << i->first << "=>" << i->second << ")\n";  
 //      }  
 //      cout << "End map\n";  
 //  
 //  }  
 //  
 //  DataTagged::DataMapType::const_iterator i,j;  
 //  for (i=mapleft.begin(),j=mapright.begin();i!=mapleft.end() && j!=mapright.end();++i,++j) {  
 //     if ((i->first!=j->first) || (i->second!=j->second))  
 //     {  
 //      cout << "(" << i->first << "=>" << i->second << ")";  
 //      cout << ":(" << j->first << "=>" << j->second << ") ";  
 //      mismatch=true;  
 //            }  
 //  }  
 //  if (mismatch)  
 //  {  
 //      cout << "#Mismatch\n";  
 //  }  
 // }  
   
485  escriptDataC  escriptDataC
486  Data::getDataC()  Data::getDataC()
487  {  {
# Line 410  Data::getDataC() const Line 498  Data::getDataC() const
498    return temp;    return temp;
499  }  }
500    
501    size_t
502    Data::getSampleBufferSize() const
503    {
504      return m_data->getSampleBufferSize();
505    }
506    
507    
508  const boost::python::tuple  const boost::python::tuple
509  Data::getShapeTuple() const  Data::getShapeTuple() const
510  {  {
# Line 435  Data::getShapeTuple() const Line 530  Data::getShapeTuple() const
530  // 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.
531  // 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
532  // 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
533  Data*  Data
534  Data::copySelf()  Data::copySelf()
535  {  {
536     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
537     return new Data(temp);     return Data(temp);
538  }  }
539    
540  void  void
# Line 447  Data::copy(const Data& other) Line 542  Data::copy(const Data& other)
542  {  {
543    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
544    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
545    m_data=p;  //   m_data=p;
546      set_m_data(p);
547  }  }
548    
549    
# Line 463  Data::delaySelf() Line 559  Data::delaySelf()
559  {  {
560    if (!isLazy())    if (!isLazy())
561    {    {
562      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
563        set_m_data((new DataLazy(m_data))->getPtr());
564    }    }
565  }  }
566    
567    
568    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
569    // to zero, all the tags from all the DataTags would be in the result.
570    // However since they all have the same value (0) whether they are there or not should not matter.
571    // So I have decided that for all types this method will create a constant 0.
572    // It can be promoted up as required.
573    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
574    // but we can deal with that if it arises.
575    //
576  void  void
577  Data::setToZero()  Data::setToZero()
578  {  {
# Line 474  Data::setToZero() Line 580  Data::setToZero()
580    {    {
581       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
582    }    }
583    m_data->setToZero();    if (isLazy())
584      {
585         DataTypes::ValueType v(getNoValues(),0);
586         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
587         DataLazy* dl=new DataLazy(dc->getPtr());
588         set_m_data(dl->getPtr());
589      }
590      else
591      {
592         exclusiveWrite();
593         m_data->setToZero();
594      }
595  }  }
596    
597    
598  void  void
599  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
600                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 656  Data::copyWithMask(const Data& other,
656    {    {
657      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
658    }    }
659      unsigned int selfrank=getDataPointRank();
660      unsigned int otherrank=other2.getDataPointRank();
661      unsigned int maskrank=mask2.getDataPointRank();
662      if ((selfrank==0) && (otherrank>0 || maskrank>0))
663      {
664        // to get here we must be copying from a large object into a scalar
665        // I am not allowing this.
666        // If you are calling copyWithMask then you are considering keeping some existing values
667        // and so I'm going to assume that you don't want your data objects getting a new shape.
668        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
669      }
670      exclusiveWrite();
671    // Now we iterate over the elements    // Now we iterate over the elements
672    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
673    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
674    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
675    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
676      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
677    {    {
678      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // Not allowing this combination.
679        // it is not clear what the rank of the target should be.
680        // Should it be filled with the scalar (rank stays the same);
681        // or should the target object be reshaped to be a scalar as well.
682        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
683      }
684      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
685      {
686        if (mvec[0]>0)      // copy whole object if scalar is >0
687        {
688            copy(other);
689        }
690        return;
691      }
692      if (isTagged())       // so all objects involved will also be tagged
693      {
694        // note the !
695        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
696            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
697        {
698            throw DataException("copyWithMask, shape mismatch.");
699        }
700    
701        // We need to consider the possibility that tags are missing or in the wrong order
702        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
703        // all tags which are not explicitly defined
704    
705        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
706        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
707        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
708    
709        // first, add any tags missing from other or mask
710        const DataTagged::DataMapType& olookup=optr->getTagLookup();
711            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
712        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
713        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
714        for (i=olookup.begin();i!=olookup.end();i++)
715        {
716               tptr->addTag(i->first);
717            }
718            for (i=mlookup.begin();i!=mlookup.end();i++) {
719               tptr->addTag(i->first);
720            }
721        // now we know that *this has all the required tags but they aren't guaranteed to be in
722        // the same order
723    
724        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
725        if ((selfrank==otherrank) && (otherrank==maskrank))
726        {
727            for (i=tlookup.begin();i!=tlookup.end();i++)
728            {
729                // get the target offset
730                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
731                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
732                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
733                for (int j=0;j<getDataPointSize();++j)
734                {
735                    if (mvec[j+moff]>0)
736                    {
737                        self[j+toff]=ovec[j+ooff];
738                    }
739                }
740                }
741            // now for the default value
742            for (int j=0;j<getDataPointSize();++j)
743            {
744                if (mvec[j+mptr->getDefaultOffset()]>0)
745                {
746                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
747                }
748            }
749        }
750        else    // other is a scalar
751        {
752            for (i=tlookup.begin();i!=tlookup.end();i++)
753            {
754                // get the target offset
755                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
756                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
757                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
758                for (int j=0;j<getDataPointSize();++j)
759                {
760                    if (mvec[j+moff]>0)
761                    {
762                        self[j+toff]=ovec[ooff];
763                    }
764                }
765                }
766            // now for the default value
767            for (int j=0;j<getDataPointSize();++j)
768            {
769                if (mvec[j+mptr->getDefaultOffset()]>0)
770                {
771                    self[j+tptr->getDefaultOffset()]=ovec[0];
772                }
773            }
774        }
775    
776        return;         // ugly
777      }
778      // mixed scalar and non-scalar operation
779      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
780      {
781            size_t num_points=self.size();
782        // OPENMP 3.0 allows unsigned loop vars.
783        #if defined(_OPENMP) && (_OPENMP < 200805)
784        long i;
785        #else
786        size_t i;
787        #endif  
788        size_t psize=getDataPointSize();    
789        #pragma omp parallel for private(i) schedule(static)
790        for (i=0;i<num_points;++i)
791        {
792            if (mvec[i]>0)
793            {
794                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
795            }                   // dest point
796        }
797        return;         // ugly!
798      }
799      // tagged data is already taken care of so we only need to worry about shapes
800      // special cases with scalars are already dealt with so all we need to worry about is shape
801      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
802      {
803        ostringstream oss;
804        oss <<"Error - size mismatch in arguments to copyWithMask.";
805        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
806        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
807        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
808        throw DataException(oss.str());
809    }    }
810    size_t num_points=self.size();    size_t num_points=self.size();
811    
# Line 564  Data::copyWithMask(const Data& other, Line 825  Data::copyWithMask(const Data& other,
825    }    }
826  }  }
827    
   
   
828  bool  bool
829  Data::isExpanded() const  Data::isExpanded() const
830  {  {
# Line 574  Data::isExpanded() const Line 833  Data::isExpanded() const
833  }  }
834    
835  bool  bool
836    Data::actsExpanded() const
837    {
838      return m_data->actsExpanded();
839    
840    }
841    
842    
843    bool
844  Data::isTagged() const  Data::isTagged() const
845  {  {
846    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 864  Data::isConstant() const
864  bool  bool
865  Data::isLazy() const  Data::isLazy() const
866  {  {
867    return m_data->isLazy();    return m_lazy;    // not asking m_data because we need to be able to ask this while m_data is changing
868  }  }
869    
870  // at the moment this is synonymous with !isLazy() but that could change  // at the moment this is synonymous with !isLazy() but that could change
# Line 628  Data::expand() Line 895  Data::expand()
895    if (isConstant()) {    if (isConstant()) {
896      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
897      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
898  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
899  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
900    } else if (isTagged()) {    } else if (isTagged()) {
901      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
902      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
903  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
904  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
905    } else if (isExpanded()) {    } else if (isExpanded()) {
906      //      //
907      // do nothing      // do nothing
# Line 656  Data::tag() Line 921  Data::tag()
921    if (isConstant()) {    if (isConstant()) {
922      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
923      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
924  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
925  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
926    } else if (isTagged()) {    } else if (isTagged()) {
927      // do nothing      // do nothing
928    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 935  Data::tag()
935       {       {
936      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
937       }       }
938       m_data=res;      //      m_data=res;
939         set_m_data(res);
940       tag();       tag();
941    } else {    } else {
942      throw DataException("Error - Tagging not implemented for this Data type.");      throw DataException("Error - Tagging not implemented for this Data type.");
# Line 683  Data::resolve() Line 948  Data::resolve()
948  {  {
949    if (isLazy())    if (isLazy())
950    {    {
951       m_data=m_data->resolve();  //      m_data=m_data->resolve();
952        set_m_data(m_data->resolve());
953    }    }
954  }  }
955    
956    void
957    Data::requireWrite()
958    {
959      resolve();
960      exclusiveWrite();
961    }
962    
963  Data  Data
964  Data::oneOver() const  Data::oneOver() const
965  {  {
966    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
967    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
968  }  }
969    
970  Data  Data
971  Data::wherePositive() const  Data::wherePositive() const
972  {  {
973    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
974    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
975  }  }
976    
977  Data  Data
978  Data::whereNegative() const  Data::whereNegative() const
979  {  {
980    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
981    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
982  }  }
983    
984  Data  Data
985  Data::whereNonNegative() const  Data::whereNonNegative() const
986  {  {
987    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
988    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
989  }  }
990    
991  Data  Data
992  Data::whereNonPositive() const  Data::whereNonPositive() const
993  {  {
994    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
995    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
996  }  }
997    
998  Data  Data
999  Data::whereZero(double tol) const  Data::whereZero(double tol) const
1000  {  {
1001    Data dataAbs=abs();  //   Data dataAbs=abs();
1002    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1003       MAKELAZYOPOFF(EZ,tol)
1004       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1005    
1006  }  }
1007    
1008  Data  Data
1009  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
1010  {  {
1011    Data dataAbs=abs();  //   Data dataAbs=abs();
1012    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1013      MAKELAZYOPOFF(NEZ,tol)
1014      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1015    
1016  }  }
1017    
1018  Data  Data
# Line 767  bool Line 1025  bool
1025  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
1026  {  {
1027    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());  
 //     }  
 //   }  
1028  }  }
1029    
1030  Data  Data
# Line 813  Data::getDataPointSize() const Line 1061  Data::getDataPointSize() const
1061    return m_data->getNoValues();    return m_data->getNoValues();
1062  }  }
1063    
1064    
1065  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1066  Data::getLength() const  Data::getLength() const
1067  {  {
1068    return m_data->getLength();    return m_data->getLength();
1069  }  }
1070    
1071  const  // does not return tuples for scalars
1072  boost::python::numeric::array  // There is no parallelism here ... elements need to be added in the correct order.
1073  Data:: getValueOfDataPoint(int dataPointNo)  //   If we could presize the list and then fill in the elements it might work
1074    //   This would need setting elements to be threadsafe.
1075    //   Having mulitple C threads calling into one interpreter is aparently a no-no.
1076    const boost::python::object
1077    Data::toListOfTuples(bool scalarastuple)
1078  {  {
1079    int i, j, k, l;      using namespace boost::python;
1080        if (get_MPISize()>1)
1081    FORCERESOLVE;      {
1082            throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1083    //      }
1084    // determine the rank and shape of each data point      unsigned int rank=getDataPointRank();
1085    int dataPointRank = getDataPointRank();      unsigned int size=getDataPointSize();
1086    const DataTypes::ShapeType& dataPointShape = getDataPointShape();      boost::python::list res;
1087        int npoints=getNumDataPoints();
1088    //      expand();           // This will also resolve if required
1089    // create the numeric array to be returned      const DataTypes::ValueType& vec=getReady()->getVectorRO();
1090    boost::python::numeric::array numArray(0.0);      if (rank==0)
1091        {
1092    //          long count;
1093    // the shape of the returned numeric array will be the same          if (scalarastuple)
1094    // as that of the data point          {
1095    int arrayRank = dataPointRank;              for (count=0;count<npoints;++count)
1096    const DataTypes::ShapeType& arrayShape = dataPointShape;              {
1097                    res.append(make_tuple(vec[count]));
1098    //              }
1099    // resize the numeric array to the shape just calculated          }
1100    if (arrayRank==0) {          else
1101      numArray.resize(1);          {
1102    }              for (count=0;count<npoints;++count)
1103    if (arrayRank==1) {              {
1104      numArray.resize(arrayShape[0]);                  res.append(vec[count]);
1105    }              }
1106    if (arrayRank==2) {          }
1107      numArray.resize(arrayShape[0],arrayShape[1]);      }
1108    }      else if (rank==1)
1109    if (arrayRank==3) {      {
1110      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);          size_t count;
1111    }          size_t offset=0;
1112    if (arrayRank==4) {          for (count=0;count<npoints;++count,offset+=size)
1113      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);          {
1114    }          // need to pull a pointer to the start of the element (also need a way to skip to the next element)
1115                res.append(pointToTuple1(getDataPointShape(), vec, offset));
1116            }
1117        }
1118        else if (rank==2)
1119        {
1120            size_t count;
1121            size_t offset=0;
1122            for (count=0;count<npoints;++count,offset+=size)
1123            {
1124            // need to pull a pointer to the start of the element (also need a way to skip to the next element)
1125                res.append(pointToTuple2(getDataPointShape(), vec, offset));
1126            }
1127        }
1128        else if (rank==3)
1129        {
1130            size_t count;
1131            size_t offset=0;
1132            for (count=0;count<npoints;++count,offset+=size)
1133            {
1134            // need to pull a pointer to the start of the element (also need a way to skip to the next element)
1135                res.append(pointToTuple3(getDataPointShape(), vec, offset));
1136            }
1137        }
1138        else if (rank==4)
1139        {
1140            size_t count;
1141            size_t offset=0;
1142            for (count=0;count<npoints;++count,offset+=size)
1143            {
1144            // need to pull a pointer to the start of the element (also need a way to skip to the next element)
1145                res.append(pointToTuple4(getDataPointShape(), vec, offset));
1146            }
1147        }
1148        else
1149        {
1150            throw DataException("Unknown rank in ::toListOfTuples()");
1151        }
1152        return res;
1153    }
1154    
1155    if (getNumDataPointsPerSample()>0) {  const boost::python::object
1156    Data::getValueOfDataPointAsTuple(int dataPointNo)
1157    {
1158       forceResolve();
1159       if (getNumDataPointsPerSample()>0) {
1160         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1161         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1162         //         //
1163         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1164         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1165             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1166         }         }
1167    
1168         //         //
1169         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1170         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1171             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1172         }         }
1173         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1174         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1175           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1176      }
1177         switch( dataPointRank ){    else
1178              case 0 :    {
1179                  numArray[0] = getDataAtOffset(offset);      // The old method would return an empty numArray of the given shape
1180                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1181              case 1 :      throw DataException("Error - need at least 1 datapoint per sample.");
                 for( i=0; i<dataPointShape[0]; i++ )  
                     numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));  
                 break;  
             case 2 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++)  
                         numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));  
                 break;  
             case 3 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++)  
                             numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));  
                 break;  
             case 4 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++ )  
                             for( l=0; l<dataPointShape[3]; l++)  
                                 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));  
                 break;  
     }  
1182    }    }
   //  
   // return the array  
   return numArray;  
1183    
1184  }  }
1185    
1186    
1187  void  void
1188  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1189  {  {
1190      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1191      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1192  }  }
1193    
1194  void  void
1195  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1196  {  {
1197    if (isProtected()) {    if (isProtected()) {
1198          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1199    }    }
1200    FORCERESOLVE;    forceResolve();
1201    
1202      WrappedArray w(obj);
1203    //    //
1204    // check rank    // check rank
1205    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1206        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1207    
1208    //    //
1209    // check shape of num_array    // check shape of array
1210    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1211      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1212         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1213    }    }
1214    //    //
1215    // make sure data is expanded:    // make sure data is expanded:
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1220  Data::setValueOfDataPointToArray(int dat
1220    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1221         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1222         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1223         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1224    } else {    } else {
1225         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1226    }    }
1227  }  }
1228    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1234  Data::setValueOfDataPoint(int dataPointN
1234    }    }
1235    //    //
1236    // make sure data is expanded:    // make sure data is expanded:
1237    FORCERESOLVE;    forceResolve();
1238    if (!isExpanded()) {    if (!isExpanded()) {
1239      expand();      expand();
1240    }    }
# Line 977  Data::setValueOfDataPoint(int dataPointN Line 1248  Data::setValueOfDataPoint(int dataPointN
1248  }  }
1249    
1250  const  const
1251  boost::python::numeric::array  boost::python::object
1252  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1253  {  {
1254    size_t length=0;    // This could be lazier than it is now
1255    int i, j, k, l, pos;    forceResolve();
   FORCERESOLVE;  
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
1256    
1257    //    // copy datapoint into a buffer
1258    // create the numeric array to be returned    // broadcast buffer to all nodes
1259    boost::python::numeric::array numArray(0.0);    // convert buffer to tuple
1260      // return tuple
   //  
   // the shape of the returned numeric array will be the same  
   // as that of the data point  
   int arrayRank = dataPointRank;  
   const DataTypes::ShapeType& arrayShape = dataPointShape;  
1261    
1262    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1263    // resize the numeric array to the shape just calculated    size_t length=DataTypes::noValues(dataPointShape);
   if (arrayRank==0) {  
     numArray.resize(1);  
   }  
   if (arrayRank==1) {  
     numArray.resize(arrayShape[0]);  
   }  
   if (arrayRank==2) {  
     numArray.resize(arrayShape[0],arrayShape[1]);  
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
   }  
1264    
1265    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1266    double *tmpData = new double[length];    double *tmpData = new double[length];
1267    
1268    //    // updated for the MPI case
1269    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1270          if (getNumDataPointsPerSample()>0) {
1271      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1272      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1273               if (getNumDataPointsPerSample()>0) {      //
1274                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1275                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1276                  //          throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
                 // Check a valid sample number has been supplied  
                 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {  
                   throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
                 }  
   
                 //  
                 // Check a valid data point number has been supplied  
                 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {  
                   throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");  
                 }  
                 // TODO: global error handling  
         // create a view of the data if it is stored locally  
         //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);  
   
         // pack the data from the view into tmpData for MPI communication  
         pos=0;  
         switch( dataPointRank ){  
             case 0 :  
                 tmpData[0] = getDataAtOffset(offset);  
                 break;  
             case 1 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));  
                 break;  
             case 2 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++, pos++ )  
                         tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));  
                 break;  
             case 3 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++, pos++ )  
                             tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));  
                 break;  
             case 4 :  
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++ )  
                         for( k=0; k<dataPointShape[2]; k++ )  
                             for( l=0; l<dataPointShape[3]; l++, pos++ )  
                                 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));  
                 break;  
         }  
             }  
1277      }      }
1278          #ifdef PASO_MPI  
1279          // broadcast the data to all other processes      //
1280      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1281          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1282            throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
     // unpack the data  
     switch( dataPointRank ){  
         case 0 :  
             numArray[0]=tmpData[0];  
             break;  
         case 1 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 numArray[i]=tmpData[i];  
             break;  
         case 2 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                    numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];  
             break;  
         case 3 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                     for( k=0; k<dataPointShape[2]; k++ )  
                         numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];  
             break;  
         case 4 :  
             for( i=0; i<dataPointShape[0]; i++ )  
                 for( j=0; j<dataPointShape[1]; j++ )  
                     for( k=0; k<dataPointShape[2]; k++ )  
                         for( l=0; l<dataPointShape[3]; l++ )  
                                 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];  
             break;  
1283      }      }
1284        // TODO: global error handling
1285        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1286    
1287      delete [] tmpData;      memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1288         }
1289      }
1290    #ifdef PASO_MPI
1291      // broadcast the data to all other processes
1292      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1293    #endif
1294    
1295      boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1296      delete [] tmpData;
1297    //    //
1298    // return the loaded array    // return the loaded array
1299    return numArray;    return t;
1300    
1301  }  }
1302    
1303    
1304  boost::python::numeric::array  boost::python::object
1305  Data::integrate_const() const  Data::integrateToTuple_const() const
1306  {  {
1307    if (isLazy())    if (isLazy())
1308    {    {
# Line 1127  Data::integrate_const() const Line 1311  Data::integrate_const() const
1311    return integrateWorker();    return integrateWorker();
1312  }  }
1313    
1314  boost::python::numeric::array  boost::python::object
1315  Data::integrate()  Data::integrateToTuple()
1316  {  {
1317    if (isLazy())    if (isLazy())
1318    {    {
1319      expand();      expand();
1320    }    }
1321    return integrateWorker();    return integrateWorker();
 }  
   
1322    
1323    }
1324    
1325  boost::python::numeric::array  boost::python::object
1326  Data::integrateWorker() const  Data::integrateWorker() const
1327  {  {
   int index;  
   int rank = getDataPointRank();  
1328    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1329    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1330    
# Line 1151  Data::integrateWorker() const Line 1332  Data::integrateWorker() const
1332    // calculate the integral values    // calculate the integral values
1333    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1334    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1335      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1336      if (dom==0)
1337      {            
1338        throw DataException("Can not integrate over non-continuous domains.");
1339      }
1340  #ifdef PASO_MPI  #ifdef PASO_MPI
1341    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1342    // 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
1343    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1344    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1345    for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }    for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1346    MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );    MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1347    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1348      tuple result=pointToTuple(shape,tmp);
1349    delete[] tmp;    delete[] tmp;
1350    delete[] tmp_local;    delete[] tmp_local;
1351  #else  #else
1352    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1353    /*  double *tmp = new double[dataPointSize];
1354      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1355      tuple result=pointToTuple(shape,integrals);
1356    //   delete tmp;
1357  #endif  #endif
1358    
   //  
   // create the numeric array to be returned  
   // and load the array with the integral values  
   boost::python::numeric::array bp_array(1.0);  
   if (rank==0) {  
     bp_array.resize(1);  
     index = 0;  
     bp_array[0] = integrals[index];  
   }  
   if (rank==1) {  
     bp_array.resize(shape[0]);  
     for (int i=0; i<shape[0]; i++) {  
       index = i;  
       bp_array[i] = integrals[index];  
     }  
   }  
   if (rank==2) {  
        bp_array.resize(shape[0],shape[1]);  
        for (int i=0; i<shape[0]; i++) {  
          for (int j=0; j<shape[1]; j++) {  
            index = i + shape[0] * j;  
            bp_array[make_tuple(i,j)] = integrals[index];  
          }  
        }  
   }  
   if (rank==3) {  
     bp_array.resize(shape[0],shape[1],shape[2]);  
     for (int i=0; i<shape[0]; i++) {  
       for (int j=0; j<shape[1]; j++) {  
         for (int k=0; k<shape[2]; k++) {  
           index = i + shape[0] * ( j + shape[1] * k );  
           bp_array[make_tuple(i,j,k)] = integrals[index];  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);  
     for (int i=0; i<shape[0]; i++) {  
       for (int j=0; j<shape[1]; j++) {  
         for (int k=0; k<shape[2]; k++) {  
           for (int l=0; l<shape[3]; l++) {  
             index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );  
             bp_array[make_tuple(i,j,k,l)] = integrals[index];  
           }  
         }  
       }  
     }  
   }  
1359    
1360    //    return result;
   // return the loaded array  
   return bp_array;  
1361  }  }
1362    
1363  Data  Data
1364  Data::sin() const  Data::sin() const
1365  {  {
1366    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1367    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1368  }  }
1369    
1370  Data  Data
1371  Data::cos() const  Data::cos() const
1372  {  {
1373    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1374    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1375  }  }
1376    
1377  Data  Data
1378  Data::tan() const  Data::tan() const
1379  {  {
1380    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1381    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1382  }  }
1383    
1384  Data  Data
1385  Data::asin() const  Data::asin() const
1386  {  {
1387    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1388    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1389  }  }
1390    
1391  Data  Data
1392  Data::acos() const  Data::acos() const
1393  {  {
1394    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1395    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1396  }  }
1397    
# Line 1279  Data::acos() const Line 1399  Data::acos() const
1399  Data  Data
1400  Data::atan() const  Data::atan() const
1401  {  {
1402    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1403    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1404  }  }
1405    
1406  Data  Data
1407  Data::sinh() const  Data::sinh() const
1408  {  {
1409    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1410    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1411  }  }
1412    
1413  Data  Data
1414  Data::cosh() const  Data::cosh() const
1415  {  {
1416    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1417    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1418  }  }
1419    
1420  Data  Data
1421  Data::tanh() const  Data::tanh() const
1422  {  {
1423    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1424    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1425  }  }
1426    
# Line 1327  Data::erf() const Line 1431  Data::erf() const
1431  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1432    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1433  #else  #else
1434    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1435    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1436  #endif  #endif
1437  }  }
# Line 1339  Data::erf() const Line 1439  Data::erf() const
1439  Data  Data
1440  Data::asinh() const  Data::asinh() const
1441  {  {
1442    if (isLazy())    MAKELAZYOP(ASINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
1443  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1444    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1445  #else  #else
# Line 1354  Data::asinh() const Line 1450  Data::asinh() const
1450  Data  Data
1451  Data::acosh() const  Data::acosh() const
1452  {  {
1453    if (isLazy())    MAKELAZYOP(ACOSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
1454  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1455    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1456  #else  #else
# Line 1369  Data::acosh() const Line 1461  Data::acosh() const
1461  Data  Data
1462  Data::atanh() const  Data::atanh() const
1463  {  {
1464    if (isLazy())    MAKELAZYOP(ATANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
1465  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1466    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1467  #else  #else
# Line 1383  Data::atanh() const Line 1471  Data::atanh() const
1471    
1472  Data  Data
1473  Data::log10() const  Data::log10() const
1474  {  if (isLazy())  {
1475    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1476    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1477  }  }
1478    
1479  Data  Data
1480  Data::log() const  Data::log() const
1481  {  {
1482    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1483    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1484  }  }
1485    
1486  Data  Data
1487  Data::sign() const  Data::sign() const
1488  {  {
1489    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1490    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1491  }  }
1492    
1493  Data  Data
1494  Data::abs() const  Data::abs() const
1495  {  {
1496    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1497    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1498  }  }
1499    
1500  Data  Data
1501  Data::neg() const  Data::neg() const
1502  {  {
1503    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1504    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1505  }  }
1506    
# Line 1448  Data::pos() const Line 1517  Data::pos() const
1517    
1518  Data  Data
1519  Data::exp() const  Data::exp() const
1520  {    {
1521    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1522    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1523  }  }
1524    
1525  Data  Data
1526  Data::sqrt() const  Data::sqrt() const
1527  {  {
1528    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1529    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1530  }  }
1531    
# Line 1483  Data::Lsup() Line 1544  Data::Lsup()
1544  {  {
1545     if (isLazy())     if (isLazy())
1546     {     {
1547      expand();      resolve();
1548     }     }
1549     return LsupWorker();     return LsupWorker();
1550  }  }
# Line 1503  Data::sup() Line 1564  Data::sup()
1564  {  {
1565     if (isLazy())     if (isLazy())
1566     {     {
1567      expand();      resolve();
1568     }     }
1569     return supWorker();     return supWorker();
1570  }  }
# Line 1523  Data::inf() Line 1584  Data::inf()
1584  {  {
1585     if (isLazy())     if (isLazy())
1586     {     {
1587      expand();      resolve();
1588     }     }
1589     return infWorker();     return infWorker();
1590  }  }
# Line 1615  Data::minval() const Line 1676  Data::minval() const
1676  Data  Data
1677  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1678  {  {
1679         if (isLazy())
1680         {
1681        Data temp(*this);
1682        temp.resolve();
1683        return temp.swapaxes(axis0,axis1);
1684         }
1685       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1686       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1687       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1739  Data::symmetric() const
1739       else {       else {
1740          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.");
1741       }       }
1742       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1743       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1744       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1745       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1749  Data::symmetric() const
1749  Data  Data
1750  Data::nonsymmetric() const  Data::nonsymmetric() const
1751  {  {
1752       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1753       // check input       // check input
1754       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1755       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1781  Data::nonsymmetric() const
1781       }       }
1782  }  }
1783    
   
 // 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).  
1784  Data  Data
1785  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1786  {  {    
1787       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1788         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1789       {       {
1790      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);  
1791       }       }
1792       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1793       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1839  Data::trace(int axis_offset) const
1839  Data  Data
1840  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1841  {      {    
1842       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);  
      }  
1843       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1844       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1845       // 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 1941  Data::calc_minGlobalDataPoint(int& ProcN
1941    double next,local_min;    double next,local_min;
1942    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1943    
1944    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1945    {    {
1946      local_min=min;      local_min=min;
1947      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1948      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1949        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1950          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1951          if (next<local_min) {          if (next<local_min) {
1952            local_min=next;            local_min=next;
1953            local_lowi=i;            local_lowi=i;
# Line 1918  Data::calc_minGlobalDataPoint(int& ProcN Line 1965  Data::calc_minGlobalDataPoint(int& ProcN
1965    
1966  #ifdef PASO_MPI  #ifdef PASO_MPI
1967      // determine the processor on which the minimum occurs      // determine the processor on which the minimum occurs
1968      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPointRO(lowi,lowj);
1969      int lowProc = 0;      int lowProc = 0;
1970      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1971      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1972        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1973    
1974      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1975          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1977  Data::saveVTK(std::string fileName) cons Line 2025  Data::saveVTK(std::string fileName) cons
2025    }    }
2026    boost::python::dict args;    boost::python::dict args;
2027    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2028    getDomain()->saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
2029    return;    return;
2030  }  }
2031    
# Line 1987  Data::operator+=(const Data& right) Line 2035  Data::operator+=(const Data& right)
2035    if (isProtected()) {    if (isProtected()) {
2036          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2037    }    }
2038    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2039    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
2040      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
2041          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
2042  }  }
2043    
2044  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2048  Data::operator+=(const boost::python::ob
2048          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2049    }    }
2050    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2051    if (isLazy())    (*this)+=tmp;
2052    {    return *this;
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD);   // for lazy + is equivalent to +=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,plus<double>());  
     return (*this);  
   }  
2053  }  }
2054    
2055  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
2056  Data&  Data&
2057  Data::operator=(const Data& other)  Data::operator=(const Data& other)
2058  {  {
2059    #if defined ASSIGNMENT_MEANS_DEEPCOPY  
2060    // This should not be used.
2061    // Just leaving this here until I have completed transition
2062    copy(other);    copy(other);
2063    #else
2064      m_protected=false;        // since any changes should be caught by exclusiveWrite();
2065    //   m_data=other.m_data;
2066      set_m_data(other.m_data);
2067    #endif
2068    return (*this);    return (*this);
2069  }  }
2070    
# Line 2034  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,SUB)
2078    {    exclusiveWrite();
2079      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
2080          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
2081  }  }
2082    
2083  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2087  Data::operator-=(const boost::python::ob
2087          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2088    }    }
2089    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2090    if (isLazy())    (*this)-=tmp;
2091    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB);   // for lazy - is equivalent to -=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,minus<double>());  
     return (*this);  
   }  
2092  }  }
2093    
2094  Data&  Data&
# Line 2073  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,MUL)
2101    {    exclusiveWrite();
2102      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
2103          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
2104  }  }
2105    
2106  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2110  Data::operator*=(const boost::python::ob
2110          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2111    }    }
2112    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2113    if (isLazy())    (*this)*=tmp;
2114    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL);   // for lazy * is equivalent to *=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,multiplies<double>());  
     return (*this);  
   }  
2115  }  }
2116    
2117  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2120  Data::operator/=(const Data& right)
2120    if (isProtected()) {    if (isProtected()) {
2121          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2122    }    }
2123    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2124    {    exclusiveWrite();
2125      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
2126          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
2127  }  }
2128    
2129  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2133  Data::operator/=(const boost::python::ob
2133          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2134    }    }
2135    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2136    if (isLazy())    (*this)/=tmp;
2137    {    return (*this);
     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV);   // for lazy / is equivalent to /=  
         m_data=c->getPtr();  
     return (*this);  
   }  
   else  
   {  
     binaryOp(tmp,divides<double>());  
     return (*this);  
   }  
2138  }  }
2139    
2140  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2154  Data::powO(const boost::python::object&
2154  Data  Data
2155  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2156  {  {
2157    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2158    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2159  }  }
2160    
# Line 2175  Data::powD(const Data& right) const Line 2163  Data::powD(const Data& right) const
2163  Data  Data
2164  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2165  {  {
2166    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2167    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2168  }  }
2169    
# Line 2188  escript::operator+(const Data& left, con Line 2172  escript::operator+(const Data& left, con
2172  Data  Data
2173  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2174  {  {
2175    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2176    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2177  }  }
2178    
# Line 2201  escript::operator-(const Data& left, con Line 2181  escript::operator-(const Data& left, con
2181  Data  Data
2182  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2183  {  {
2184    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2185    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2186  }  }
2187    
# Line 2214  escript::operator*(const Data& left, con Line 2190  escript::operator*(const Data& left, con
2190  Data  Data
2191  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2192  {  {
2193    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2194    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2195  }  }
2196    
# Line 2227  escript::operator/(const Data& left, con Line 2199  escript::operator/(const Data& left, con
2199  Data  Data
2200  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2201  {  {
2202    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2203    {    MAKELAZYBIN2(left,tmp,ADD)
2204      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);  
2205  }  }
2206    
2207  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2209  escript::operator+(const Data& left, con
2209  Data  Data
2210  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2211  {  {
2212    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2213    {    MAKELAZYBIN2(left,tmp,SUB)
2214      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);  
2215  }  }
2216    
2217  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2219  escript::operator-(const Data& left, con
2219  Data  Data
2220  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2221  {  {
2222    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2223    {    MAKELAZYBIN2(left,tmp,MUL)
2224      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);  
2225  }  }
2226    
2227  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2229  escript::operator*(const Data& left, con
2229  Data  Data
2230  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2231  {  {
2232    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2233    {    MAKELAZYBIN2(left,tmp,DIV)
2234      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);  
2235  }  }
2236    
2237  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2239  escript::operator/(const Data& left, con
2239  Data  Data
2240  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2241  {  {
2242    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2243    {    MAKELAZYBIN2(tmp,right,ADD)
2244      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;  
2245  }  }
2246    
2247  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2249  escript::operator+(const boost::python::
2249  Data  Data
2250  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2251  {  {
2252    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2253    {    MAKELAZYBIN2(tmp,right,SUB)
2254      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;  
2255  }  }
2256    
2257  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2259  escript::operator-(const boost::python::
2259  Data  Data
2260  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2261  {  {
2262    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2263    {    MAKELAZYBIN2(tmp,right,MUL)
2264      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;  
2265  }  }
2266    
2267  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2269  escript::operator*(const boost::python::
2269  Data  Data
2270  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2271  {  {
2272    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2273    {    MAKELAZYBIN2(tmp,right,DIV)
2274      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;  
2275  }  }
2276    
2277    
# Line 2364  void Line 2312  void
2312  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2313                 const Data& value)                 const Data& value)
2314  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2315    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2316    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2317      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2318    }    }
2319      exclusiveWrite();
2320    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2321       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2322    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2331  Data::setSlice(const Data& value,
2331    if (isProtected()) {    if (isProtected()) {
2332          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2333    }    }
2334    FORCERESOLVE;    forceResolve();
2335  /*  if (isLazy())    exclusiveWrite();     // In case someone finds a way to call this without going through setItemD
   {  
     throw DataException("Error - setSlice not permitted on lazy data.");  
   }*/  
2336    Data tempValue(value);    Data tempValue(value);
2337    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2338    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2431  Data::typeMatchRight(const Data& right) Line 2375  Data::typeMatchRight(const Data& right)
2375    }    }
2376  }  }
2377    
2378    // The normal TaggedValue adds the tag if it is not already present
2379    // This form does not. It throws instead.
2380    // This is because the names are maintained by the domain and cannot be added
2381    // without knowing the tag number to map it to.
2382  void  void
2383  Data::setTaggedValueByName(std::string name,  Data::setTaggedValueByName(std::string name,
2384                             const boost::python::object& value)                             const boost::python::object& value)
2385  {  {
2386       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2387      FORCERESOLVE;      forceResolve();
2388        exclusiveWrite();
2389          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2390          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2391       }       }
2392         else
2393         {                  // The
2394        throw DataException("Error - unknown tag in setTaggedValueByName.");
2395         }
2396  }  }
2397    
2398  void  void
2399  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2400                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2404  Data::setTaggedValue(int tagKey,
2404    }    }
2405    //    //
2406    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2407    FORCERESOLVE;    forceResolve();
2408      exclusiveWrite();
2409    if (isConstant()) tag();    if (isConstant()) tag();
2410    numeric::array asNumArray(value);    WrappedArray w(value);
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape;  
   for (int i=0; i < asNumArray.getrank(); i++) {  
     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
   }  
2411    
2412    DataVector temp_data2;    DataVector temp_data2;
2413    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2414    
2415    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2416  }  }
2417    
2418    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2427  Data::setTaggedValueFromCPP(int tagKey,
2427    }    }
2428    //    //
2429    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2430    FORCERESOLVE;    forceResolve();
2431    if (isConstant()) tag();    if (isConstant()) tag();
2432      exclusiveWrite();
2433    //    //
2434    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2435    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2462  escript::C_GeneralTensorProduct(Data& ar
2462    // 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().
2463    
2464    // deal with any lazy data    // deal with any lazy data
2465    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2466    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2467      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2468      {
2469        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2470        return Data(c);
2471      }
2472    
2473    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2474    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2545  escript::C_GeneralTensorProduct(Data& ar
2545       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
2546    }    }
2547    
2548      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2549      {
2550         ostringstream os;
2551         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2552         throw DataException(os.str());
2553      }
2554    
2555    // Declare output Data object    // Declare output Data object
2556    Data res;    Data res;
2557    
2558    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2559      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2560      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2561      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2562      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2563      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2564    }    }
2565    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
# Line 2618  escript::C_GeneralTensorProduct(Data& ar Line 2580  escript::C_GeneralTensorProduct(Data& ar
2580    
2581      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2582      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2583      double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
     // Get the views  
 //     DataArrayView view_1 = tmp_1->getDefaultValue();  
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_1 = &((view_1.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_1 = &(tmp_1->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
2584    
2585        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2586        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2587    
2588      // Compute an MVP for the default      // Compute an MVP for the default
2589      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2637  escript::C_GeneralTensorProduct(Data& ar Line 2592  escript::C_GeneralTensorProduct(Data& ar
2592      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2593      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2594        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
 //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);  
 //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);  
 //       double *ptr_1 = &view_1.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
2595    
2596        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2597        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2598            
2599        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2600      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2618  escript::C_GeneralTensorProduct(Data& ar
2618        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2619          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2620          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2621          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2622          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2623          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2624          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2625        }        }
2626      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2644  escript::C_GeneralTensorProduct(Data& ar
2644    
2645      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2646      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2647      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2648      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2649  //     DataArrayView view_0 = tmp_0->getDefaultValue();      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_0 = &((view_0.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_0 = &(tmp_0->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
2650    
2651      // Compute an MVP for the default      // Compute an MVP for the default
2652      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2710  escript::C_GeneralTensorProduct(Data& ar Line 2654  escript::C_GeneralTensorProduct(Data& ar
2654      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2655      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2656      for (i=lookup_0.begin();i!=lookup_0.end();i++) {      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
 //      tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());  
 //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);  
 //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);  
 //       double *ptr_0 = &view_0.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
2657    
2658        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2659        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2660        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2661        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2662      }      }
2663    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2678  escript::C_GeneralTensorProduct(Data& ar
2678      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2679      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2680    
2681  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2682  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2683  //     DataArrayView view_1 = tmp_1->getDefaultValue();      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
 //     DataArrayView view_2 = tmp_2->getDefaultValue();  
 //     // Get the pointers to the actual data  
 //     double *ptr_0 = &((view_0.getData())[0]);  
 //     double *ptr_1 = &((view_1.getData())[0]);  
 //     double *ptr_2 = &((view_2.getData())[0]);  
   
     double *ptr_0 = &(tmp_0->getDefaultValue(0));  
     double *ptr_1 = &(tmp_1->getDefaultValue(0));  
     double *ptr_2 = &(tmp_2->getDefaultValue(0));  
   
2684    
2685      // Compute an MVP for the default      // Compute an MVP for the default
2686      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
# Line 2768  escript::C_GeneralTensorProduct(Data& ar Line 2697  escript::C_GeneralTensorProduct(Data& ar
2697      // Compute an MVP for each tag      // Compute an MVP for each tag
2698      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2699      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2700  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2701  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2702  //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
 //       double *ptr_0 = &view_0.getData(0);  
 //       double *ptr_1 = &view_1.getData(0);  
 //       double *ptr_2 = &view_2.getData(0);  
   
       double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));  
       double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));  
       double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));  
2703    
2704        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2705      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2721  escript::C_GeneralTensorProduct(Data& ar
2721      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2722      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2723        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2724        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2725        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2726          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2727          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2728          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2729          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2730          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2731        }        }
2732      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2750  escript::C_GeneralTensorProduct(Data& ar
2750        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2751          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2752          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2753          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2754          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2755          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2756          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2757        }        }
2758      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2775  escript::C_GeneralTensorProduct(Data& ar
2775      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2776      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2777        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2778        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2779        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2780          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2781          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2782          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2783          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2784          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2785        }        }
2786      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2805  escript::C_GeneralTensorProduct(Data& ar
2805          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2806          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2807          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2808          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2809          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2810          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2811          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2812        }        }
2813      }      }
# Line 2923  Data::borrowReadyPtr() const Line 2845  Data::borrowReadyPtr() const
2845  std::string  std::string
2846  Data::toString() const  Data::toString() const
2847  {  {
2848      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2849      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2850        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2851      {      {
2852      stringstream temp;      stringstream temp;
2853      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
# Line 2934  Data::toString() const Line 2857  Data::toString() const
2857  }  }
2858    
2859    
2860    // This method is not thread-safe
2861    DataTypes::ValueType::reference
2862    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
2863    {
2864        checkExclusiveWrite();
2865        return getReady()->getDataAtOffsetRW(i);
2866    }
2867    
2868    // This method is not thread-safe
2869  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2870  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
2871  {  {
2872      if (isLazy())      forceResolve();
2873      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
2874  }  }
2875    
2876    
2877  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
2878  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
2879  {  // {
2880  //     if (isLazy())  //     if (isLazy())
2881  //     {  //     {
2882  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
2883  //     }  //     }
2884      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
2885      return getReady()->getDataAtOffset(i);  // }
2886  }  
2887    
2888  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2889  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
2890  {  {
2891      forceResolve();
2892    if (!isReady())    if (!isReady())
2893    {    {
2894      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data (object is const which prevents resolving).");      throw DataException("Programmer error -getDataPointRO() not permitted on Lazy Data.");
2895    }    }
2896    else    else
2897    {    {
2898      const DataReady* dr=getReady();      const DataReady* dr=getReady();
2899      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
2900    }    }
2901  }  }
2902    
2903    
2904    
2905    
2906  DataTypes::ValueType::reference  DataTypes::ValueType::reference
2907  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
2908  {  {
2909    FORCERESOLVE;    checkExclusiveWrite();
2910    if (!isReady())    DataReady* dr=getReady();
2911    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
2912      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
2913    }  
2914    else  BufferGroup*
2915    {  Data::allocSampleBuffer() const
2916      DataReady* dr=getReady();  {
2917      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
2918    }       {
2919        #ifdef _OPENMP
2920        int tnum=omp_get_max_threads();
2921        #else
2922        int tnum=1;
2923        #endif
2924        return new BufferGroup(getSampleBufferSize(),tnum);
2925         }
2926         else
2927         {
2928        return NULL;
2929         }
2930    }
2931    
2932    void
2933    Data::freeSampleBuffer(BufferGroup* bufferg)
2934    {
2935         if (bufferg!=0)
2936         {
2937        delete bufferg;
2938         }
2939  }  }
2940    
2941    
# Line 3000  Data::print() Line 2951  Data::print()
2951    {    {
2952      printf( "[%6d]", i );      printf( "[%6d]", i );
2953      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
2954        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
2955      printf( "\n" );      printf( "\n" );
2956    }    }
2957  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 2971  Data::dump(const std::string fileName) c
2971            return m_data->dump(fileName);            return m_data->dump(fileName);
2972      }      }
2973    }    }
2974    catch (exception& e)    catch (std::exception& e)
2975    {    {
2976          cout << e.what() << endl;          cout << e.what() << endl;
2977    }    }

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

  ViewVC Help
Powered by ViewVC 1.1.26