/[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 2458 by jfenwick, Wed Jun 3 06:18:21 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    // This should be safer once the DataC RO changes have been brought in
92    template <class ARR>
93    boost::python::tuple
94    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
95    {
96       using namespace boost::python;
97       using boost::python::list;
98       int rank=shape.size();
99       if (rank==0)
100       {
101        return make_tuple(v[0]);
102       }
103       else if (rank==1)
104       {
105        list l;
106        for (size_t i=0;i<shape[0];++i)
107        {
108           l.append(v[i]);
109        }
110        return tuple(l);
111       }
112       else if (rank==2)
113       {
114        list lj;
115        for (size_t j=0;j<shape[0];++j)
116        {
117           list li;
118           for (size_t i=0;i<shape[1];++i)
119           {
120              li.append(v[DataTypes::getRelIndex(shape,j,i)]);
121           }
122           lj.append(tuple(li));
123        }
124        return tuple(lj);
125       }
126       else if (rank==3)
127       {
128        list lk;
129        for (size_t k=0;k<shape[0];++k)
130        {
131           list lj;
132           for (size_t j=0;j<shape[1];++j)
133           {
134            list li;
135            for (size_t i=0;i<shape[2];++i)
136            {
137               li.append(v[DataTypes::getRelIndex(shape,k,j,i)]);
138            }
139            lj.append(tuple(li));
140           }
141           lk.append(tuple(lj));
142        }
143        return tuple(lk);
144       }
145       else if (rank==4)
146       {
147        list ll;
148        for (size_t l=0;l<shape[0];++l)
149        {
150           list lk;
151           for (size_t k=0;k<shape[1];++k)
152           {
153            list lj;
154            for (size_t j=0;j<shape[2];++j)
155            {
156                list li;
157                for (size_t i=0;i<shape[3];++i)
158                {
159                   li.append(v[DataTypes::getRelIndex(shape,l,k,j,i)]);
160                }
161                lj.append(tuple(li));
162            }
163            lk.append(tuple(lj));
164           }
165           ll.append(tuple(lk));
166        }
167        return tuple(ll);
168       }
169       else
170         throw DataException("Unknown rank in pointToTuple.");
171    }
172    
173    }  // anonymous namespace
174    
175  Data::Data()  Data::Data()
176        : m_shared(false), m_lazy(false)
177  {  {
178    //    //
179    // Default data is type DataEmpty    // Default data is type DataEmpty
180    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
181    m_data=temp->getPtr();  //   m_data=temp->getPtr();
182      set_m_data(temp->getPtr());
183    m_protected=false;    m_protected=false;
184  }  }
185    
# Line 60  Data::Data(double value, Line 187  Data::Data(double value,
187             const tuple& shape,             const tuple& shape,
188             const FunctionSpace& what,             const FunctionSpace& what,
189             bool expanded)             bool expanded)
190        : m_shared(false), m_lazy(false)
191  {  {
192    DataTypes::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
193    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 204  Data::Data(double value,
204         const DataTypes::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
205         const FunctionSpace& what,         const FunctionSpace& what,
206             bool expanded)             bool expanded)
207        : m_shared(false), m_lazy(false)
208  {  {
209    int len = DataTypes::noValues(dataPointShape);    int len = DataTypes::noValues(dataPointShape);
210    
# Line 89  Data::Data(double value, Line 218  Data::Data(double value,
218  }  }
219    
220  Data::Data(const Data& inData)  Data::Data(const Data& inData)
221        : m_shared(false), m_lazy(false)
222  {  {
223    m_data=inData.m_data;  //   m_data=inData.m_data;
224      set_m_data(inData.m_data);
225    m_protected=inData.isProtected();    m_protected=inData.isProtected();
226  }  }
227    
228    
229  Data::Data(const Data& inData,  Data::Data(const Data& inData,
230             const DataTypes::RegionType& region)             const DataTypes::RegionType& region)
231        : m_shared(false), m_lazy(false)
232  {  {
233    DataAbstract_ptr dat=inData.m_data;    DataAbstract_ptr dat=inData.m_data;
234    if (inData.isLazy())    if (inData.isLazy())
# Line 110  Data::Data(const Data& inData, Line 242  Data::Data(const Data& inData,
242    //    //
243    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
244    DataAbstract* tmp = dat->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
245    m_data=DataAbstract_ptr(tmp);  //   m_data=DataAbstract_ptr(tmp);
246      set_m_data(DataAbstract_ptr(tmp));
247    m_protected=false;    m_protected=false;
248    
249  }  }
250    
251  Data::Data(const Data& inData,  Data::Data(const Data& inData,
252             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
253        : m_shared(false), m_lazy(false)
254  {  {
255    if (inData.isEmpty())    if (inData.isEmpty())
256    {    {
257      throw DataException("Error - will not interpolate for instances of DataEmpty.");      throw DataException("Error - will not interpolate for instances of DataEmpty.");
258    }    }
259    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
260      m_data=inData.m_data;  //     m_data=inData.m_data;
261        set_m_data(inData.m_data);
262    }    }
263    else    else
264    {    {
# Line 130  Data::Data(const Data& inData, Line 266  Data::Data(const Data& inData,
266      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
267        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
268        {           // 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
269      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
270        }        }
271        // 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
272        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
273        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());          DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
274        m_data=DataAbstract_ptr(dc);  //       m_data=DataAbstract_ptr(dc);
275          set_m_data(DataAbstract_ptr(dc));
276      } else {      } else {
277        Data tmp(0,inData.getDataPointShape(),functionspace,true);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
278        // 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 286  Data::Data(const Data& inData,
286        } else {        } else {
287          inDataDomain->interpolateACross(tmp,inData);          inDataDomain->interpolateACross(tmp,inData);
288        }        }
289        m_data=tmp.m_data;  //       m_data=tmp.m_data;
290          set_m_data(tmp.m_data);
291      }      }
292    }    }
293    m_protected=false;    m_protected=false;
294  }  }
295    
296  Data::Data(DataAbstract* underlyingdata)  Data::Data(DataAbstract* underlyingdata)
297        : m_shared(false), m_lazy(false)
298  {  {
299  //  m_data=shared_ptr<DataAbstract>(underlyingdata);      set_m_data(underlyingdata->getPtr());
     m_data=underlyingdata->getPtr();  
300      m_protected=false;      m_protected=false;
301  }  }
302    
303  Data::Data(DataAbstract_ptr underlyingdata)  Data::Data(DataAbstract_ptr underlyingdata)
304        : m_shared(false), m_lazy(false)
305  {  {
306      m_data=underlyingdata;      set_m_data(underlyingdata);
307      m_protected=false;      m_protected=false;
308  }  }
309    
   
 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;  
 }*/  
   
