/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2037 by jfenwick, Thu Nov 13 06:17:12 2008 UTC revision 2271 by jfenwick, Mon Feb 16 05:08:29 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);    //numeric::array asNumArray(value);
327    initialise(asNumArray,what,expanded);    WrappedArray w(value);
328      initialise(w,what,expanded);
329    m_protected=false;    m_protected=false;
330  }  }
331    
332    
333  Data::Data(const object& value,  Data::Data(const object& value,
334             const Data& other)             const Data& other)
335        : m_shared(false), m_lazy(false)
336  {  {
337    numeric::array asNumArray(value);    WrappedArray w(value);
338    
339    // extract the shape of the numarray    // extract the shape of the numarray
340    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);    const DataTypes::ShapeType& tempShape=w.getShape();
341  // /*  for (int i=0; i < asNumArray.getrank(); i++) {    if (w.getRank()==0) {
 //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
 //   }*/  
 //   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 // /*  DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(asNumArray);*/  
 //   temp_data.copyFromNumArray(asNumArray);  
   
   //  
   // Create DataConstant using the given value and all other parameters  
   // copied from other. If value is a rank 0 object this Data  
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
342    
343    
344      // get the space for the data vector      // get the space for the data vector
345      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
346      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
347      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
348    
349      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
350    
351      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);  
   
352      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
353  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
354  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
355    
356    } else {    } else {
357      //      //
358      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
359  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
360      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
361  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
362    }    }
363    m_protected=false;    m_protected=false;
364  }  }
365    
366  Data::~Data()  Data::~Data()
367  {  {
368      set_m_data(DataAbstract_ptr());
369  }  }
370    
371    
372    // only call in thread safe contexts.
373    // This method should be atomic
374    void Data::set_m_data(DataAbstract_ptr p)
375    {
376      if (m_data.get()!=0)  // release old ownership
377      {
378        m_data->removeOwner(this);
379      }
380      if (p.get()!=0)
381      {
382        m_data=p;
383        m_data->addOwner(this);
384        m_shared=m_data->isShared();
385        m_lazy=m_data->isLazy();
386      }
387    }
388    
389  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
390                   const FunctionSpace& what,                   const FunctionSpace& what,
391                   bool expanded)                   bool expanded)
392  {  {
# Line 277  Data::initialise(const boost::python::nu Line 397  Data::initialise(const boost::python::nu
397    // within the shared_ptr constructor.    // within the shared_ptr constructor.
398    if (expanded) {    if (expanded) {
399      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
400  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
401  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
402    } else {    } else {
403      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
404  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
405  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
406    }    }
407  }  }
408    
# Line 302  Data::initialise(const DataTypes::ValueT Line 420  Data::initialise(const DataTypes::ValueT
420    // within the shared_ptr constructor.    // within the shared_ptr constructor.
421    if (expanded) {    if (expanded) {
422      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
423  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
424  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
425    } else {    } else {
426      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
427  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
428  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
429    }    }
430  }  }
431    
432    
 // 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";  
 //  }  
 // }  
   
