/[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 2260 by jfenwick, Tue Feb 10 04:50:10 2009 UTC revision 2271 by jfenwick, Mon Feb 16 05:08:29 2009 UTC
# Line 38  extern "C" { Line 38  extern "C" {
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 46  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()  #define AUTOLAZYON escriptParams.getAUTOLAZY()
53  #define MAKELAZYOP(X)   if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \  #define MAKELAZYOP(X)   if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
54    {\    {\
# Line 62  using namespace escript; Line 64  using namespace escript;
64  #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \  #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
65    {\    {\
66      DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\      DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
67          m_data=c->getPtr();\  /*         m_data=c->getPtr();*/     set_m_data(c->getPtr());\
68      return (*this);\      return (*this);\
69    }    }
70    
# Line 79  using namespace escript; Line 81  using namespace escript;
81      return Data(c);\      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 92  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 108  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 121  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 142  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 166  Data::Data(const Data& inData, Line 270  Data::Data(const Data& inData,
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 181  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;  
 }  
   
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 222  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    if (DataTypes::getRank(tempShape)==0) {    if (w.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,1);      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);
352      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
353      m_data=DataAbstract_ptr(t);  //     m_data=DataAbstract_ptr(t);
354        set_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      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
360  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
361  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     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 280  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 305  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 413  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 450  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 466  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 477  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 552  Data::copyWithMask(const Data& other, Line 615  Data::copyWithMask(const Data& other,
615      // and so I'm going to assume that you don't want your data objects getting a new shape.      // 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.");      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    
624    if ((selfrank>0) && (otherrank==0) &&(maskrank==0))    if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
625    {    {
# Line 709  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 719  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 742  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 773  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 801  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 816  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 828  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
# Line 941  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)));
   
        switch( dataPointRank ){  
             case 0 :  
                 numArray[0] = getDataAtOffset(offset);  
                 break;  
             case 1 :  
                 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;  
     }  
1046    }    }
1047    //    else
1048    // return the array    {
1049    return numArray;      // The old method would return an empty numArray of the given shape
1050        // I'm going to throw an exception because if we have zero points per sample we have problems
1051        throw DataException("Error - need at least 1 datapoint per sample.");
1052      }
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 1069  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 1083  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 1100  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 1247  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::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::numeric::array  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 1284  Data::integrateWorker() const Line 1245  Data::integrateWorker() const
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    dom->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
# Line 1932  Data::calc_minGlobalDataPoint(int& ProcN Line 1847  Data::calc_minGlobalDataPoint(int& ProcN
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 1950  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;      int error;
# Line 2021  Data::operator+=(const Data& right) Line 1936  Data::operator+=(const Data& right)
1936          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1937    }    }
1938    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=    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    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1941    return (*this);    return (*this);
1942  }  }
# Line 2032  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    MAKELAZYBINSELF(tmp,ADD)    (*this)+=tmp;
1952    binaryOp(tmp,plus<double>());    return *this;
   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 2052  Data::operator-=(const Data& right) Line 1975  Data::operator-=(const Data& right)
1975          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1976    }    }
1977    MAKELAZYBINSELF(right,SUB)    MAKELAZYBINSELF(right,SUB)
1978      exclusiveWrite();
1979    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1980    return (*this);    return (*this);
1981  }  }
# Line 2063  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    MAKELAZYBINSELF(tmp,SUB)    (*this)-=tmp;
   binaryOp(tmp,minus<double>());  
1991    return (*this);    return (*this);
1992  }  }
1993    
# Line 2075  Data::operator*=(const Data& right) Line 1998  Data::operator*=(const Data& right)
1998          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1999    }    }
2000    MAKELAZYBINSELF(right,MUL)    MAKELAZYBINSELF(right,MUL)
2001      exclusiveWrite();
2002    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
2003    return (*this);    return (*this);
2004  }  }
# Line 2086  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    MAKELAZYBINSELF(tmp,MUL)    (*this)*=tmp;
   binaryOp(tmp,multiplies<double>());  
2014    return (*this);    return (*this);
2015  }  }
2016    
# Line 2098  Data::operator/=(const Data& right) Line 2021  Data::operator/=(const Data& right)
2021          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2022    }    }
2023    MAKELAZYBINSELF(right,DIV)    MAKELAZYBINSELF(right,DIV)
2024      exclusiveWrite();
2025    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
2026    return (*this);    return (*this);
2027  }  }
# Line 2109  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    MAKELAZYBINSELF(tmp,DIV)    (*this)/=tmp;
   binaryOp(tmp,divides<double>());  
2037    return (*this);    return (*this);
2038  }  }
2039    
# Line 2289  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 2309  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 2361  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 2375  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,1);    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 2403  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 2532  escript::C_GeneralTensorProduct(Data& ar Line 2449  escript::C_GeneralTensorProduct(Data& ar
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 2555  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 2574  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 2604  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 2630  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 2647  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 2676  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 2705  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 2736  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 2765  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 2790  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 2820  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 2872  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 2938  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  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26