310  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
311           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
312                   const FunctionSpace& what,                   const FunctionSpace& what,
313                   bool expanded)                   bool expanded)
314        : m_shared(false), m_lazy(false)
315  {  {
316     initialise(value,shape,what,expanded);     initialise(value,shape,what,expanded);
317     m_protected=false;     m_protected=false;
# Line 198  Data::Data(const DataTypes::ValueType& v Line 321  Data::Data(const DataTypes::ValueType& v
321  Data::Data(const object& value,  Data::Data(const object& value,
322         const FunctionSpace& what,         const FunctionSpace& what,
323             bool expanded)             bool expanded)
324        : m_shared(false), m_lazy(false)
325  {  {
326    numeric::array asNumArray(value);    WrappedArray w(value);
327    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
328    m_protected=false;    m_protected=false;
329  }  }
330    
331    
332  Data::Data(const object& value,  Data::Data(const object& value,
333             const Data& other)             const Data& other)
334        : m_shared(false), m_lazy(false)
335  {  {
336    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);  
337    
338    //    // extract the shape of the array
339    // Create DataConstant using the given value and all other parameters    const DataTypes::ShapeType& tempShape=w.getShape();
340    // copied from other. If value is a rank 0 object this Data    if (w.getRank()==0) {
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
341    
342    
343      // get the space for the data vector      // get the space for the data vector
344      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
345      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
346      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
347    
348      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
349    
350      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);  
   
351      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
352  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
353  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
354    
355    } else {    } else {
356      //      //
357      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
358  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
359      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
360  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
361    }    }
362    m_protected=false;    m_protected=false;
363  }  }
364    
365  Data::~Data()  Data::~Data()
366  {  {
367      set_m_data(DataAbstract_ptr());
368  }  }
369    
370    
371    // only call in thread safe contexts.
372    // This method should be atomic
373    void Data::set_m_data(DataAbstract_ptr p)
374    {
375      if (m_data.get()!=0)  // release old ownership
376      {
377        m_data->removeOwner(this);
378      }
379      if (p.get()!=0)
380      {
381        m_data=p;
382        m_data->addOwner(this);
383        m_shared=m_data->isShared();
384        m_lazy=m_data->isLazy();
385      }
386    }
387    
388  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
389                   const FunctionSpace& what,                   const FunctionSpace& what,
390                   bool expanded)                   bool expanded)
391  {  {
# Line 277  Data::initialise(const boost::python::nu Line 396  Data::initialise(const boost::python::nu
396    // within the shared_ptr constructor.    // within the shared_ptr constructor.
397    if (expanded) {    if (expanded) {
398      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
399  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
400  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
401    } else {    } else {
402      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
403  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
404  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
405    }    }
406  }  }
407    
# Line 302  Data::initialise(const DataTypes::ValueT Line 419  Data::initialise(const DataTypes::ValueT
419    // within the shared_ptr constructor.    // within the shared_ptr constructor.
420    if (expanded) {    if (expanded) {
421      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
422  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
423  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
424    } else {    } else {
425      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
426  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
427  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
428    }    }
429  }  }
430    
431    
 // 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";  
 //  }  
 // }  
   