433  escriptDataC  escriptDataC
434  Data::getDataC()  Data::getDataC()
435  {  {
# Line 410  Data::getDataC() const Line 446  Data::getDataC() const
446    return temp;    return temp;
447  }  }
448    
449    size_t
450    Data::getSampleBufferSize() const
451    {
452      return m_data->getSampleBufferSize();
453    }
454    
455    
456  const boost::python::tuple  const boost::python::tuple
457  Data::getShapeTuple() const  Data::getShapeTuple() const
458  {  {
# Line 435  Data::getShapeTuple() const Line 478  Data::getShapeTuple() const
478  // 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.
479  // 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
480  // 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
481  Data*  Data
482  Data::copySelf()  Data::copySelf()
483  {  {
484     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
485     return new Data(temp);     return Data(temp);
486  }  }
487    
488  void  void
# Line 447  Data::copy(const Data& other) Line 490  Data::copy(const Data& other)
490  {  {
491    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
492    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
493    m_data=p;  //   m_data=p;
494      set_m_data(p);
495  }  }
496    
497    
# Line 463  Data::delaySelf() Line 507  Data::delaySelf()
507  {  {
508    if (!isLazy())    if (!isLazy())
509    {    {
510      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
511        set_m_data((new DataLazy(m_data))->getPtr());
512    }    }
513  }  }
514    
515    
516    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
517    // to zero, all the tags from all the DataTags would be in the result.
518    // However since they all have the same value (0) whether they are there or not should not matter.
519    // So I have decided that for all types this method will create a constant 0.
520    // It can be promoted up as required.
521    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
522    // but we can deal with that if it arises.
523    //
524  void  void
525  Data::setToZero()  Data::setToZero()
526  {  {
# Line 474  Data::setToZero() Line 528  Data::setToZero()
528    {    {
529       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
530    }    }
531    m_data->setToZero();    if (isLazy())
532      {
533         DataTypes::ValueType v(getNoValues(),0);
534         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
535         DataLazy* dl=new DataLazy(dc->getPtr());
536         set_m_data(dl->getPtr());
537      }
538      else
539      {
540         exclusiveWrite();
541         m_data->setToZero();
542      }
543  }  }
544    
545    
546  void  void
547  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
548                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 604  Data::copyWithMask(const Data& other,
604    {    {
605      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
606    }    }
607      unsigned int selfrank=getDataPointRank();
608      unsigned int otherrank=other2.getDataPointRank();
609      unsigned int maskrank=mask2.getDataPointRank();
610      if ((selfrank==0) && (otherrank>0 || maskrank>0))
611      {
612        // to get here we must be copying from a large object into a scalar
613        // I am not allowing this.
614        // If you are calling copyWithMask then you are considering keeping some existing values
615        // and so I'm going to assume that you don't want your data objects getting a new shape.
616        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
617      }
618      exclusiveWrite();
619    // Now we iterate over the elements    // Now we iterate over the elements
620    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
621    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
622    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
623    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
624      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
625    {    {
626      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // Not allowing this combination.
627        // it is not clear what the rank of the target should be.
628        // Should it be filled with the scalar (rank stays the same);
629        // or should the target object be reshaped to be a scalar as well.
630        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
631      }
632      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
633      {
634        if (mvec[0]>0)      // copy whole object if scalar is >0
635        {
636            copy(other);
637        }
638        return;
639      }
640      if (isTagged())       // so all objects involved will also be tagged
641      {
642        // note the !
643        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
644            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
645        {
646            throw DataException("copyWithMask, shape mismatch.");
647        }
648    
649        // We need to consider the possibility that tags are missing or in the wrong order
650        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
651        // all tags which are not explicitly defined
652    
653        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
654        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
655        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
656    
657        // first, add any tags missing from other or mask
658        const DataTagged::DataMapType& olookup=optr->getTagLookup();
659            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
660        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
661        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
662        for (i=olookup.begin();i!=olookup.end();i++)
663        {
664               tptr->addTag(i->first);
665            }
666            for (i=mlookup.begin();i!=mlookup.end();i++) {
667               tptr->addTag(i->first);
668            }
669        // now we know that *this has all the required tags but they aren't guaranteed to be in
670        // the same order
671    
672        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
673        if ((selfrank==otherrank) && (otherrank==maskrank))
674        {
675            for (i=tlookup.begin();i!=tlookup.end();i++)
676            {
677                // get the target offset
678                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
679                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
680                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
681                for (int j=0;j<getDataPointSize();++j)
682                {
683                    if (mvec[j+moff]>0)
684                    {
685                        self[j+toff]=ovec[j+ooff];
686                    }
687                }
688                }
689            // now for the default value
690            for (int j=0;j<getDataPointSize();++j)
691            {
692                if (mvec[j+mptr->getDefaultOffset()]>0)
693                {
694                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
695                }
696            }
697        }
698        else    // other is a scalar
699        {
700            for (i=tlookup.begin();i!=tlookup.end();i++)
701            {
702                // get the target offset
703                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
704                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
705                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
706                for (int j=0;j<getDataPointSize();++j)
707                {
708                    if (mvec[j+moff]>0)
709                    {
710                        self[j+toff]=ovec[ooff];
711                    }
712                }
713                }
714            // now for the default value
715            for (int j=0;j<getDataPointSize();++j)
716            {
717                if (mvec[j+mptr->getDefaultOffset()]>0)
718                {
719                    self[j+tptr->getDefaultOffset()]=ovec[0];
720                }
721            }
722        }
723    
724        return;         // ugly
725      }
726      // mixed scalar and non-scalar operation
727      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
728      {
729            size_t num_points=self.size();
730        // OPENMP 3.0 allows unsigned loop vars.
731        #if defined(_OPENMP) && (_OPENMP < 200805)
732        long i;
733        #else
734        size_t i;
735        #endif  
736        size_t psize=getDataPointSize();    
737        #pragma omp parallel for private(i) schedule(static)
738        for (i=0;i<num_points;++i)
739        {
740            if (mvec[i]>0)
741            {
742                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
743            }                   // dest point
744        }
745        return;         // ugly!
746      }
747      // tagged data is already taken care of so we only need to worry about shapes
748      // special cases with scalars are already dealt with so all we need to worry about is shape
749      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
750      {
751        ostringstream oss;
752        oss <<"Error - size mismatch in arguments to copyWithMask.";
753        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
754        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
755        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
756        throw DataException(oss.str());
757    }    }
758    size_t num_points=self.size();    size_t num_points=self.size();
759    
# Line 564  Data::copyWithMask(const Data& other, Line 773  Data::copyWithMask(const Data& other,
773    }    }
774  }  }
775    
   
   
776  bool  bool
777  Data::isExpanded() const  Data::isExpanded() const
778  {  {
# Line 574  Data::isExpanded() const Line 781  Data::isExpanded() const
781  }  }
782    
783  bool  bool
784    Data::actsExpanded() const
785    {
786      return m_data->actsExpanded();
787    
788    }
789    
790    
791    bool
792  Data::isTagged() const  Data::isTagged() const
793  {  {
794    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 812  Data::isConstant() const
812  bool  bool
813  Data::isLazy() const  Data::isLazy() const
814  {  {
815    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
816  }  }
817    
818  // 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 843  Data::expand()
843    if (isConstant()) {    if (isConstant()) {
844      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
845      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
846  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
847  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
848    } else if (isTagged()) {    } else if (isTagged()) {
849      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
850      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
851  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
852  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
853    } else if (isExpanded()) {    } else if (isExpanded()) {
854      //      //
855      // do nothing      // do nothing
# Line 656  Data::tag() Line 869  Data::tag()
869    if (isConstant()) {    if (isConstant()) {
870      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
871      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
872  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
873  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
874    } else if (isTagged()) {    } else if (isTagged()) {
875      // do nothing      // do nothing
876    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 883  Data::tag()
883       {       {
884      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
885       }       }
886       m_data=res;      //      m_data=res;
887         set_m_data(res);
888       tag();       tag();
889    } else {    } else {
890      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 896  Data::resolve()
896  {  {
897    if (isLazy())    if (isLazy())
898    {    {
899       m_data=m_data->resolve();  //      m_data=m_data->resolve();
900        set_m_data(m_data->resolve());
901    }    }
902  }  }
903    
904    void
905    Data::requireWrite()
906    {
907      resolve();
908      exclusiveWrite();
909    }
910    
911  Data  Data
912  Data::oneOver() const  Data::oneOver() const
913  {  {
914    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
915    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
916  }  }
917    
918  Data  Data
919  Data::wherePositive() const  Data::wherePositive() const
920  {  {
921    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
922    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
923  }  }
924    
925  Data  Data
926  Data::whereNegative() const  Data::whereNegative() const
927  {  {
928    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
929    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
930  }  }
931    
932  Data  Data
933  Data::whereNonNegative() const  Data::whereNonNegative() const
934  {  {
935    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
936    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
937  }  }
938    
939  Data  Data
940  Data::whereNonPositive() const  Data::whereNonPositive() const
941  {  {
942    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
943    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
944  }  }
945    
946  Data  Data
947  Data::whereZero(double tol) const  Data::whereZero(double tol) const
948  {  {
949    Data dataAbs=abs();  //   Data dataAbs=abs();
950    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
951       MAKELAZYOPOFF(EZ,tol)
952       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
953    
954  }  }
955    
956  Data  Data
957  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
958  {  {
959    Data dataAbs=abs();  //   Data dataAbs=abs();
960    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
961      MAKELAZYOPOFF(NEZ,tol)
962      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
963    
964  }  }
965    
966  Data  Data
# Line 767  bool Line 973  bool
973  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
974  {  {
975    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());  
 //     }  
 //   }  
976  }  }
977    
978  Data  Data
# Line 820  Data::getLength() const Line 1016  Data::getLength() const
1016  }  }
1017    
1018  const  const
1019  boost::python::numeric::array  boost::python::numeric :: array
1020  Data:: getValueOfDataPoint(int dataPointNo)  Data:: getValueOfDataPoint(int dataPointNo)
1021  {  {
1022    int i, j, k, l;      return boost::python::numeric::array(getValueOfDataPointAsTuple(dataPointNo));
1023    }
   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]);  
   }  
1024    
1025    if (getNumDataPointsPerSample()>0) {  const boost::python::object
1026    Data::getValueOfDataPointAsTuple(int dataPointNo)
1027    {
1028       forceResolve();
1029       if (getNumDataPointsPerSample()>0) {
1030         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1031         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1032         //         //
1033         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1034         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1035             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1036         }         }
1037    
1038         //         //
1039         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1040         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1041             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1042         }         }
1043         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1044         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1045           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1046      }
1047         switch( dataPointRank ){    else
1048              case 0 :    {
1049                  numArray[0] = getDataAtOffset(offset);      // The old method would return an empty numArray of the given shape
1050                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1051              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;  
     }  
1052    }    }
   //  
   // return the array  
   return numArray;  
1053    
1054  }  }
1055    
1056    
1057  void  void
1058  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1059  {  {
1060      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1061      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1062  }  }
1063    
1064  void  void
1065  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1066  {  {
1067    if (isProtected()) {    if (isProtected()) {
1068          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1069    }    }
1070    FORCERESOLVE;    forceResolve();
1071    
1072      WrappedArray w(obj);
1073    //    //
1074    // check rank    // check rank
1075    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1076        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of numarray does not match Data object rank");
1077    
1078    //    //
1079    // check shape of num_array    // check shape of num_array
1080    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1081      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1082         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of numarray does not match Data object rank");
1083    }    }
1084    //    //
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1090  Data::setValueOfDataPointToArray(int dat
1090    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1091         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1092         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1093         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1094    } else {    } else {
1095         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1096    }    }
1097  }  }
1098    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1104  Data::setValueOfDataPoint(int dataPointN
1104    }    }
1105    //    //
1106    // make sure data is expanded:    // make sure data is expanded:
1107    FORCERESOLVE;    forceResolve();
1108    if (!isExpanded()) {    if (!isExpanded()) {
1109      expand();      expand();
1110    }    }
# Line 980  const Line 1121  const
1121  boost::python::numeric::array  boost::python::numeric::array
1122  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
1123  {  {
1124    size_t length=0;     return boost::python::numeric::array(getValueOfGlobalDataPointAsTuple(procNo, dataPointNo));
1125    int i, j, k, l, pos;  }
   FORCERESOLVE;  
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
1126    
1127    //  const
1128    // create the numeric array to be returned  boost::python::object
1129    boost::python::numeric::array numArray(0.0);  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1130    {
1131      // This could be lazier than it is now
1132      forceResolve();
1133    
1134    //    // copy datapoint into a buffer
1135    // the shape of the returned numeric array will be the same    // broadcast buffer to all nodes
1136    // as that of the data point    // convert buffer to tuple
1137    int arrayRank = dataPointRank;    // return tuple
   const DataTypes::ShapeType& arrayShape = dataPointShape;  
1138    
1139    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1140    // 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]);  
   }  
1141    
1142    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1143    double *tmpData = new double[length];    double *tmpData = new double[length];
1144    
1145    //    // updated for the MPI case
1146    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1147          if (getNumDataPointsPerSample()>0) {
1148      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1149      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1150               if (getNumDataPointsPerSample()>0) {      //
1151                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1152                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1153                  //          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;  
         }  
             }  
1154      }      }
1155          #ifdef PASO_MPI  
1156          // broadcast the data to all other processes      //
1157      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1158          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1159            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;  
1160      }      }
1161        // TODO: global error handling
1162        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1163    
1164      delete [] tmpData;      memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1165         }
1166      }
1167    #ifdef PASO_MPI
1168      // broadcast the data to all other processes
1169      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1170    #endif
1171    
1172      boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1173      delete [] tmpData;
1174    //    //
1175    // return the loaded array    // return the loaded array
1176    return numArray;    return t;
1177    
1178  }  }
1179    
1180    
1181  boost::python::numeric::array  boost::python::object
1182  Data::integrate_const() const  Data::integrateToTuple_const() const
1183  {  {
1184    if (isLazy())    if (isLazy())
1185    {    {
# Line 1127  Data::integrate_const() const Line 1188  Data::integrate_const() const
1188    return integrateWorker();    return integrateWorker();
1189  }  }
1190    
1191  boost::python::numeric::array  boost::python::object
1192  Data::integrate()  Data::integrateToTuple()
1193  {  {
1194    if (isLazy())    if (isLazy())
1195    {    {
1196      expand();      expand();
1197    }    }
1198    return integrateWorker();    return integrateWorker();
1199    
1200  }  }
1201    
1202    
1203    boost::python::object
1204    Data::integrate_const() const
1205    {
1206      if (isLazy())
1207      {
1208        throw DataException("Error - cannot integrate for constant lazy data.");
1209      }
1210      return boost::python::numeric::array(integrateWorker());
1211    }
1212    
1213  boost::python::numeric::array  boost::python::object
1214    Data::integrate()
1215    {
1216      if (isLazy())
1217      {
1218        expand();
1219      }
1220      return boost::python::numeric::array(integrateWorker());
1221    }
1222    
1223    
1224    
1225    boost::python::object
1226  Data::integrateWorker() const  Data::integrateWorker() const
1227  {  {
   int index;  
   int rank = getDataPointRank();  
1228    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1229    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1230    
# Line 1151  Data::integrateWorker() const Line 1232  Data::integrateWorker() const
1232    // calculate the integral values    // calculate the integral values
1233    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1234    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1235      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1236      if (dom==0)
1237      {            
1238        throw DataException("Can not integrate over non-continuous domains.");
1239      }
1240  #ifdef PASO_MPI  #ifdef PASO_MPI
1241    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1242    // 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
1243    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1244    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1245    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]; }
1246    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 );
1247    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1248      tuple result=pointToTuple(shape,tmp);
1249    delete[] tmp;    delete[] tmp;
1250    delete[] tmp_local;    delete[] tmp_local;
1251  #else  #else
1252    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1253    /*  double *tmp = new double[dataPointSize];
1254      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1255      tuple result=pointToTuple(shape,integrals);
1256    //   delete tmp;
1257  #endif  #endif
1258    
   //  
   // 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];  
           }  
         }  
       }  
     }  
   }  
1259    
1260    //    return result;
   // return the loaded array  
   return bp_array;  
1261  }  }
1262    
1263  Data  Data
1264  Data::sin() const  Data::sin() const
1265  {  {
1266    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1267    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1268  }  }
1269    
1270  Data  Data
1271  Data::cos() const  Data::cos() const
1272  {  {
1273    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1274    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1275  }  }
1276    
1277  Data  Data
1278  Data::tan() const  Data::tan() const
1279  {  {
1280    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1281    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1282  }  }
1283    
1284  Data  Data
1285  Data::asin() const  Data::asin() const
1286  {  {
1287    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1288    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1289  }  }
1290    
1291  Data  Data
1292  Data::acos() const  Data::acos() const
1293  {  {
1294    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1295    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1296  }  }
1297    
# Line 1279  Data::acos() const Line 1299  Data::acos() const
1299  Data  Data
1300  Data::atan() const  Data::atan() const
1301  {  {
1302    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1303    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1304  }  }
1305    
1306  Data  Data
1307  Data::sinh() const  Data::sinh() const
1308  {  {
1309    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1310    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1311  }  }
1312    
1313  Data  Data
1314  Data::cosh() const  Data::cosh() const
1315  {  {
1316    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1317    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1318  }  }
1319    
1320  Data  Data
1321  Data::tanh() const  Data::tanh() const
1322  {  {
1323    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1324    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1325  }  }
1326    
# Line 1324  Data::tanh() const Line 1328  Data::tanh() const
1328  Data  Data
1329  Data::erf() const  Data::erf() const
1330  {  {
1331  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1332    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1333  #else  #else
1334    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1335    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1336  #endif  #endif
1337  }  }
# Line 1339  Data::erf() const Line 1339  Data::erf() const
1339  Data  Data
1340  Data::asinh() const  Data::asinh() const
1341  {  {
1342    if (isLazy())    MAKELAZYOP(ASINH)
1343    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1344    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1345  #else  #else
1346    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1354  Data::asinh() const Line 1350  Data::asinh() const
1350  Data  Data
1351  Data::acosh() const  Data::acosh() const
1352  {  {
1353    if (isLazy())    MAKELAZYOP(ACOSH)
1354    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1355    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1356  #else  #else
1357    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1369  Data::acosh() const Line 1361  Data::acosh() const
1361  Data  Data
1362  Data::atanh() const  Data::atanh() const
1363  {  {
1364    if (isLazy())    MAKELAZYOP(ATANH)
1365    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1366    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1367  #else  #else
1368    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1383  Data::atanh() const Line 1371  Data::atanh() const
1371    
1372  Data  Data
1373  Data::log10() const  Data::log10() const
1374  {  if (isLazy())  {
1375    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1376    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1377  }  }
1378    
1379  Data  Data
1380  Data::log() const  Data::log() const
1381  {  {
1382    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1383    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1384  }  }
1385    
1386  Data  Data
1387  Data::sign() const  Data::sign() const
1388  {  {
1389    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1390    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1391  }  }
1392    
1393  Data  Data
1394  Data::abs() const  Data::abs() const
1395  {  {
1396    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1397    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1398  }  }
1399    
1400  Data  Data
1401  Data::neg() const  Data::neg() const
1402  {  {
1403    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1404    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1405  }  }
1406    
# Line 1448  Data::pos() const Line 1417  Data::pos() const
1417    
1418  Data  Data
1419  Data::exp() const  Data::exp() const
1420  {    {
1421    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1422    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1423  }  }
1424    
1425  Data  Data
1426  Data::sqrt() const  Data::sqrt() const
1427  {  {
1428    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1429    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1430  }  }
1431    
# Line 1483  Data::Lsup() Line 1444  Data::Lsup()
1444  {  {
1445     if (isLazy())     if (isLazy())
1446     {     {
1447      expand();      resolve();
1448     }     }
1449     return LsupWorker();     return LsupWorker();
1450  }  }
# Line 1503  Data::sup() Line 1464  Data::sup()
1464  {  {
1465     if (isLazy())     if (isLazy())
1466     {     {
1467      expand();      resolve();
1468     }     }
1469     return supWorker();     return supWorker();
1470  }  }
# Line 1523  Data::inf() Line 1484  Data::inf()
1484  {  {
1485     if (isLazy())     if (isLazy())
1486     {     {
1487      expand();      resolve();
1488     }     }
1489     return infWorker();     return infWorker();
1490  }  }
# Line 1615  Data::minval() const Line 1576  Data::minval() const
1576  Data  Data
1577  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1578  {  {
1579         if (isLazy())
1580         {
1581        Data temp(*this);
1582        temp.resolve();
1583        return temp.swapaxes(axis0,axis1);
1584         }
1585       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1586       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1587       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1672  Data::symmetric() const Line 1639  Data::symmetric() const
1639       else {       else {
1640          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.");
1641       }       }
1642       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1643       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1644       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1645       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1649  Data::symmetric() const
1649  Data  Data
1650  Data::nonsymmetric() const  Data::nonsymmetric() const
1651  {  {
1652       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1653       // check input       // check input
1654       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1655       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1681  Data::nonsymmetric() const
1681       }       }
1682  }  }
1683    
   
 // 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).  
1684  Data  Data
1685  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1686  {  {    
1687       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1688         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1689       {       {
1690      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);  
1691       }       }
1692       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1693       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1739  Data::trace(int axis_offset) const
1739  Data  Data
1740  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1741  {      {    
1742       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);  
      }  
1743       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1744       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1745       // 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 1841  Data::calc_minGlobalDataPoint(int& ProcN
1841    double next,local_min;    double next,local_min;
1842    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1843    
1844    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1845    {    {
1846      local_min=min;      local_min=min;
1847      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1848      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1849        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1850          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1851          if (next<local_min) {          if (next<local_min) {
1852            local_min=next;            local_min=next;
1853            local_lowi=i;            local_lowi=i;
# Line 1918  Data::calc_minGlobalDataPoint(int& ProcN Line 1865  Data::calc_minGlobalDataPoint(int& ProcN
1865    
1866  #ifdef PASO_MPI  #ifdef PASO_MPI
1867      // determine the processor on which the minimum occurs      // determine the processor on which the minimum occurs
1868      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPointRO(lowi,lowj);
1869      int lowProc = 0;      int lowProc = 0;
1870      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1871      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1872        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1873    
1874      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1875          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 1987  Data::operator+=(const Data& right) Line 1935  Data::operator+=(const Data& right)
1935    if (isProtected()) {    if (isProtected()) {
1936          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1937    }    }
1938    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
1939    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
1940      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
1941          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
1942  }  }
1943    
1944  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 1948  Data::operator+=(const boost::python::ob
1948          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1949    }    }
1950    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1951    if (isLazy())    (*this)+=tmp;
1952    {    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);  
   }  
1953  }  }
1954    
1955  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
1956  Data&  Data&
1957  Data::operator=(const Data& other)  Data::operator=(const Data& other)
1958  {  {
1959    #if defined ASSIGNMENT_MEANS_DEEPCOPY  
1960    // This should not be used.
1961    // Just leaving this here until I have completed transition
1962    copy(other);    copy(other);
1963    #else
1964      m_protected=false;        // since any changes should be caught by exclusiveWrite();
1965    //   m_data=other.m_data;
1966      set_m_data(other.m_data);
1967    #endif
1968    return (*this);    return (*this);
1969  }  }
1970    
# Line 2034  Data::operator-=(const Data& right) Line 1974  Data::operator-=(const Data& right)
1974    if (isProtected()) {    if (isProtected()) {
1975          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1976    }    }
1977    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
1978    {    exclusiveWrite();
1979      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
1980          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
1981  }  }
1982    
1983  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 1987  Data::operator-=(const boost::python::ob
1987          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1988    }    }
1989    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
1990    if (isLazy())    (*this)-=tmp;
1991    {    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);  
   }  
1992  }  }
1993    
1994  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 1997  Data::operator*=(const Data& right)
1997    if (isProtected()) {    if (isProtected()) {
1998          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1999    }    }
2000    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2001    {    exclusiveWrite();
2002      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
2003          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
2004  }  }
2005    
2006  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2010  Data::operator*=(const boost::python::ob
2010          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2011    }    }
2012    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2013    if (isLazy())    (*this)*=tmp;
2014    {    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);  
   }  
2015  }  }
2016    
2017  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2020  Data::operator/=(const Data& right)
2020    if (isProtected()) {    if (isProtected()) {
2021          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2022    }    }
2023    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2024    {    exclusiveWrite();
2025      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
2026          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
2027  }  }
2028    
2029  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2033  Data::operator/=(const boost::python::ob
2033          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2034    }    }
2035    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2036    if (isLazy())    (*this)/=tmp;
2037    {    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);  
   }  
2038  }  }
2039    
2040  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2054  Data::powO(const boost::python::object&
2054  Data  Data
2055  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2056  {  {
2057    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2058    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2059  }  }
2060    
# Line 2175  Data::powD(const Data& right) const Line 2063  Data::powD(const Data& right) const
2063  Data  Data
2064  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2065  {  {
2066    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2067    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2068  }  }
2069    
# Line 2188  escript::operator+(const Data& left, con Line 2072  escript::operator+(const Data& left, con
2072  Data  Data
2073  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2074  {  {
2075    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2076    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2077  }  }
2078    
# Line 2201  escript::operator-(const Data& left, con Line 2081  escript::operator-(const Data& left, con
2081  Data  Data
2082  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2083  {  {
2084    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2085    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2086  }  }
2087    
# Line 2214  escript::operator*(const Data& left, con Line 2090  escript::operator*(const Data& left, con
2090  Data  Data
2091  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2092  {  {
2093    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2094    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2095  }  }
2096    
# Line 2227  escript::operator/(const Data& left, con Line 2099  escript::operator/(const Data& left, con
2099  Data  Data
2100  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2101  {  {
2102    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2103    {    MAKELAZYBIN2(left,tmp,ADD)
2104      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);  
2105  }  }
2106    
2107  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2109  escript::operator+(const Data& left, con
2109  Data  Data
2110  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2111  {  {
2112    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2113    {    MAKELAZYBIN2(left,tmp,SUB)
2114      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);  
2115  }  }
2116    
2117  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2119  escript::operator-(const Data& left, con
2119  Data  Data
2120  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2121  {  {
2122    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2123    {    MAKELAZYBIN2(left,tmp,MUL)
2124      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);  
2125  }  }
2126    
2127  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2129  escript::operator*(const Data& left, con
2129  Data  Data
2130  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2131  {  {
2132    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2133    {    MAKELAZYBIN2(left,tmp,DIV)
2134      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);  
2135  }  }
2136    
2137  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2139  escript::operator/(const Data& left, con
2139  Data  Data
2140  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2141  {  {
2142    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2143    {    MAKELAZYBIN2(tmp,right,ADD)
2144      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;  
2145  }  }
2146    
2147  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2149  escript::operator+(const boost::python::
2149  Data  Data
2150  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2151  {  {
2152    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2153    {    MAKELAZYBIN2(tmp,right,SUB)
2154      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;  
2155  }  }
2156    
2157  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2159  escript::operator-(const boost::python::
2159  Data  Data
2160  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2161  {  {
2162    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2163    {    MAKELAZYBIN2(tmp,right,MUL)
2164      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;  
2165  }  }
2166    
2167  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2169  escript::operator*(const boost::python::
2169  Data  Data
2170  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2171  {  {
2172    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2173    {    MAKELAZYBIN2(tmp,right,DIV)
2174      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;  
2175  }  }
2176    
2177    
# Line 2364  void Line 2212  void
2212  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2213                 const Data& value)                 const Data& value)
2214  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2215    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2216    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2217      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2218    }    }
2219      exclusiveWrite();
2220    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2221       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2222    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2231  Data::setSlice(const Data& value,
2231    if (isProtected()) {    if (isProtected()) {
2232          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2233    }    }
2234    FORCERESOLVE;    forceResolve();
2235  /*  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.");  
   }*/  
2236    Data tempValue(value);    Data tempValue(value);
2237    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2238    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2436  Data::setTaggedValueByName(std::string n Line 2280  Data::setTaggedValueByName(std::string n
2280                             const boost::python::object& value)                             const boost::python::object& value)
2281  {  {
2282       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2283      FORCERESOLVE;      forceResolve();
2284        exclusiveWrite();
2285          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2286          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2287       }       }
2288  }  }
2289    
2290  void  void
2291  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2292                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2296  Data::setTaggedValue(int tagKey,
2296    }    }
2297    //    //
2298    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2299    FORCERESOLVE;    forceResolve();
2300      exclusiveWrite();
2301    if (isConstant()) tag();    if (isConstant()) tag();
2302    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]));  
   }  