432  escriptDataC  escriptDataC
433  Data::getDataC()  Data::getDataC()
434  {  {
# Line 410  Data::getDataC() const Line 445  Data::getDataC() const
445    return temp;    return temp;
446  }  }
447    
448    size_t
449    Data::getSampleBufferSize() const
450    {
451      return m_data->getSampleBufferSize();
452    }
453    
454    
455  const boost::python::tuple  const boost::python::tuple
456  Data::getShapeTuple() const  Data::getShapeTuple() const
457  {  {
# Line 435  Data::getShapeTuple() const Line 477  Data::getShapeTuple() const
477  // 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.
478  // 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
479  // 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
480  Data*  Data
481  Data::copySelf()  Data::copySelf()
482  {  {
483     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
484     return new Data(temp);     return Data(temp);
485  }  }
486    
487  void  void
# Line 447  Data::copy(const Data& other) Line 489  Data::copy(const Data& other)
489  {  {
490    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
491    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
492    m_data=p;  //   m_data=p;
493      set_m_data(p);
494  }  }
495    
496    
# Line 463  Data::delaySelf() Line 506  Data::delaySelf()
506  {  {
507    if (!isLazy())    if (!isLazy())
508    {    {
509      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
510        set_m_data((new DataLazy(m_data))->getPtr());
511    }    }
512  }  }
513    
514    
515    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
516    // to zero, all the tags from all the DataTags would be in the result.
517    // However since they all have the same value (0) whether they are there or not should not matter.
518    // So I have decided that for all types this method will create a constant 0.
519    // It can be promoted up as required.
520    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
521    // but we can deal with that if it arises.
522    //
523  void  void
524  Data::setToZero()  Data::setToZero()
525  {  {
# Line 474  Data::setToZero() Line 527  Data::setToZero()
527    {    {
528       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
529    }    }
530    m_data->setToZero();    if (isLazy())
531      {
532         DataTypes::ValueType v(getNoValues(),0);
533         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
534         DataLazy* dl=new DataLazy(dc->getPtr());
535         set_m_data(dl->getPtr());
536      }
537      else
538      {
539         exclusiveWrite();
540         m_data->setToZero();
541      }
542  }  }
543    
544    
545  void  void
546  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
547                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 603  Data::copyWithMask(const Data& other,
603    {    {
604      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
605    }    }
606      unsigned int selfrank=getDataPointRank();
607      unsigned int otherrank=other2.getDataPointRank();
608      unsigned int maskrank=mask2.getDataPointRank();
609      if ((selfrank==0) && (otherrank>0 || maskrank>0))
610      {
611        // to get here we must be copying from a large object into a scalar
612        // I am not allowing this.
613        // If you are calling copyWithMask then you are considering keeping some existing values
614        // and so I'm going to assume that you don't want your data objects getting a new shape.
615        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
616      }
617      exclusiveWrite();
618    // Now we iterate over the elements    // Now we iterate over the elements
619    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
620    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
621    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
622    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
623      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
624    {    {
625      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // Not allowing this combination.
626        // it is not clear what the rank of the target should be.
627        // Should it be filled with the scalar (rank stays the same);
628        // or should the target object be reshaped to be a scalar as well.
629        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
630      }
631      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
632      {
633        if (mvec[0]>0)      // copy whole object if scalar is >0
634        {
635            copy(other);
636        }
637        return;
638      }
639      if (isTagged())       // so all objects involved will also be tagged
640      {
641        // note the !
642        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
643            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
644        {
645            throw DataException("copyWithMask, shape mismatch.");
646        }
647    
648        // We need to consider the possibility that tags are missing or in the wrong order
649        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
650        // all tags which are not explicitly defined
651    
652        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
653        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
654        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
655    
656        // first, add any tags missing from other or mask
657        const DataTagged::DataMapType& olookup=optr->getTagLookup();
658            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
659        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
660        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
661        for (i=olookup.begin();i!=olookup.end();i++)
662        {
663               tptr->addTag(i->first);
664            }
665            for (i=mlookup.begin();i!=mlookup.end();i++) {
666               tptr->addTag(i->first);
667            }
668        // now we know that *this has all the required tags but they aren't guaranteed to be in
669        // the same order
670    
671        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
672        if ((selfrank==otherrank) && (otherrank==maskrank))
673        {
674            for (i=tlookup.begin();i!=tlookup.end();i++)
675            {
676                // get the target offset
677                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
678                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
679                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
680                for (int j=0;j<getDataPointSize();++j)
681                {
682                    if (mvec[j+moff]>0)
683                    {
684                        self[j+toff]=ovec[j+ooff];
685                    }
686                }
687                }
688            // now for the default value
689            for (int j=0;j<getDataPointSize();++j)
690            {
691                if (mvec[j+mptr->getDefaultOffset()]>0)
692                {
693                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
694                }
695            }
696        }
697        else    // other is a scalar
698        {
699            for (i=tlookup.begin();i!=tlookup.end();i++)
700            {
701                // get the target offset
702                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
703                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
704                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
705                for (int j=0;j<getDataPointSize();++j)
706                {
707                    if (mvec[j+moff]>0)
708                    {
709                        self[j+toff]=ovec[ooff];
710                    }
711                }
712                }
713            // now for the default value
714            for (int j=0;j<getDataPointSize();++j)
715            {
716                if (mvec[j+mptr->getDefaultOffset()]>0)
717                {
718                    self[j+tptr->getDefaultOffset()]=ovec[0];
719                }
720            }
721        }
722    
723        return;         // ugly
724      }
725      // mixed scalar and non-scalar operation
726      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
727      {
728            size_t num_points=self.size();
729        // OPENMP 3.0 allows unsigned loop vars.
730        #if defined(_OPENMP) && (_OPENMP < 200805)
731        long i;
732        #else
733        size_t i;
734        #endif  
735        size_t psize=getDataPointSize();    
736        #pragma omp parallel for private(i) schedule(static)
737        for (i=0;i<num_points;++i)
738        {
739            if (mvec[i]>0)
740            {
741                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
742            }                   // dest point
743        }
744        return;         // ugly!
745      }
746      // tagged data is already taken care of so we only need to worry about shapes
747      // special cases with scalars are already dealt with so all we need to worry about is shape
748      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
749      {
750        ostringstream oss;
751        oss <<"Error - size mismatch in arguments to copyWithMask.";
752        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
753        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
754        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
755        throw DataException(oss.str());
756    }    }
757    size_t num_points=self.size();    size_t num_points=self.size();
758    
# Line 564  Data::copyWithMask(const Data& other, Line 772  Data::copyWithMask(const Data& other,
772    }    }
773  }  }
774    
   
   
775  bool  bool
776  Data::isExpanded() const  Data::isExpanded() const
777  {  {
# Line 574  Data::isExpanded() const Line 780  Data::isExpanded() const
780  }  }
781    
782  bool  bool
783    Data::actsExpanded() const
784    {
785      return m_data->actsExpanded();
786    
787    }
788    
789    
790    bool
791  Data::isTagged() const  Data::isTagged() const
792  {  {
793    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 811  Data::isConstant() const
811  bool  bool
812  Data::isLazy() const  Data::isLazy() const
813  {  {
814    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
815  }  }
816    
817  // 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 842  Data::expand()
842    if (isConstant()) {    if (isConstant()) {
843      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
844      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
845  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
846  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
847    } else if (isTagged()) {    } else if (isTagged()) {
848      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
849      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
850  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
851  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
852    } else if (isExpanded()) {    } else if (isExpanded()) {
853      //      //
854      // do nothing      // do nothing
# Line 656  Data::tag() Line 868  Data::tag()
868    if (isConstant()) {    if (isConstant()) {
869      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
870      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
871  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
872  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
873    } else if (isTagged()) {    } else if (isTagged()) {
874      // do nothing      // do nothing
875    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 882  Data::tag()
882       {       {
883      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
884       }       }
885       m_data=res;      //      m_data=res;
886         set_m_data(res);
887       tag();       tag();
888    } else {    } else {
889      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 895  Data::resolve()
895  {  {
896    if (isLazy())    if (isLazy())
897    {    {
898       m_data=m_data->resolve();  //      m_data=m_data->resolve();
899        set_m_data(m_data->resolve());
900    }    }
901  }  }
902    
903    void
904    Data::requireWrite()
905    {
906      resolve();
907      exclusiveWrite();
908    }
909    
910  Data  Data
911  Data::oneOver() const  Data::oneOver() const
912  {  {
913    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
914    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
915  }  }
916    
917  Data  Data
918  Data::wherePositive() const  Data::wherePositive() const
919  {  {
920    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
921    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
922  }  }
923    
924  Data  Data
925  Data::whereNegative() const  Data::whereNegative() const
926  {  {
927    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
928    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
929  }  }
930    
931  Data  Data
932  Data::whereNonNegative() const  Data::whereNonNegative() const
933  {  {
934    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
935    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
936  }  }
937    
938  Data  Data
939  Data::whereNonPositive() const  Data::whereNonPositive() const
940  {  {
941    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
942    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
943  }  }
944    
945  Data  Data
946  Data::whereZero(double tol) const  Data::whereZero(double tol) const
947  {  {
948    Data dataAbs=abs();  //   Data dataAbs=abs();
949    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
950       MAKELAZYOPOFF(EZ,tol)
951       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
952    
953  }  }
954    
955  Data  Data
956  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
957  {  {
958    Data dataAbs=abs();  //   Data dataAbs=abs();
959    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
960      MAKELAZYOPOFF(NEZ,tol)
961      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
962    
963  }  }
964    
965  Data  Data
# Line 767  bool Line 972  bool
972  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
973  {  {
974    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());  
 //     }  
 //   }  
975  }  }
976    
977  Data  Data
# Line 819  Data::getLength() const Line 1014  Data::getLength() const
1014    return m_data->getLength();    return m_data->getLength();
1015  }  }
1016    
1017  const  const boost::python::object
1018  boost::python::numeric::array  Data::getValueOfDataPointAsTuple(int dataPointNo)
 Data:: getValueOfDataPoint(int dataPointNo)  
1019  {  {
1020    int i, j, k, l;     forceResolve();
1021       if (getNumDataPointsPerSample()>0) {
   FORCERESOLVE;  
   
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
   
   //  
   // create the numeric array to be returned  
   boost::python::numeric::array numArray(0.0);  
   
   //  
   // 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;  
   
   //  
   // resize the numeric array to the shape just calculated  
   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]);  
   }  
   
   if (getNumDataPointsPerSample()>0) {  
1022         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1023         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1024         //         //
1025         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1026         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1027             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1028         }         }
1029    
1030         //         //
1031         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1032         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1033             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1034         }         }
1035         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1036         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1037           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1038      }
1039         switch( dataPointRank ){    else
1040              case 0 :    {
1041                  numArray[0] = getDataAtOffset(offset);      // The old method would return an empty numArray of the given shape
1042                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1043              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;  
     }  
1044    }    }
   //  
   // return the array  
   return numArray;  
1045    
1046  }  }
1047    
1048    
1049  void  void
1050  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1051  {  {
1052      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1053      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1054  }  }
1055    
1056  void  void
1057  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1058  {  {
1059    if (isProtected()) {    if (isProtected()) {
1060          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1061    }    }
1062    FORCERESOLVE;    forceResolve();
1063    
1064      WrappedArray w(obj);
1065    //    //
1066    // check rank    // check rank
1067    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1068        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1069    
1070    //    //
1071    // check shape of num_array    // check shape of array
1072    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1073      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1074         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1075    }    }
1076    //    //
1077    // make sure data is expanded:    // make sure data is expanded:
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1082  Data::setValueOfDataPointToArray(int dat
1082    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1083         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1084         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1085         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1086    } else {    } else {
1087         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1088    }    }
1089  }  }
1090    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1096  Data::setValueOfDataPoint(int dataPointN
1096    }    }
1097    //    //
1098    // make sure data is expanded:    // make sure data is expanded:
1099    FORCERESOLVE;    forceResolve();
1100    if (!isExpanded()) {    if (!isExpanded()) {
1101      expand();      expand();
1102    }    }
# Line 977  Data::setValueOfDataPoint(int dataPointN Line 1110  Data::setValueOfDataPoint(int dataPointN
1110  }  }
1111    
1112  const  const
1113  boost::python::numeric::array  boost::python::object
1114  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1115  {  {
1116    size_t length=0;    // This could be lazier than it is now
1117    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();  
1118    
1119    //    // copy datapoint into a buffer
1120    // create the numeric array to be returned    // broadcast buffer to all nodes
1121    boost::python::numeric::array numArray(0.0);    // convert buffer to tuple
1122      // 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;  
1123    
1124    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1125    // 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]);  
   }  
1126    
1127    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1128    double *tmpData = new double[length];    double *tmpData = new double[length];
1129    
1130    //    // updated for the MPI case
1131    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1132          if (getNumDataPointsPerSample()>0) {
1133      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1134      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1135               if (getNumDataPointsPerSample()>0) {      //
1136                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1137                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1138                  //          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;  
         }  
             }  
1139      }      }
1140          #ifdef PASO_MPI  
1141          // broadcast the data to all other processes      //
1142      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1143          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1144            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;  
1145      }      }
1146        // TODO: global error handling
1147        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1148    
1149        memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1150         }
1151      }
1152    #ifdef PASO_MPI
1153      // broadcast the data to all other processes
1154      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1155    #endif
1156    
1157      delete [] tmpData;    boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1158      delete [] tmpData;
1159    //    //
1160    // return the loaded array    // return the loaded array
1161    return numArray;    return t;
1162    
1163  }  }
1164    
1165    
1166  boost::python::numeric::array  boost::python::object
1167  Data::integrate_const() const  Data::integrateToTuple_const() const
1168  {  {
1169    if (isLazy())    if (isLazy())
1170    {    {
# Line 1127  Data::integrate_const() const Line 1173  Data::integrate_const() const
1173    return integrateWorker();    return integrateWorker();
1174  }  }
1175    
1176  boost::python::numeric::array  boost::python::object
1177  Data::integrate()  Data::integrateToTuple()
1178  {  {
1179    if (isLazy())    if (isLazy())
1180    {    {
1181      expand();      expand();
1182    }    }
1183    return integrateWorker();    return integrateWorker();
 }  
   
1184    
1185    }
1186    
1187  boost::python::numeric::array  boost::python::object
1188  Data::integrateWorker() const  Data::integrateWorker() const
1189  {  {
   int index;  
   int rank = getDataPointRank();  
1190    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1191    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1192    
# Line 1151  Data::integrateWorker() const Line 1194  Data::integrateWorker() const
1194    // calculate the integral values    // calculate the integral values
1195    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1196    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1197      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1198      if (dom==0)
1199      {            
1200        throw DataException("Can not integrate over non-continuous domains.");
1201      }
1202  #ifdef PASO_MPI  #ifdef PASO_MPI
1203    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1204    // 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
1205    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1206    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1207    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]; }
1208    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 );
1209    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1210      tuple result=pointToTuple(shape,tmp);
1211    delete[] tmp;    delete[] tmp;
1212    delete[] tmp_local;    delete[] tmp_local;
1213  #else  #else
1214    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1215    /*  double *tmp = new double[dataPointSize];
1216      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1217      tuple result=pointToTuple(shape,integrals);
1218    //   delete tmp;
1219  #endif  #endif
1220    
   //  
   // 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];  
           }  
         }  
       }  
     }  
   }  
1221    
1222    //    return result;
   // return the loaded array  
   return bp_array;  
1223  }  }
1224    
1225  Data  Data
1226  Data::sin() const  Data::sin() const
1227  {  {
1228    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1229    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1230  }  }
1231    
1232  Data  Data
1233  Data::cos() const  Data::cos() const
1234  {  {
1235    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1236    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1237  }  }
1238    
1239  Data  Data
1240  Data::tan() const  Data::tan() const
1241  {  {
1242    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1243    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1244  }  }
1245    
1246  Data  Data
1247  Data::asin() const  Data::asin() const
1248  {  {
1249    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1250    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1251  }  }
1252    
1253  Data  Data
1254  Data::acos() const  Data::acos() const
1255  {  {
1256    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1257    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1258  }  }
1259    
# Line 1279  Data::acos() const Line 1261  Data::acos() const
1261  Data  Data
1262  Data::atan() const  Data::atan() const
1263  {  {
1264    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1265    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1266  }  }
1267    
1268  Data  Data
1269  Data::sinh() const  Data::sinh() const
1270  {  {
1271    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1272    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1273  }  }
1274    
1275  Data  Data
1276  Data::cosh() const  Data::cosh() const
1277  {  {
1278    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1279    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1280  }  }
1281    
1282  Data  Data
1283  Data::tanh() const  Data::tanh() const
1284  {  {
1285    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1286    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1287  }  }
1288    
# Line 1327  Data::erf() const Line 1293  Data::erf() const
1293  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1294    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1295  #else  #else
1296    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1297    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1298  #endif  #endif
1299  }  }
# Line 1339  Data::erf() const Line 1301  Data::erf() const
1301  Data  Data
1302  Data::asinh() const  Data::asinh() const
1303  {  {
1304    if (isLazy())    MAKELAZYOP(ASINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
1305  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1306    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1307  #else  #else
# Line 1354  Data::asinh() const Line 1312  Data::asinh() const
1312  Data  Data
1313  Data::acosh() const  Data::acosh() const
1314  {  {
1315    if (isLazy())    MAKELAZYOP(ACOSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
1316  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1317    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1318  #else  #else
# Line 1369  Data::acosh() const Line 1323  Data::acosh() const
1323  Data  Data
1324  Data::atanh() const  Data::atanh() const
1325  {  {
1326    if (isLazy())    MAKELAZYOP(ATANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
1327  #if defined (_WIN32) && !defined(__INTEL_COMPILER)  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1328    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1329  #else  #else
# Line 1383  Data::atanh() const Line 1333  Data::atanh() const
1333    
1334  Data  Data
1335  Data::log10() const  Data::log10() const
1336  {  if (isLazy())  {
1337    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1338    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1339  }  }
1340    
1341  Data  Data
1342  Data::log() const  Data::log() const
1343  {  {
1344    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1345    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1346  }  }
1347    
1348  Data  Data
1349  Data::sign() const  Data::sign() const
1350  {  {
1351    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1352    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1353  }  }
1354    
1355  Data  Data
1356  Data::abs() const  Data::abs() const
1357  {  {
1358    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1359    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1360  }  }
1361    
1362  Data  Data
1363  Data::neg() const  Data::neg() const
1364  {  {
1365    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1366    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1367  }  }
1368    
# Line 1448  Data::pos() const Line 1379  Data::pos() const
1379    
1380  Data  Data
1381  Data::exp() const  Data::exp() const
1382  {    {
1383    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1384    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1385  }  }
1386    
1387  Data  Data
1388  Data::sqrt() const  Data::sqrt() const
1389  {  {
1390    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1391    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1392  }  }
1393    
# Line 1483  Data::Lsup() Line 1406  Data::Lsup()
1406  {  {
1407     if (isLazy())     if (isLazy())
1408     {     {
1409      expand();      resolve();
1410     }     }
1411     return LsupWorker();     return LsupWorker();
1412  }  }
# Line 1503  Data::sup() Line 1426  Data::sup()
1426  {  {
1427     if (isLazy())     if (isLazy())
1428     {     {
1429      expand();      resolve();
1430     }     }
1431     return supWorker();     return supWorker();
1432  }  }
# Line 1523  Data::inf() Line 1446  Data::inf()
1446  {  {
1447     if (isLazy())     if (isLazy())
1448     {     {
1449      expand();      resolve();
1450     }     }
1451     return infWorker();     return infWorker();
1452  }  }
# Line 1615  Data::minval() const Line 1538  Data::minval() const
1538  Data  Data
1539  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1540  {  {
1541         if (isLazy())
1542         {
1543        Data temp(*this);
1544        temp.resolve();
1545        return temp.swapaxes(axis0,axis1);
1546         }
1547       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1548       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1549       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1601  Data::symmetric() const
1601       else {       else {
1602          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.");
1603       }       }
1604       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1605       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1606       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1607       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1611  Data::symmetric() const
1611  Data  Data
1612  Data::nonsymmetric() const  Data::nonsymmetric() const
1613  {  {
1614       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1615       // check input       // check input
1616       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1617       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1643  Data::nonsymmetric() const
1643       }       }
1644  }  }
1645    
   
 // 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).  
1646  Data  Data
1647  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1648  {  {    
1649       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1650         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1651       {       {
1652      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);  
1653       }       }
1654       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1655       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1701  Data::trace(int axis_offset) const
1701  Data  Data
1702  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1703  {      {    
1704       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);  
      }  
1705       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1706       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1707       // 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 1803  Data::calc_minGlobalDataPoint(int& ProcN
1803    double next,local_min;    double next,local_min;
1804    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1805    
1806    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1807    {    {
1808      local_min=min;      local_min=min;
1809      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1810      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1811        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1812          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1813          if (next<local_min) {          if (next<local_min) {
1814            local_min=next;            local_min=next;
1815            local_lowi=i;            local_lowi=i;
# Line 1918  Data::calc_minGlobalDataPoint(int& ProcN Line 1827  Data::calc_minGlobalDataPoint(int& ProcN
1827    
1828  #ifdef PASO_MPI  #ifdef PASO_MPI
1829      // determine the processor on which the minimum occurs      // determine the processor on which the minimum occurs
1830      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPointRO(lowi,lowj);
1831      int lowProc = 0;      int lowProc = 0;
1832      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1833      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1834        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1835    
1836      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1837          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1977  Data::saveVTK(std::string fileName) cons Line 1887  Data::saveVTK(std::string fileName) cons
1887    }    }
1888    boost::python::dict args;    boost::python::dict args;
1889    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
1890    getDomain()->saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
1891    return;    return;
1892  }  }
1893    
# Line 1987  Data::operator+=(const Data& right) Line 1897  Data::operator+=(const Data& right)
1897    if (isProtected()) {    if (isProtected()) {
1898          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1899    }    }
1900    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
1901    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
1902      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
1903          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
1904  }  }
1905    
1906  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 1910  Data::operator+=(const boost::python::ob
1910          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1911    }    }
1912    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1913    if (isLazy())    (*this)+=tmp;
1914    {    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);  
   }  
1915  }  }
1916    
1917  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
1918  Data&  Data&
1919  Data::operator=(const Data& other)  Data::operator=(const Data& other)
1920  {  {
1921    #if defined ASSIGNMENT_MEANS_DEEPCOPY  
1922    // This should not be used.
1923    // Just leaving this here until I have completed transition
1924    copy(other);    copy(other);
1925    #else
1926      m_protected=false;        // since any changes should be caught by exclusiveWrite();
1927    //   m_data=other.m_data;
1928      set_m_data(other.m_data);
1929    #endif
1930    return (*this);    return (*this);
1931  }  }
1932    
# Line 2034  Data::operator-=(const Data& right) Line 1936  Data::operator-=(const Data& right)
1936    if (isProtected()) {    if (isProtected()) {
1937          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1938    }    }
1939    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
1940    {    exclusiveWrite();
1941      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
1942          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
1943  }  }
1944    
1945  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 1949  Data::operator-=(const boost::python::ob
1949          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1950    }    }
1951    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1952    if (isLazy())    (*this)-=tmp;
1953    {    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);  
   }  
1954  }  }
1955    
1956  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 1959  Data::operator*=(const Data& right)
1959    if (isProtected()) {    if (isProtected()) {
1960          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1961    }    }
1962    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
1963    {    exclusiveWrite();
1964      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
1965          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
1966  }  }
1967    
1968  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 1972  Data::operator*=(const boost::python::ob
1972          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1973    }    }
1974    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1975    if (isLazy())    (*this)*=tmp;
1976    {    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);  
   }  
1977  }  }
1978    
1979  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 1982  Data::operator/=(const Data& right)
1982    if (isProtected()) {    if (isProtected()) {
1983          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1984    }    }
1985    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
1986    {    exclusiveWrite();
1987      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
1988          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
1989  }  }
1990    
1991  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 1995  Data::operator/=(const boost::python::ob
1995          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1996    }    }
1997    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1998    if (isLazy())    (*this)/=tmp;
1999    {    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);  
   }  
2000  }  }
2001    
2002  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2016  Data::powO(const boost::python::object&
2016  Data  Data
2017  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2018  {  {
2019    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2020    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2021  }  }
2022    
# Line 2175  Data::powD(const Data& right) const Line 2025  Data::powD(const Data& right) const
2025  Data  Data
2026  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2027  {  {
2028    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2029    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2030  }  }
2031    
# Line 2188  escript::operator+(const Data& left, con Line 2034  escript::operator+(const Data& left, con
2034  Data  Data
2035  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2036  {  {
2037    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2038    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2039  }  }
2040    
# Line 2201  escript::operator-(const Data& left, con Line 2043  escript::operator-(const Data& left, con
2043  Data  Data
2044  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2045  {  {
2046    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2047    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2048  }  }
2049    
# Line 2214  escript::operator*(const Data& left, con Line 2052  escript::operator*(const Data& left, con
2052  Data  Data
2053  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2054  {  {
2055    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2056    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2057  }  }
2058    
# Line 2227  escript::operator/(const Data& left, con Line 2061  escript::operator/(const Data& left, con
2061  Data  Data
2062  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2063  {  {
2064    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2065    {    MAKELAZYBIN2(left,tmp,ADD)
2066      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);  
2067  }  }
2068    
2069  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2071  escript::operator+(const Data& left, con
2071  Data  Data
2072  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2073  {  {
2074    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2075    {    MAKELAZYBIN2(left,tmp,SUB)
2076      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);  
2077  }  }
2078    
2079  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2081  escript::operator-(const Data& left, con
2081  Data  Data
2082  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2083  {  {
2084    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2085    {    MAKELAZYBIN2(left,tmp,MUL)
2086      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);  
2087  }  }
2088    
2089  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2091  escript::operator*(const Data& left, con
2091  Data  Data
2092  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2093  {  {
2094    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2095    {    MAKELAZYBIN2(left,tmp,DIV)
2096      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);  
2097  }  }
2098    
2099  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2101  escript::operator/(const Data& left, con
2101  Data  Data
2102  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2103  {  {
2104    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2105    {    MAKELAZYBIN2(tmp,right,ADD)
2106      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;  
2107  }  }
2108    
2109  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2111  escript::operator+(const boost::python::
2111  Data  Data
2112  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2113  {  {
2114    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2115    {    MAKELAZYBIN2(tmp,right,SUB)
2116      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;  
2117  }  }
2118    
2119  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2121  escript::operator-(const boost::python::
2121  Data  Data
2122  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2123  {  {
2124    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2125    {    MAKELAZYBIN2(tmp,right,MUL)
2126      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;  
2127  }  }
2128    
2129  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2131  escript::operator*(const boost::python::
2131  Data  Data
2132  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2133  {  {
2134    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2135    {    MAKELAZYBIN2(tmp,right,DIV)
2136      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;  
2137  }  }
2138    
2139    
# Line 2364  void Line 2174  void
2174  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2175                 const Data& value)                 const Data& value)
2176  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2177    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2178    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2179      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2180    }    }
2181      exclusiveWrite();
2182    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2183       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2184    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2193  Data::setSlice(const Data& value,
2193    if (isProtected()) {    if (isProtected()) {
2194          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2195    }    }
2196    FORCERESOLVE;    forceResolve();
2197  /*  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.");  
   }*/  
2198    Data tempValue(value);    Data tempValue(value);
2199    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2200    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2431  Data::typeMatchRight(const Data& right) Line 2237  Data::typeMatchRight(const Data& right)
2237    }    }
2238  }  }
2239    
2240    // The normal TaggedValue adds the tag if it is not already present
2241    // This form does not. It throws instead.
2242    // This is because the names are maintained by the domain and cannot be added
2243    // without knowing the tag number to map it to.
2244  void  void
2245  Data::setTaggedValueByName(std::string name,  Data::setTaggedValueByName(std::string name,
2246                             const boost::python::object& value)                             const boost::python::object& value)
2247  {  {
2248       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2249      FORCERESOLVE;      forceResolve();
2250        exclusiveWrite();
2251          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2252          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2253       }       }
2254         else
2255         {                  // The
2256        throw DataException("Error - unknown tag in setTaggedValueByName.");
2257         }
2258  }  }
2259    
2260  void  void
2261  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2262                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2266  Data::setTaggedValue(int tagKey,
2266    }    }
2267    //    //
2268    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2269    FORCERESOLVE;    forceResolve();
2270      exclusiveWrite();
2271    if (isConstant()) tag();    if (isConstant()) tag();
2272    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]));  
   }  