2303    
2304    DataVector temp_data2;    DataVector temp_data2;
2305    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2306    
2307    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2308  }  }
2309    
2310    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2319  Data::setTaggedValueFromCPP(int tagKey,
2319    }    }
2320    //    //
2321    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2322    FORCERESOLVE;    forceResolve();
2323    if (isConstant()) tag();    if (isConstant()) tag();
2324      exclusiveWrite();
2325    //    //
2326    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2327    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2354  escript::C_GeneralTensorProduct(Data& ar
2354    // 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().
2355    
2356    // deal with any lazy data    // deal with any lazy data
2357    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2358    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2359      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2360      {
2361        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2362        return Data(c);
2363      }
2364    
2365    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2366    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2437  escript::C_GeneralTensorProduct(Data& ar
2437       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
2438    }    }
2439    
2440      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2441      {
2442         ostringstream os;
2443         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2444         throw DataException(os.str());
2445      }
2446    
2447    // Declare output Data object    // Declare output Data object
2448    Data res;    Data res;
2449    
2450    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2451      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2452      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2453      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2454      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2455      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);
2456    }    }
2457    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 2472  escript::C_GeneralTensorProduct(Data& ar
2472    
2473      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2474      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2475      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));  
2476    
2477        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2478        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2479    
2480      // Compute an MVP for the default      // Compute an MVP for the default
2481      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 2484  escript::C_GeneralTensorProduct(Data& ar
2484      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
2485      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2486        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);  
2487    
2488        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2489        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2490            
2491        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);
2492      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2510  escript::C_GeneralTensorProduct(Data& ar
2510        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2511          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2512          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2513          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2514          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2515          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2516          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);
2517        }        }
2518      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2536  escript::C_GeneralTensorProduct(Data& ar
2536    
2537      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2538      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2539      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2540      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2541  //     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));  
2542    
2543      // Compute an MVP for the default      // Compute an MVP for the default
2544      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 2546  escript::C_GeneralTensorProduct(Data& ar
2546      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2547      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
2548      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);  
2549    
2550        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2551        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2552        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2553        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);
2554      }      }
2555    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2570  escript::C_GeneralTensorProduct(Data& ar
2570      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2571      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2572    
2573  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2574  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2575  //     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));  
   
2576    
2577      // Compute an MVP for the default      // Compute an MVP for the default
2578      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 2589  escript::C_GeneralTensorProduct(Data& ar
2589      // Compute an MVP for each tag      // Compute an MVP for each tag
2590      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2591      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2592  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2593  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2594  //       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));  
2595    
2596        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);
2597      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2613  escript::C_GeneralTensorProduct(Data& ar
2613      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2614      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2615        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
2616        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2617        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2618          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2619          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2620          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2621          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2622          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);
2623        }        }
2624      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2642  escript::C_GeneralTensorProduct(Data& ar
2642        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2643          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2644          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2645          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2646          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2647          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2648          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);
2649        }        }
2650      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2667  escript::C_GeneralTensorProduct(Data& ar
2667      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2668      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2669        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2670        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2671        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2672          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2673          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2674          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2675          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2676          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);
2677        }        }
2678      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2697  escript::C_GeneralTensorProduct(Data& ar
2697          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2698          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2699          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2700          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2701          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2702          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2703          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);
2704        }        }
2705      }      }
# Line 2923  Data::borrowReadyPtr() const Line 2737  Data::borrowReadyPtr() const
2737  std::string  std::string
2738  Data::toString() const  Data::toString() const
2739  {  {
2740      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2741      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2742        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2743      {      {
2744      stringstream temp;      stringstream temp;
2745      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 2749  Data::toString() const
2749  }  }
2750    
2751    
2752    // This method is not thread-safe
2753    DataTypes::ValueType::reference
2754    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
2755    {
2756        checkExclusiveWrite();
2757        return getReady()->getDataAtOffsetRW(i);
2758    }
2759    
2760    // This method is not thread-safe
2761  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2762  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
2763  {  {
2764      if (isLazy())      forceResolve();
2765      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
2766  }  }
2767    
2768    
2769  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
2770  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
2771  {  // {
2772  //     if (isLazy())  //     if (isLazy())
2773  //     {  //     {
2774  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
2775  //     }  //     }
2776      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
2777      return getReady()->getDataAtOffset(i);  // }
2778  }  
2779    
2780  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2781  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
2782  {  {
2783      forceResolve();
2784    if (!isReady())    if (!isReady())
2785    {    {
2786      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.");
2787    }    }
2788    else    else
2789    {    {
2790      const DataReady* dr=getReady();      const DataReady* dr=getReady();
2791      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
2792    }    }
2793  }  }
2794    
2795    
2796    
2797    
2798  DataTypes::ValueType::reference  DataTypes::ValueType::reference
2799  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
2800  {  {
2801    FORCERESOLVE;    checkExclusiveWrite();
2802    if (!isReady())    DataReady* dr=getReady();
2803    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
2804      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
2805    }  
2806    else  BufferGroup*
2807    {  Data::allocSampleBuffer() const
2808      DataReady* dr=getReady();  {
2809      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
2810    }       {
2811        #ifdef _OPENMP
2812        int tnum=omp_get_max_threads();
2813        #else
2814        int tnum=1;
2815        #endif
2816        return new BufferGroup(getSampleBufferSize(),tnum);
2817         }
2818         else
2819         {
2820        return NULL;
2821         }
2822    }
2823    
2824    void
2825    Data::freeSampleBuffer(BufferGroup* bufferg)
2826    {
2827         if (bufferg!=0)
2828         {
2829        delete bufferg;
2830         }
2831  }  }
2832    
2833    
# Line 3000  Data::print() Line 2843  Data::print()
2843    {    {
2844      printf( "[%6d]", i );      printf( "[%6d]", i );
2845      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
2846        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
2847      printf( "\n" );      printf( "\n" );
2848    }    }
2849  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 2863  Data::dump(const std::string fileName) c
2863            return m_data->dump(fileName);            return m_data->dump(fileName);
2864      }      }
2865    }    }
2866    catch (exception& e)    catch (std::exception& e)
2867    {    {
2868          cout << e.what() << endl;          cout << e.what() << endl;
2869    }    }

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

  ViewVC Help
Powered by ViewVC 1.1.26