2273    
2274    DataVector temp_data2;    DataVector temp_data2;
2275    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2276    
2277    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2278  }  }
2279    
2280    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2289  Data::setTaggedValueFromCPP(int tagKey,
2289    }    }
2290    //    //
2291    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2292    FORCERESOLVE;    forceResolve();
2293    if (isConstant()) tag();    if (isConstant()) tag();
2294      exclusiveWrite();
2295    //    //
2296    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2297    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2324  escript::C_GeneralTensorProduct(Data& ar
2324    // 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().
2325    
2326    // deal with any lazy data    // deal with any lazy data
2327    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2328    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2329      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2330      {
2331        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2332        return Data(c);
2333      }
2334    
2335    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2336    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2407  escript::C_GeneralTensorProduct(Data& ar
2407       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
2408    }    }
2409    
2410      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2411      {
2412         ostringstream os;
2413         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2414         throw DataException(os.str());
2415      }
2416    
2417    // Declare output Data object    // Declare output Data object
2418    Data res;    Data res;
2419    
2420    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2421      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2422      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2423      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2424      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2425      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);
2426    }    }
2427    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 2442  escript::C_GeneralTensorProduct(Data& ar
2442    
2443      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2444      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2445      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));  
2446    
2447        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2448        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2449    
2450      // Compute an MVP for the default      // Compute an MVP for the default
2451      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 2454  escript::C_GeneralTensorProduct(Data& ar
2454      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
2455      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2456        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);  
2457    
2458        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2459        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2460            
2461        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);
2462      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2480  escript::C_GeneralTensorProduct(Data& ar
2480        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2481          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2482          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2483          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2484          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2485          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2486          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);
2487        }        }
2488      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2506  escript::C_GeneralTensorProduct(Data& ar
2506    
2507      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2508      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2509      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2510      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2511  //     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));  
2512    
2513      // Compute an MVP for the default      // Compute an MVP for the default
2514      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 2516  escript::C_GeneralTensorProduct(Data& ar
2516      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2517      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
2518      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);  
2519    
2520        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2521        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2522        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2523        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);
2524      }      }
2525    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2540  escript::C_GeneralTensorProduct(Data& ar
2540      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2541      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2542    
2543  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2544  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2545  //     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));  
   
2546    
2547      // Compute an MVP for the default      // Compute an MVP for the default
2548      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 2559  escript::C_GeneralTensorProduct(Data& ar
2559      // Compute an MVP for each tag      // Compute an MVP for each tag
2560      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2561      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2562  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2563  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2564  //       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));  
2565    
2566        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);
2567      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2583  escript::C_GeneralTensorProduct(Data& ar
2583      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2584      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2585        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
2586        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2587        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2588          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2589          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2590          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2591          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2592          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);
2593        }        }
2594      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2612  escript::C_GeneralTensorProduct(Data& ar
2612        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2613          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2614          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2615          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2616          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2617          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2618          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);
2619        }        }
2620      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2637  escript::C_GeneralTensorProduct(Data& ar
2637      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2638      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2639        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2640        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2641        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2642          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2643          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2644          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2645          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2646          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);
2647        }        }
2648      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2667  escript::C_GeneralTensorProduct(Data& ar
2667          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2668          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2669          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2670          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2671          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2672          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2673          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);
2674        }        }
2675      }      }
# Line 2923  Data::borrowReadyPtr() const Line 2707  Data::borrowReadyPtr() const
2707  std::string  std::string
2708  Data::toString() const  Data::toString() const
2709  {  {
2710      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2711      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2712        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2713      {      {
2714      stringstream temp;      stringstream temp;
2715      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 2719  Data::toString() const
2719  }  }
2720    
2721    
2722    // This method is not thread-safe
2723    DataTypes::ValueType::reference
2724    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
2725    {
2726        checkExclusiveWrite();
2727        return getReady()->getDataAtOffsetRW(i);
2728    }
2729    
2730    // This method is not thread-safe
2731  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2732  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
2733  {  {
2734      if (isLazy())      forceResolve();
2735      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
2736  }  }
2737    
2738    
2739  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
2740  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
2741  {  // {
2742  //     if (isLazy())  //     if (isLazy())
2743  //     {  //     {
2744  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
2745  //     }  //     }
2746      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
2747      return getReady()->getDataAtOffset(i);  // }
2748  }  
2749    
2750  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2751  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
2752  {  {
2753      forceResolve();
2754    if (!isReady())    if (!isReady())
2755    {    {
2756      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.");
2757    }    }
2758    else    else
2759    {    {
2760      const DataReady* dr=getReady();      const DataReady* dr=getReady();
2761      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
2762    }    }
2763  }  }
2764    
2765    
2766    
2767    
2768  DataTypes::ValueType::reference  DataTypes::ValueType::reference
2769  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
2770  {  {
2771    FORCERESOLVE;    checkExclusiveWrite();
2772    if (!isReady())    DataReady* dr=getReady();
2773    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
2774      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
2775    }  
2776    else  BufferGroup*
2777    {  Data::allocSampleBuffer() const
2778      DataReady* dr=getReady();  {
2779      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
2780    }       {
2781        #ifdef _OPENMP
2782        int tnum=omp_get_max_threads();
2783        #else
2784        int tnum=1;
2785        #endif
2786        return new BufferGroup(getSampleBufferSize(),tnum);
2787         }
2788         else
2789         {
2790        return NULL;
2791         }
2792    }
2793    
2794    void
2795    Data::freeSampleBuffer(BufferGroup* bufferg)
2796    {
2797         if (bufferg!=0)
2798         {
2799        delete bufferg;
2800         }
2801  }  }
2802    
2803    
# Line 3000  Data::print() Line 2813  Data::print()
2813    {    {
2814      printf( "[%6d]", i );      printf( "[%6d]", i );
2815      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
2816        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
2817      printf( "\n" );      printf( "\n" );
2818    }    }
2819  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 2833  Data::dump(const std::string fileName) c
2833            return m_data->dump(fileName);            return m_data->dump(fileName);
2834      }      }
2835    }    }
2836    catch (exception& e)    catch (std::exception& e)
2837    {    {
2838          cout << e.what() << endl;          cout << e.what() << endl;
2839    }    }

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

  ViewVC Help
Powered by ViewVC 1.1.26