/[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 2644 by jfenwick, Wed Sep 2 04:14:03 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# 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 MAKELAZYOP2(X,Y,Z) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
65      {\
66        DataLazy* c=new DataLazy(borrowDataPtr(),X,Y,Z);\
67        return Data(c);\
68      }
69    
70    #define MAKELAZYBINSELF(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
71      {\
72        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
73    /*         m_data=c->getPtr();*/     set_m_data(c->getPtr());\
74        return (*this);\
75      }
76    
77    // like the above but returns a new data rather than *this
78    #define MAKELAZYBIN(R,X)   if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
79      {\
80        DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
81        return Data(c);\
82      }
83    
84    #define MAKELAZYBIN2(L,R,X)   if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
85      {\
86        DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
87        return Data(c);\
88      }
89    
90    // Do not use the following unless you want to make copies on assignment rather than
91    // share data.  CopyOnWrite should make this unnescessary.
92    // #define ASSIGNMENT_MEANS_DEEPCOPY
93    
94    namespace
95    {
96    
97    template <class ARR>
98    inline
99    boost::python::tuple
100    pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
101    {
102        using namespace boost::python;
103        using boost::python::tuple;
104        using boost::python::list;
105    
106        list l;
107        unsigned int dim=shape[0];
108        for (size_t i=0;i<dim;++i)
109        {
110        l.append(v[i+offset]);
111        }
112        return tuple(l);
113    }
114    
115    template <class ARR>
116    inline
117    boost::python::tuple
118    pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
119    {
120        using namespace boost::python;
121        using boost::python::tuple;
122        using boost::python::list;
123    
124        unsigned int shape0=shape[0];
125        unsigned int shape1=shape[1];
126        list lj;
127        for (size_t j=0;j<shape0;++j)
128        {
129            list li;
130        for (size_t i=0;i<shape1;++i)
131        {
132            li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
133        }
134        lj.append(tuple(li));
135        }
136        return tuple(lj);
137    }
138    
139    template <class ARR>
140    inline
141    boost::python::tuple
142    pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
143    {
144        using namespace boost::python;
145        using boost::python::tuple;
146        using boost::python::list;
147    
148        unsigned int shape0=shape[0];
149        unsigned int shape1=shape[1];
150        unsigned int shape2=shape[2];
151    
152        list lk;
153        for (size_t k=0;k<shape0;++k)
154        {
155            list lj;
156        for (size_t j=0;j<shape1;++j)
157        {
158            list li;
159            for (size_t i=0;i<shape2;++i)
160            {
161                    li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
162                }
163            lj.append(tuple(li));
164            }
165            lk.append(tuple(lj));
166        }
167        return tuple(lk);
168    }
169    
170    template <class ARR>
171    inline
172    boost::python::tuple
173    pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
174    {
175        using namespace boost::python;
176        using boost::python::tuple;
177        using boost::python::list;
178    
179        unsigned int shape0=shape[0];
180        unsigned int shape1=shape[1];
181        unsigned int shape2=shape[2];
182        unsigned int shape3=shape[3];
183    
184        list ll;
185        for (size_t l=0;l<shape0;++l)
186        {
187            list lk;
188        for (size_t k=0;k<shape1;++k)
189        {
190                list lj;
191                for (size_t j=0;j<shape2;++j)
192                {
193                    list li;
194                    for (size_t i=0;i<shape3;++i)
195                    {
196                        li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
197                    }
198                    lj.append(tuple(li));
199                }
200                lk.append(tuple(lj));
201        }
202            ll.append(tuple(lk));
203        }
204        return tuple(ll);
205    }
206    
207    
208    // This should be safer once the DataC RO changes have been brought in
209    template <class ARR>
210    boost::python::tuple
211    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
212    {
213       using namespace boost::python;
214       using boost::python::list;
215       int rank=shape.size();
216       if (rank==0)
217       {
218        return make_tuple(v[0]);
219       }
220       else if (rank==1)
221       {
222            return pointToTuple1(shape,v,0);
223       }
224       else if (rank==2)
225       {
226        return pointToTuple2(shape,v,0);
227       }
228       else if (rank==3)
229       {
230        return pointToTuple3(shape,v,0);
231       }
232       else if (rank==4)
233       {
234        return pointToTuple4(shape,v,0);
235       }
236       else
237         throw DataException("Unknown rank in pointToTuple.");
238    }
239    
240    }  // anonymous namespace
241    
242  Data::Data()  Data::Data()
243        : m_shared(false), m_lazy(false)
244  {  {
245    //    //
246    // Default data is type DataEmpty    // Default data is type DataEmpty
247    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
248    m_data=temp->getPtr();  //   m_data=temp->getPtr();
249      set_m_data(temp->getPtr());
250    m_protected=false;    m_protected=false;
251  }  }
252    
# Line 60  Data::Data(double value, Line 254  Data::Data(double value,
254             const tuple& shape,             const tuple& shape,
255             const FunctionSpace& what,             const FunctionSpace& what,
256             bool expanded)             bool expanded)
257        : m_shared(false), m_lazy(false)
258  {  {
259    DataTypes::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
260    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 271  Data::Data(double value,
271         const DataTypes::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
272         const FunctionSpace& what,         const FunctionSpace& what,
273             bool expanded)             bool expanded)
274        : m_shared(false), m_lazy(false)
275  {  {
276    int len = DataTypes::noValues(dataPointShape);    int len = DataTypes::noValues(dataPointShape);
   
277    DataVector temp_data(len,value,len);    DataVector temp_data(len,value,len);
 //   DataArrayView temp_dataView(temp_data, dataPointShape);  
   
 //   initialise(temp_dataView, what, expanded);  
278    initialise(temp_data, dataPointShape, what, expanded);    initialise(temp_data, dataPointShape, what, expanded);
   
279    m_protected=false;    m_protected=false;
280  }  }
281    
282  Data::Data(const Data& inData)  Data::Data(const Data& inData)
283        : m_shared(false), m_lazy(false)
284  {  {
285    m_data=inData.m_data;    set_m_data(inData.m_data);
286    m_protected=inData.isProtected();    m_protected=inData.isProtected();
287  }  }
288    
289    
290  Data::Data(const Data& inData,  Data::Data(const Data& inData,
291             const DataTypes::RegionType& region)             const DataTypes::RegionType& region)
292        : m_shared(false), m_lazy(false)
293  {  {
294    DataAbstract_ptr dat=inData.m_data;    DataAbstract_ptr dat=inData.m_data;
295    if (inData.isLazy())    if (inData.isLazy())
# Line 110  Data::Data(const Data& inData, Line 303  Data::Data(const Data& inData,
303    //    //
304    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
305    DataAbstract* tmp = dat->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
306    m_data=DataAbstract_ptr(tmp);    set_m_data(DataAbstract_ptr(tmp));
307    m_protected=false;    m_protected=false;
308    
309  }  }
310    
311  Data::Data(const Data& inData,  Data::Data(const Data& inData,
312             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
313        : m_shared(false), m_lazy(false)
314  {  {
315    if (inData.isEmpty())    if (inData.isEmpty())
316    {    {
317      throw DataException("Error - will not interpolate for instances of DataEmpty.");      throw DataException("Error - will not interpolate for instances of DataEmpty.");
318    }    }
319    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
320      m_data=inData.m_data;      set_m_data(inData.m_data);
321    }    }
322    else    else
323    {    {
# Line 130  Data::Data(const Data& inData, Line 325  Data::Data(const Data& inData,
325      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
326        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
327        {           // 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
328      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
329        }        }
330        // 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
331        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
332        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());          DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
333        m_data=DataAbstract_ptr(dc);  //       m_data=DataAbstract_ptr(dc);
334          set_m_data(DataAbstract_ptr(dc));
335      } else {      } else {
336        Data tmp(0,inData.getDataPointShape(),functionspace,true);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
337        // 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 345  Data::Data(const Data& inData,
345        } else {        } else {
346          inDataDomain->interpolateACross(tmp,inData);          inDataDomain->interpolateACross(tmp,inData);
347        }        }
348        m_data=tmp.m_data;  //       m_data=tmp.m_data;
349          set_m_data(tmp.m_data);
350      }      }
351    }    }
352    m_protected=false;    m_protected=false;
353  }  }
354    
355  Data::Data(DataAbstract* underlyingdata)  Data::Data(DataAbstract* underlyingdata)
356        : m_shared(false), m_lazy(false)
357  {  {
358  //  m_data=shared_ptr<DataAbstract>(underlyingdata);      set_m_data(underlyingdata->getPtr());
     m_data=underlyingdata->getPtr();  
359      m_protected=false;      m_protected=false;
360  }  }
361    
362  Data::Data(DataAbstract_ptr underlyingdata)  Data::Data(DataAbstract_ptr underlyingdata)
363        : m_shared(false), m_lazy(false)
364  {  {
365      m_data=underlyingdata;      set_m_data(underlyingdata);
366      m_protected=false;      m_protected=false;
367  }  }
368    
   
 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;  
 }*/  
   
369  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
370           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
371                   const FunctionSpace& what,                   const FunctionSpace& what,
372                   bool expanded)                   bool expanded)
373        : m_shared(false), m_lazy(false)
374  {  {
375     initialise(value,shape,what,expanded);     initialise(value,shape,what,expanded);
376     m_protected=false;     m_protected=false;
# Line 198  Data::Data(const DataTypes::ValueType& v Line 380  Data::Data(const DataTypes::ValueType& v
380  Data::Data(const object& value,  Data::Data(const object& value,
381         const FunctionSpace& what,         const FunctionSpace& what,
382             bool expanded)             bool expanded)
383        : m_shared(false), m_lazy(false)
384  {  {
385    numeric::array asNumArray(value);    WrappedArray w(value);
386    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
387    m_protected=false;    m_protected=false;
388  }  }
389    
390    
391  Data::Data(const object& value,  Data::Data(const object& value,
392             const Data& other)             const Data& other)
393        : m_shared(false), m_lazy(false)
394  {  {
395    numeric::array asNumArray(value);    WrappedArray w(value);
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);  
 // /*  for (int i=0; i < asNumArray.getrank(); i++) {  
 //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
 //   }*/  
 //   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 // /*  DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(asNumArray);*/  
 //   temp_data.copyFromNumArray(asNumArray);  
   
   //  
   // Create DataConstant using the given value and all other parameters  
   // copied from other. If value is a rank 0 object this Data  
   // will assume the point data shape of other.  
396    
397    if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {    // extract the shape of the array
398      const DataTypes::ShapeType& tempShape=w.getShape();
399      if (w.getRank()==0) {
400    
401    
402      // get the space for the data vector      // get the space for the data vector
403      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
404      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
405      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromArray(w,1);
406    
407      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
408    
409      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);  
   
410      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
411  //     boost::shared_ptr<DataAbstract> sp(t);  //     m_data=DataAbstract_ptr(t);
412  //     m_data=sp;      set_m_data(DataAbstract_ptr(t));
     m_data=DataAbstract_ptr(t);  
413    
414    } else {    } else {
415      //      //
416      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
417  //     initialise(temp_dataView, other.getFunctionSpace(), false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
418      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());  //     m_data=DataAbstract_ptr(t);
419  //     boost::shared_ptr<DataAbstract> sp(t);      set_m_data(DataAbstract_ptr(t));
 //     m_data=sp;  
     m_data=DataAbstract_ptr(t);  
420    }    }
421    m_protected=false;    m_protected=false;
422  }  }
423    
424  Data::~Data()  Data::~Data()
425  {  {
426      set_m_data(DataAbstract_ptr());
427  }  }
428    
429    
430    // only call in thread safe contexts.
431    // This method should be atomic
432    void Data::set_m_data(DataAbstract_ptr p)
433    {
434      if (m_data.get()!=0)  // release old ownership
435      {
436        m_data->removeOwner(this);
437      }
438      if (p.get()!=0)
439      {
440        m_data=p;
441        m_data->addOwner(this);
442        m_shared=m_data->isShared();
443        m_lazy=m_data->isLazy();
444      }
445    }
446    
447  void  void Data::initialise(const WrappedArray& value,
 Data::initialise(const boost::python::numeric::array& value,  
448                   const FunctionSpace& what,                   const FunctionSpace& what,
449                   bool expanded)                   bool expanded)
450  {  {
# Line 277  Data::initialise(const boost::python::nu Line 455  Data::initialise(const boost::python::nu
455    // within the shared_ptr constructor.    // within the shared_ptr constructor.
456    if (expanded) {    if (expanded) {
457      DataAbstract* temp=new DataExpanded(value, what);      DataAbstract* temp=new DataExpanded(value, what);
458  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
459  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
460    } else {    } else {
461      DataAbstract* temp=new DataConstant(value, what);      DataAbstract* temp=new DataConstant(value, what);
462  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
463  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
464    }    }
465  }  }
466    
# Line 302  Data::initialise(const DataTypes::ValueT Line 478  Data::initialise(const DataTypes::ValueT
478    // within the shared_ptr constructor.    // within the shared_ptr constructor.
479    if (expanded) {    if (expanded) {
480      DataAbstract* temp=new DataExpanded(what, shape, value);      DataAbstract* temp=new DataExpanded(what, shape, value);
481  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
482  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
483    } else {    } else {
484      DataAbstract* temp=new DataConstant(what, shape, value);      DataAbstract* temp=new DataConstant(what, shape, value);
485  //     boost::shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
486  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
487    }    }
488  }  }
489    
490    
 // 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";  
 //  }  
 // }  
   
491  escriptDataC  escriptDataC
492  Data::getDataC()  Data::getDataC()
493  {  {
# Line 410  Data::getDataC() const Line 504  Data::getDataC() const
504    return temp;    return temp;
505  }  }
506    
507    size_t
508    Data::getSampleBufferSize() const
509    {
510      return m_data->getSampleBufferSize();
511    }
512    
513    
514  const boost::python::tuple  const boost::python::tuple
515  Data::getShapeTuple() const  Data::getShapeTuple() const
516  {  {
# Line 435  Data::getShapeTuple() const Line 536  Data::getShapeTuple() const
536  // 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.
537  // 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
538  // 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
539  Data*  Data
540  Data::copySelf()  Data::copySelf()
541  {  {
542     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
543     return new Data(temp);     return Data(temp);
544  }  }
545    
546  void  void
# Line 447  Data::copy(const Data& other) Line 548  Data::copy(const Data& other)
548  {  {
549    DataAbstract* temp=other.m_data->deepCopy();    DataAbstract* temp=other.m_data->deepCopy();
550    DataAbstract_ptr p=temp->getPtr();    DataAbstract_ptr p=temp->getPtr();
551    m_data=p;  //   m_data=p;
552      set_m_data(p);
553  }  }
554    
555    
# Line 463  Data::delaySelf() Line 565  Data::delaySelf()
565  {  {
566    if (!isLazy())    if (!isLazy())
567    {    {
568      m_data=(new DataLazy(m_data))->getPtr();  //  m_data=(new DataLazy(m_data))->getPtr();
569        set_m_data((new DataLazy(m_data))->getPtr());
570    }    }
571  }  }
572    
573    
574    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
575    // to zero, all the tags from all the DataTags would be in the result.
576    // However since they all have the same value (0) whether they are there or not should not matter.
577    // So I have decided that for all types this method will create a constant 0.
578    // It can be promoted up as required.
579    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
580    // but we can deal with that if it arises.
581    //
582  void  void
583  Data::setToZero()  Data::setToZero()
584  {  {
# Line 474  Data::setToZero() Line 586  Data::setToZero()
586    {    {
587       throw DataException("Error - Operations not permitted on instances of DataEmpty.");       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
588    }    }
589    m_data->setToZero();    if (isLazy())
590      {
591         DataTypes::ValueType v(getNoValues(),0);
592         DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
593         DataLazy* dl=new DataLazy(dc->getPtr());
594         set_m_data(dl->getPtr());
595      }
596      else
597      {
598         exclusiveWrite();
599         m_data->setToZero();
600      }
601  }  }
602    
603    
604  void  void
605  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
606                     const Data& mask)                     const Data& mask)
# Line 538  Data::copyWithMask(const Data& other, Line 662  Data::copyWithMask(const Data& other,
662    {    {
663      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");      throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
664    }    }
665      unsigned int selfrank=getDataPointRank();
666      unsigned int otherrank=other2.getDataPointRank();
667      unsigned int maskrank=mask2.getDataPointRank();
668      if ((selfrank==0) && (otherrank>0 || maskrank>0))
669      {
670        // to get here we must be copying from a large object into a scalar
671        // I am not allowing this.
672        // If you are calling copyWithMask then you are considering keeping some existing values
673        // and so I'm going to assume that you don't want your data objects getting a new shape.
674        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
675      }
676      exclusiveWrite();
677    // Now we iterate over the elements    // Now we iterate over the elements
678    DataVector& self=getReadyPtr()->getVector();    DataVector& self=getReady()->getVectorRW();;
679    const DataVector& ovec=other2.getReadyPtr()->getVector();    const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
680    const DataVector& mvec=mask2.getReadyPtr()->getVector();    const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
681    if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))  
682      if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
683    {    {
684      throw DataException("Error - size mismatch in arguments to copyWithMask.");      // Not allowing this combination.
685        // it is not clear what the rank of the target should be.
686        // Should it be filled with the scalar (rank stays the same);
687        // or should the target object be reshaped to be a scalar as well.
688        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
689      }
690      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
691      {
692        if (mvec[0]>0)      // copy whole object if scalar is >0
693        {
694            copy(other);
695        }
696        return;
697      }
698      if (isTagged())       // so all objects involved will also be tagged
699      {
700        // note the !
701        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
702            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
703        {
704            throw DataException("copyWithMask, shape mismatch.");
705        }
706    
707        // We need to consider the possibility that tags are missing or in the wrong order
708        // My guiding assumption here is: All tagged Datas are assumed to have the default value for
709        // all tags which are not explicitly defined
710    
711        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
712        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
713        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
714    
715        // first, add any tags missing from other or mask
716        const DataTagged::DataMapType& olookup=optr->getTagLookup();
717            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
718        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
719        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
720        for (i=olookup.begin();i!=olookup.end();i++)
721        {
722               tptr->addTag(i->first);
723            }
724            for (i=mlookup.begin();i!=mlookup.end();i++) {
725               tptr->addTag(i->first);
726            }
727        // now we know that *this has all the required tags but they aren't guaranteed to be in
728        // the same order
729    
730        // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
731        if ((selfrank==otherrank) && (otherrank==maskrank))
732        {
733            for (i=tlookup.begin();i!=tlookup.end();i++)
734            {
735                // get the target offset
736                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
737                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
738                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
739                for (int j=0;j<getDataPointSize();++j)
740                {
741                    if (mvec[j+moff]>0)
742                    {
743                        self[j+toff]=ovec[j+ooff];
744                    }
745                }
746                }
747            // now for the default value
748            for (int j=0;j<getDataPointSize();++j)
749            {
750                if (mvec[j+mptr->getDefaultOffset()]>0)
751                {
752                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
753                }
754            }
755        }
756        else    // other is a scalar
757        {
758            for (i=tlookup.begin();i!=tlookup.end();i++)
759            {
760                // get the target offset
761                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
762                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
763                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
764                for (int j=0;j<getDataPointSize();++j)
765                {
766                    if (mvec[j+moff]>0)
767                    {
768                        self[j+toff]=ovec[ooff];
769                    }
770                }
771                }
772            // now for the default value
773            for (int j=0;j<getDataPointSize();++j)
774            {
775                if (mvec[j+mptr->getDefaultOffset()]>0)
776                {
777                    self[j+tptr->getDefaultOffset()]=ovec[0];
778                }
779            }
780        }
781    
782        return;         // ugly
783      }
784      // mixed scalar and non-scalar operation
785      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
786      {
787            size_t num_points=self.size();
788        // OPENMP 3.0 allows unsigned loop vars.
789        #if defined(_OPENMP) && (_OPENMP < 200805)
790        long i;
791        #else
792        size_t i;
793        #endif  
794        size_t psize=getDataPointSize();    
795        #pragma omp parallel for private(i) schedule(static)
796        for (i=0;i<num_points;++i)
797        {
798            if (mvec[i]>0)
799            {
800                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
801            }                   // dest point
802        }
803        return;         // ugly!
804      }
805      // tagged data is already taken care of so we only need to worry about shapes
806      // special cases with scalars are already dealt with so all we need to worry about is shape
807      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
808      {
809        ostringstream oss;
810        oss <<"Error - size mismatch in arguments to copyWithMask.";
811        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
812        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
813        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
814        throw DataException(oss.str());
815    }    }
816    size_t num_points=self.size();    size_t num_points=self.size();
817    
# Line 564  Data::copyWithMask(const Data& other, Line 831  Data::copyWithMask(const Data& other,
831    }    }
832  }  }
833    
   
   
834  bool  bool
835  Data::isExpanded() const  Data::isExpanded() const
836  {  {
# Line 574  Data::isExpanded() const Line 839  Data::isExpanded() const
839  }  }
840    
841  bool  bool
842    Data::actsExpanded() const
843    {
844      return m_data->actsExpanded();
845    
846    }
847    
848    
849    bool
850  Data::isTagged() const  Data::isTagged() const
851  {  {
852    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
# Line 597  Data::isConstant() const Line 870  Data::isConstant() const
870  bool  bool
871  Data::isLazy() const  Data::isLazy() const
872  {  {
873    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
874  }  }
875    
876  // 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 901  Data::expand()
901    if (isConstant()) {    if (isConstant()) {
902      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
903      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
904  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
905  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
906    } else if (isTagged()) {    } else if (isTagged()) {
907      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
908      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
909  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
910  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
911    } else if (isExpanded()) {    } else if (isExpanded()) {
912      //      //
913      // do nothing      // do nothing
# Line 656  Data::tag() Line 927  Data::tag()
927    if (isConstant()) {    if (isConstant()) {
928      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
929      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
930  //     shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
931  //     m_data=temp_data;      set_m_data(temp->getPtr());
     m_data=temp->getPtr();  
932    } else if (isTagged()) {    } else if (isTagged()) {
933      // do nothing      // do nothing
934    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 671  Data::tag() Line 941  Data::tag()
941       {       {
942      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");      throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
943       }       }
944       m_data=res;      //      m_data=res;
945         set_m_data(res);
946       tag();       tag();
947    } else {    } else {
948      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 954  Data::resolve()
954  {  {
955    if (isLazy())    if (isLazy())
956    {    {
957       m_data=m_data->resolve();  //      m_data=m_data->resolve();
958        set_m_data(m_data->resolve());
959    }    }
960  }  }
961    
962    void
963    Data::requireWrite()
964    {
965      resolve();
966      exclusiveWrite();
967    }
968    
969  Data  Data
970  Data::oneOver() const  Data::oneOver() const
971  {  {
972    if (isLazy())    MAKELAZYOP(RECIP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);  
     return Data(c);  
   }  
973    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));    return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
974  }  }
975    
976  Data  Data
977  Data::wherePositive() const  Data::wherePositive() const
978  {  {
979    if (isLazy())    MAKELAZYOP(GZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GZ);  
     return Data(c);  
   }  
980    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
981  }  }
982    
983  Data  Data
984  Data::whereNegative() const  Data::whereNegative() const
985  {  {
986    if (isLazy())    MAKELAZYOP(LZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LZ);  
     return Data(c);  
   }  
987    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
988  }  }
989    
990  Data  Data
991  Data::whereNonNegative() const  Data::whereNonNegative() const
992  {  {
993    if (isLazy())    MAKELAZYOP(GEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);  
     return Data(c);  
   }  
994    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
995  }  }
996    
997  Data  Data
998  Data::whereNonPositive() const  Data::whereNonPositive() const
999  {  {
1000    if (isLazy())    MAKELAZYOP(LEZ)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);  
     return Data(c);  
   }  
1001    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1002  }  }
1003    
1004  Data  Data
1005  Data::whereZero(double tol) const  Data::whereZero(double tol) const
1006  {  {
1007    Data dataAbs=abs();  //   Data dataAbs=abs();
1008    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1009       MAKELAZYOPOFF(EZ,tol)
1010       return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1011    
1012  }  }
1013    
1014  Data  Data
1015  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
1016  {  {
1017    Data dataAbs=abs();  //   Data dataAbs=abs();
1018    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1019      MAKELAZYOPOFF(NEZ,tol)
1020      return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1021    
1022  }  }
1023    
1024  Data  Data
# Line 767  bool Line 1031  bool
1031  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
1032  {  {
1033    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());  
 //     }  
 //   }  
1034  }  }
1035    
1036  Data  Data
# Line 813  Data::getDataPointSize() const Line 1067  Data::getDataPointSize() const
1067    return m_data->getNoValues();    return m_data->getNoValues();
1068  }  }
1069    
1070    
1071  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1072  Data::getLength() const  Data::getLength() const
1073  {  {
1074    return m_data->getLength();    return m_data->getLength();
1075  }  }
1076    
 const  
 boost::python::numeric::array  
 Data:: getValueOfDataPoint(int dataPointNo)  
 {  
   int i, j, k, l;  
   
   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);  
1077    
1078    //  // There is no parallelism here ... elements need to be added in the correct order.
1079    // the shape of the returned numeric array will be the same  //   If we could presize the list and then fill in the elements it might work
1080    // as that of the data point  //   This would need setting elements to be threadsafe.
1081    int arrayRank = dataPointRank;  //   Having mulitple C threads calling into one interpreter is aparently a no-no.
1082    const DataTypes::ShapeType& arrayShape = dataPointShape;  const boost::python::object
1083    Data::toListOfTuples(bool scalarastuple)
1084    {
1085        using namespace boost::python;
1086        using boost::python::list;
1087        if (get_MPISize()>1)
1088        {
1089            throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1090        }
1091        unsigned int rank=getDataPointRank();
1092        unsigned int size=getDataPointSize();
1093    
1094    //      int npoints=getNumDataPoints();
1095    // resize the numeric array to the shape just calculated      expand();           // This will also resolve if required
1096    if (arrayRank==0) {      const DataTypes::ValueType& vec=getReady()->getVectorRO();
1097      numArray.resize(1);      boost::python::list temp;
1098    }      temp.append(object());
1099    if (arrayRank==1) {      boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints"  trick
1100      numArray.resize(arrayShape[0]);      if (rank==0)
1101    }      {
1102    if (arrayRank==2) {          long count;
1103      numArray.resize(arrayShape[0],arrayShape[1]);          if (scalarastuple)
1104    }          {
1105    if (arrayRank==3) {              for (count=0;count<npoints;++count)
1106      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);              {
1107    }          res[count]=make_tuple(vec[count]);
1108    if (arrayRank==4) {              }
1109      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);          }
1110    }          else
1111            {
1112                for (count=0;count<npoints;++count)
1113                {
1114                    res[count]=vec[count];
1115                }
1116            }
1117        }
1118        else if (rank==1)
1119        {
1120            size_t count;
1121            size_t offset=0;
1122            for (count=0;count<npoints;++count,offset+=size)
1123            {
1124                res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1125            }
1126        }
1127        else if (rank==2)
1128        {
1129            size_t count;
1130            size_t offset=0;
1131            for (count=0;count<npoints;++count,offset+=size)
1132            {
1133            res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1134            }
1135        }
1136        else if (rank==3)
1137        {
1138            size_t count;
1139            size_t offset=0;
1140            for (count=0;count<npoints;++count,offset+=size)
1141            {
1142                res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1143            }
1144        }
1145        else if (rank==4)
1146        {
1147            size_t count;
1148            size_t offset=0;
1149            for (count=0;count<npoints;++count,offset+=size)
1150            {
1151                res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1152            }
1153        }
1154        else
1155        {
1156            throw DataException("Unknown rank in ::toListOfTuples()");
1157        }
1158        return res;
1159    }
1160    
1161    if (getNumDataPointsPerSample()>0) {  const boost::python::object
1162    Data::getValueOfDataPointAsTuple(int dataPointNo)
1163    {
1164       forceResolve();
1165       if (getNumDataPointsPerSample()>0) {
1166         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1167         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1168         //         //
1169         // Check a valid sample number has been supplied         // Check a valid sample number has been supplied
1170         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1171             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1172         }         }
1173    
1174         //         //
1175         // Check a valid data point number has been supplied         // Check a valid data point number has been supplied
1176         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1177             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");             throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1178         }         }
1179         // TODO: global error handling         // TODO: global error handling
        // create a view of the data if it is stored locally  
 //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);  
1180         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);         DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1181           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1182      }
1183         switch( dataPointRank ){    else
1184              case 0 :    {
1185                  numArray[0] = getDataAtOffset(offset);      // The pre-numpy method would return an empty array of the given shape
1186                  break;      // I'm going to throw an exception because if we have zero points per sample we have problems
1187              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;  
     }  
1188    }    }
   //  
   // return the array  
   return numArray;  
1189    
1190  }  }
1191    
1192    
1193  void  void
1194  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1195  {  {
1196      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
1197      boost::python::numeric::array num_array(py_object);      setValueOfDataPointToArray(dataPointNo,py_object);
     setValueOfDataPointToArray(dataPointNo,num_array);  
1198  }  }
1199    
1200  void  void
1201  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1202  {  {
1203    if (isProtected()) {    if (isProtected()) {
1204          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1205    }    }
1206    FORCERESOLVE;    forceResolve();
1207    
1208      WrappedArray w(obj);
1209    //    //
1210    // check rank    // check rank
1211    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1212        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1213    
1214    //    //
1215    // check shape of num_array    // check shape of array
1216    for (unsigned int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1217      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1218         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1219    }    }
1220    //    //
1221    // make sure data is expanded:    // make sure data is expanded:
# Line 949  Data::setValueOfDataPointToArray(int dat Line 1226  Data::setValueOfDataPointToArray(int dat
1226    if (getNumDataPointsPerSample()>0) {    if (getNumDataPointsPerSample()>0) {
1227         int sampleNo = dataPointNo/getNumDataPointsPerSample();         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1228         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1229         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1230    } else {    } else {
1231         m_data->copyToDataPoint(-1, 0,num_array);         m_data->copyToDataPoint(-1, 0,w);
1232    }    }
1233  }  }
1234    
# Line 963  Data::setValueOfDataPoint(int dataPointN Line 1240  Data::setValueOfDataPoint(int dataPointN
1240    }    }
1241    //    //
1242    // make sure data is expanded:    // make sure data is expanded:
1243    FORCERESOLVE;    forceResolve();
1244    if (!isExpanded()) {    if (!isExpanded()) {
1245      expand();      expand();
1246    }    }
# Line 977  Data::setValueOfDataPoint(int dataPointN Line 1254  Data::setValueOfDataPoint(int dataPointN
1254  }  }
1255    
1256  const  const
1257  boost::python::numeric::array  boost::python::object
1258  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1259  {  {
1260    size_t length=0;    // This could be lazier than it is now
1261    int i, j, k, l, pos;    forceResolve();
   FORCERESOLVE;  
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   const DataTypes::ShapeType& dataPointShape = getDataPointShape();  
1262    
1263    //    // copy datapoint into a buffer
1264    // create the numeric array to be returned    // broadcast buffer to all nodes
1265    boost::python::numeric::array numArray(0.0);    // convert buffer to tuple
1266      // return tuple
1267    
1268    //    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1269    // the shape of the returned numeric array will be the same    size_t length=DataTypes::noValues(dataPointShape);
   // 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]);  
   }  
1270    
1271    // added for the MPI communication    // added for the MPI communication
   length=1;  
   for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];  
1272    double *tmpData = new double[length];    double *tmpData = new double[length];
1273    
1274    //    // updated for the MPI case
1275    // load the values for the data point into the numeric array.    if( get_MPIRank()==procNo ){
1276          if (getNumDataPointsPerSample()>0) {
1277      // updated for the MPI case      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1278      if( get_MPIRank()==procNo ){      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1279               if (getNumDataPointsPerSample()>0) {      //
1280                  int sampleNo = dataPointNo/getNumDataPointsPerSample();      // Check a valid sample number has been supplied
1281                  int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1282                  //          throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: 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;  
         }  
             }  
1283      }      }
1284          #ifdef PASO_MPI  
1285          // broadcast the data to all other processes      //
1286      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      // Check a valid data point number has been supplied
1287          #endif      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1288            throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: 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;  
1289      }      }
1290        // TODO: global error handling
1291        DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1292    
1293      delete [] tmpData;      memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1294         }
1295      }
1296    #ifdef PASO_MPI
1297      // broadcast the data to all other processes
1298      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1299    #endif
1300    
1301      boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1302      delete [] tmpData;
1303    //    //
1304    // return the loaded array    // return the loaded array
1305    return numArray;    return t;
1306    
1307  }  }
1308    
1309    
1310  boost::python::numeric::array  boost::python::object
1311  Data::integrate_const() const  Data::integrateToTuple_const() const
1312  {  {
1313    if (isLazy())    if (isLazy())
1314    {    {
# Line 1127  Data::integrate_const() const Line 1317  Data::integrate_const() const
1317    return integrateWorker();    return integrateWorker();
1318  }  }
1319    
1320  boost::python::numeric::array  boost::python::object
1321  Data::integrate()  Data::integrateToTuple()
1322  {  {
1323    if (isLazy())    if (isLazy())
1324    {    {
1325      expand();      expand();
1326    }    }
1327    return integrateWorker();    return integrateWorker();
 }  
   
1328    
1329    }
1330    
1331  boost::python::numeric::array  boost::python::object
1332  Data::integrateWorker() const  Data::integrateWorker() const
1333  {  {
   int index;  
   int rank = getDataPointRank();  
1334    DataTypes::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1335    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1336    
# Line 1151  Data::integrateWorker() const Line 1338  Data::integrateWorker() const
1338    // calculate the integral values    // calculate the integral values
1339    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1340    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1341      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1342      if (dom==0)
1343      {            
1344        throw DataException("Can not integrate over non-continuous domains.");
1345      }
1346  #ifdef PASO_MPI  #ifdef PASO_MPI
1347    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1348    // 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
1349    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1350    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
1351    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]; }
1352    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 );
1353    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }    for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1354      tuple result=pointToTuple(shape,tmp);
1355    delete[] tmp;    delete[] tmp;
1356    delete[] tmp_local;    delete[] tmp_local;
1357  #else  #else
1358    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1359    /*  double *tmp = new double[dataPointSize];
1360      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1361      tuple result=pointToTuple(shape,integrals);
1362    //   delete tmp;
1363  #endif  #endif
1364    
   //  
   // 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];  
           }  
         }  
       }  
     }  
   }  
1365    
1366    //    return result;
   // return the loaded array  
   return bp_array;  
1367  }  }
1368    
1369  Data  Data
1370  Data::sin() const  Data::sin() const
1371  {  {
1372    if (isLazy())    MAKELAZYOP(SIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIN);  
     return Data(c);  
   }  
1373    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1374  }  }
1375    
1376  Data  Data
1377  Data::cos() const  Data::cos() const
1378  {  {
1379    if (isLazy())    MAKELAZYOP(COS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COS);  
     return Data(c);  
   }  
1380    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1381  }  }
1382    
1383  Data  Data
1384  Data::tan() const  Data::tan() const
1385  {  {
1386    if (isLazy())    MAKELAZYOP(TAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TAN);  
     return Data(c);  
   }  
1387    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1388  }  }
1389    
1390  Data  Data
1391  Data::asin() const  Data::asin() const
1392  {  {
1393    if (isLazy())    MAKELAZYOP(ASIN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);  
     return Data(c);  
   }  
1394    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1395  }  }
1396    
1397  Data  Data
1398  Data::acos() const  Data::acos() const
1399  {  {
1400    if (isLazy())    MAKELAZYOP(ACOS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);  
     return Data(c);  
   }  
1401    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1402  }  }
1403    
# Line 1279  Data::acos() const Line 1405  Data::acos() const
1405  Data  Data
1406  Data::atan() const  Data::atan() const
1407  {  {
1408    if (isLazy())    MAKELAZYOP(ATAN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);  
     return Data(c);  
   }  
1409    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1410  }  }
1411    
1412  Data  Data
1413  Data::sinh() const  Data::sinh() const
1414  {  {
1415    if (isLazy())    MAKELAZYOP(SINH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SINH);  
     return Data(c);  
   }  
1416    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1417  }  }
1418    
1419  Data  Data
1420  Data::cosh() const  Data::cosh() const
1421  {  {
1422    if (isLazy())    MAKELAZYOP(COSH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),COSH);  
     return Data(c);  
   }  
1423    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1424  }  }
1425    
1426  Data  Data
1427  Data::tanh() const  Data::tanh() const
1428  {  {
1429    if (isLazy())    MAKELAZYOP(TANH)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),TANH);  
     return Data(c);  
   }  
1430    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1431  }  }
1432    
# Line 1324  Data::tanh() const Line 1434  Data::tanh() const
1434  Data  Data
1435  Data::erf() const  Data::erf() const
1436  {  {
1437  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1438    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1439  #else  #else
1440    if (isLazy())    MAKELAZYOP(ERF)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ERF);  
     return Data(c);  
   }  
1441    return C_TensorUnaryOperation(*this, ::erf);    return C_TensorUnaryOperation(*this, ::erf);
1442  #endif  #endif
1443  }  }
# Line 1339  Data::erf() const Line 1445  Data::erf() const
1445  Data  Data
1446  Data::asinh() const  Data::asinh() const
1447  {  {
1448    if (isLazy())    MAKELAZYOP(ASINH)
1449    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1450    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1451  #else  #else
1452    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1354  Data::asinh() const Line 1456  Data::asinh() const
1456  Data  Data
1457  Data::acosh() const  Data::acosh() const
1458  {  {
1459    if (isLazy())    MAKELAZYOP(ACOSH)
1460    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1461    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1462  #else  #else
1463    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1369  Data::acosh() const Line 1467  Data::acosh() const
1467  Data  Data
1468  Data::atanh() const  Data::atanh() const
1469  {  {
1470    if (isLazy())    MAKELAZYOP(ATANH)
1471    {  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
     DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);  
     return Data(c);  
   }  
 #ifdef _WIN32  
1472    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1473  #else  #else
1474    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1383  Data::atanh() const Line 1477  Data::atanh() const
1477    
1478  Data  Data
1479  Data::log10() const  Data::log10() const
1480  {  if (isLazy())  {
1481    {    MAKELAZYOP(LOG10)
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);  
     return Data(c);  
   }  
1482    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1483  }  }
1484    
1485  Data  Data
1486  Data::log() const  Data::log() const
1487  {  {
1488    if (isLazy())    MAKELAZYOP(LOG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),LOG);  
     return Data(c);  
   }  
1489    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1490  }  }
1491    
1492  Data  Data
1493  Data::sign() const  Data::sign() const
1494  {  {
1495    if (isLazy())    MAKELAZYOP(SIGN)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);  
     return Data(c);  
   }  
1496    return C_TensorUnaryOperation(*this, escript::fsign);    return C_TensorUnaryOperation(*this, escript::fsign);
1497  }  }
1498    
1499  Data  Data
1500  Data::abs() const  Data::abs() const
1501  {  {
1502    if (isLazy())    MAKELAZYOP(ABS)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),ABS);  
     return Data(c);  
   }  
1503    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1504  }  }
1505    
1506  Data  Data
1507  Data::neg() const  Data::neg() const
1508  {  {
1509    if (isLazy())    MAKELAZYOP(NEG)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NEG);  
     return Data(c);  
   }  
1510    return C_TensorUnaryOperation(*this, negate<double>());    return C_TensorUnaryOperation(*this, negate<double>());
1511  }  }
1512    
# Line 1448  Data::pos() const Line 1523  Data::pos() const
1523    
1524  Data  Data
1525  Data::exp() const  Data::exp() const
1526  {    {
1527    if (isLazy())    MAKELAZYOP(EXP)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),EXP);  
     return Data(c);  
   }  
1528    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1529  }  }
1530    
1531  Data  Data
1532  Data::sqrt() const  Data::sqrt() const
1533  {  {
1534    if (isLazy())    MAKELAZYOP(SQRT)
   {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);  
     return Data(c);  
   }  
1535    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1536  }  }
1537    
# Line 1483  Data::Lsup() Line 1550  Data::Lsup()
1550  {  {
1551     if (isLazy())     if (isLazy())
1552     {     {
1553      expand();      resolve();
1554     }     }
1555     return LsupWorker();     return LsupWorker();
1556  }  }
# Line 1503  Data::sup() Line 1570  Data::sup()
1570  {  {
1571     if (isLazy())     if (isLazy())
1572     {     {
1573      expand();      resolve();
1574     }     }
1575     return supWorker();     return supWorker();
1576  }  }
# Line 1523  Data::inf() Line 1590  Data::inf()
1590  {  {
1591     if (isLazy())     if (isLazy())
1592     {     {
1593      expand();      resolve();
1594     }     }
1595     return infWorker();     return infWorker();
1596  }  }
# Line 1633  Data::swapaxes(const int axis0, const in Line 1700  Data::swapaxes(const int axis0, const in
1700       if (axis0 == axis1) {       if (axis0 == axis1) {
1701           throw DataException("Error - Data::swapaxes: axis indices must be different.");           throw DataException("Error - Data::swapaxes: axis indices must be different.");
1702       }       }
1703       if (axis0 > axis1) {       MAKELAZYOP2(SWAP,axis0,axis1)
1704           axis0_tmp=axis1;       if (axis0 > axis1)
1705           axis1_tmp=axis0;       {
1706       } else {      axis0_tmp=axis1;
1707           axis0_tmp=axis0;      axis1_tmp=axis0;
          axis1_tmp=axis1;  
1708       }       }
1709       for (int i=0; i<rank; i++) {       else
1710         if (i == axis0_tmp) {       {
1711            ev_shape.push_back(s[axis1_tmp]);      axis0_tmp=axis0;
1712         } else if (i == axis1_tmp) {      axis1_tmp=axis1;
1713            ev_shape.push_back(s[axis0_tmp]);       }
1714         } else {       for (int i=0; i<rank; i++)
1715            ev_shape.push_back(s[i]);       {
1716         }      if (i == axis0_tmp)
1717        {
1718            ev_shape.push_back(s[axis1_tmp]);
1719        }
1720        else if (i == axis1_tmp)
1721        {
1722            ev_shape.push_back(s[axis0_tmp]);
1723        }
1724        else
1725        {
1726            ev_shape.push_back(s[i]);
1727        }
1728       }       }
1729       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1730       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1731       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);       m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1732       return ev;       return ev;
   
1733  }  }
1734    
1735  Data  Data
# Line 1672  Data::symmetric() const Line 1748  Data::symmetric() const
1748       else {       else {
1749          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.");
1750       }       }
1751       if (isLazy())       MAKELAZYOP(SYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),SYM);  
     return Data(c);  
      }  
1752       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1753       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1754       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1686  Data::symmetric() const Line 1758  Data::symmetric() const
1758  Data  Data
1759  Data::nonsymmetric() const  Data::nonsymmetric() const
1760  {  {
1761       if (isLazy())       MAKELAZYOP(NSYM)
      {  
     DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);  
     return Data(c);  
      }  
1762       // check input       // check input
1763       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1764       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1722  Data::nonsymmetric() const Line 1790  Data::nonsymmetric() const
1790       }       }
1791  }  }
1792    
   
 // 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).  
1793  Data  Data
1794  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1795  {  {    
1796       if (isLazy())       MAKELAZYOPOFF(TRACE,axis_offset)
1797         if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1798       {       {
1799      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);  
1800       }       }
1801       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1802       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1787  Data::trace(int axis_offset) const Line 1848  Data::trace(int axis_offset) const
1848  Data  Data
1849  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1850  {      {    
1851       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);  
      }  
1852       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1853       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
1854       // 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 1950  Data::calc_minGlobalDataPoint(int& ProcN
1950    double next,local_min;    double next,local_min;
1951    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1952    
1953    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1954    {    {
1955      local_min=min;      local_min=min;
1956      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1957      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1958        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1959          next=temp.getDataAtOffset(temp.getDataOffset(i,j));          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1960          if (next<local_min) {          if (next<local_min) {
1961            local_min=next;            local_min=next;
1962            local_lowi=i;            local_lowi=i;
# Line 1909  Data::calc_minGlobalDataPoint(int& ProcN Line 1965  Data::calc_minGlobalDataPoint(int& ProcN
1965        }        }
1966      }      }
1967      #pragma omp critical      #pragma omp critical
1968      if (local_min<min) {      if (local_min<min) {    // If we found a smaller value than our sentinel
1969        min=local_min;        min=local_min;
1970        lowi=local_lowi;        lowi=local_lowi;
1971        lowj=local_lowj;        lowj=local_lowj;
# Line 1917  Data::calc_minGlobalDataPoint(int& ProcN Line 1973  Data::calc_minGlobalDataPoint(int& ProcN
1973    }    }
1974    
1975  #ifdef PASO_MPI  #ifdef PASO_MPI
1976      // determine the processor on which the minimum occurs    // determine the processor on which the minimum occurs
1977      next = temp.getDataPoint(lowi,lowj);    next = temp.getDataPointRO(lowi,lowj);
1978      int lowProc = 0;    int lowProc = 0;
1979      double *globalMins = new double[get_MPISize()+1];    double *globalMins = new double[get_MPISize()+1];
1980      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );    int error;
1981      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1982      if( get_MPIRank()==0 ){  
1983          next = globalMins[lowProc];    if( get_MPIRank()==0 ){
1984          for( i=1; i<get_MPISize(); i++ )      next = globalMins[lowProc];
1985              if( next>globalMins[i] ){      for( i=1; i<get_MPISize(); i++ )
1986                  lowProc = i;          if( next>globalMins[i] ){
1987                  next = globalMins[i];              lowProc = i;
1988              }              next = globalMins[i];
1989      }          }
1990      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );    }
1991      MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
1992    
1993      delete [] globalMins;    delete [] globalMins;
1994      ProcNo = lowProc;    ProcNo = lowProc;
1995  #else  #else
1996      ProcNo = 0;    ProcNo = 0;
1997  #endif  #endif
1998    DataPointNo = lowj + lowi * numDPPSample;    DataPointNo = lowj + lowi * numDPPSample;
1999  }  }
2000    
2001    
2002    const boost::python::tuple
2003    Data::maxGlobalDataPoint() const
2004    {
2005      int DataPointNo;
2006      int ProcNo;
2007      calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2008      return make_tuple(ProcNo,DataPointNo);
2009    }
2010    
2011    void
2012    Data::calc_maxGlobalDataPoint(int& ProcNo,
2013                            int& DataPointNo) const
2014    {
2015      if (isLazy())
2016      {
2017        Data temp(*this);   // to get around the fact that you can't resolve a const Data
2018        temp.resolve();
2019        return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2020      }
2021      int i,j;
2022      int highi=0,highj=0;
2023    //-------------
2024      double max=numeric_limits<double>::min();
2025    
2026      Data temp=maxval();
2027    
2028      int numSamples=temp.getNumSamples();
2029      int numDPPSample=temp.getNumDataPointsPerSample();
2030    
2031      double next,local_max;
2032      int local_highi=0,local_highj=0;  
2033    
2034      #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2035      {
2036        local_max=max;
2037        #pragma omp for private(i,j) schedule(static)
2038        for (i=0; i<numSamples; i++) {
2039          for (j=0; j<numDPPSample; j++) {
2040            next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2041            if (next>local_max) {
2042              local_max=next;
2043              local_highi=i;
2044              local_highj=j;
2045            }
2046          }
2047        }
2048        #pragma omp critical
2049        if (local_max>max) {    // If we found a larger value than our sentinel
2050          max=local_max;
2051          highi=local_highi;
2052          highj=local_highj;
2053        }
2054      }
2055    
2056    #ifdef PASO_MPI
2057      // determine the processor on which the maximum occurs
2058      next = temp.getDataPointRO(highi,highj);
2059      int highProc = 0;
2060      double *globalMaxs = new double[get_MPISize()+1];
2061      int error;
2062      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2063    
2064      if( get_MPIRank()==0 ){
2065      next = globalMaxs[highProc];
2066      for( i=1; i<get_MPISize(); i++ )
2067        if( next>globalMaxs[i] ){
2068            highProc = i;
2069            next = globalMaxs[i];
2070        }
2071      }
2072      MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2073      delete [] globalMaxs;
2074      ProcNo = highProc;
2075    #else
2076      ProcNo = 0;
2077    #endif
2078      DataPointNo = highj + highi * numDPPSample;
2079    }
2080    
2081  void  void
2082  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
2083  {  {
# Line 1977  Data::saveVTK(std::string fileName) cons Line 2114  Data::saveVTK(std::string fileName) cons
2114    }    }
2115    boost::python::dict args;    boost::python::dict args;
2116    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2117    getDomain()->saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
2118    return;    return;
2119  }  }
2120    
2121    
2122    
2123  Data&  Data&
2124  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
2125  {  {
2126    if (isProtected()) {    if (isProtected()) {
2127          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2128    }    }
2129    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2130    {    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
2131      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=    binaryOp(right,plus<double>());
2132          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,plus<double>());  
     return (*this);  
   }  
2133  }  }
2134    
2135  Data&  Data&
# Line 2007  Data::operator+=(const boost::python::ob Line 2139  Data::operator+=(const boost::python::ob
2139          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2140    }    }
2141    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2142    if (isLazy())    (*this)+=tmp;
2143    {    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);  
   }  
2144  }  }
2145    
2146  // Hmmm, operator= makes a deep copy but the copy constructor does not?  // Hmmm, operator= makes a deep copy but the copy constructor does not?
2147  Data&  Data&
2148  Data::operator=(const Data& other)  Data::operator=(const Data& other)
2149  {  {
2150    #if defined ASSIGNMENT_MEANS_DEEPCOPY  
2151    // This should not be used.
2152    // Just leaving this here until I have completed transition
2153    copy(other);    copy(other);
2154    #else
2155      m_protected=false;        // since any changes should be caught by exclusiveWrite();
2156    //   m_data=other.m_data;
2157      set_m_data(other.m_data);
2158    #endif
2159    return (*this);    return (*this);
2160  }  }
2161    
# Line 2034  Data::operator-=(const Data& right) Line 2165  Data::operator-=(const Data& right)
2165    if (isProtected()) {    if (isProtected()) {
2166          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2167    }    }
2168    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,SUB)
2169    {    exclusiveWrite();
2170      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=    binaryOp(right,minus<double>());
2171          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,minus<double>());  
     return (*this);  
   }  
2172  }  }
2173    
2174  Data&  Data&
# Line 2054  Data::operator-=(const boost::python::ob Line 2178  Data::operator-=(const boost::python::ob
2178          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2179    }    }
2180    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2181    if (isLazy())    (*this)-=tmp;
2182    {    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);  
   }  
2183  }  }
2184    
2185  Data&  Data&
# Line 2073  Data::operator*=(const Data& right) Line 2188  Data::operator*=(const Data& right)
2188    if (isProtected()) {    if (isProtected()) {
2189          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2190    }    }
2191    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,MUL)
2192    {    exclusiveWrite();
2193      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=    binaryOp(right,multiplies<double>());
2194          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,multiplies<double>());  
     return (*this);  
   }  
2195  }  }
2196    
2197  Data&  Data&
# Line 2093  Data::operator*=(const boost::python::ob Line 2201  Data::operator*=(const boost::python::ob
2201          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2202    }    }
2203    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2204    if (isLazy())    (*this)*=tmp;
2205    {    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);  
   }  
2206  }  }
2207    
2208  Data&  Data&
# Line 2112  Data::operator/=(const Data& right) Line 2211  Data::operator/=(const Data& right)
2211    if (isProtected()) {    if (isProtected()) {
2212          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2213    }    }
2214    if (isLazy() || right.isLazy())    MAKELAZYBINSELF(right,DIV)
2215    {    exclusiveWrite();
2216      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=    binaryOp(right,divides<double>());
2217          m_data=c->getPtr();    return (*this);
     return (*this);  
   }  
   else  
   {  
     binaryOp(right,divides<double>());  
     return (*this);  
   }  
2218  }  }
2219    
2220  Data&  Data&
# Line 2132  Data::operator/=(const boost::python::ob Line 2224  Data::operator/=(const boost::python::ob
2224          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2225    }    }
2226    Data tmp(right,getFunctionSpace(),false);    Data tmp(right,getFunctionSpace(),false);
2227    if (isLazy())    (*this)/=tmp;
2228    {    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);  
   }  
2229  }  }
2230    
2231  Data  Data
# Line 2162  Data::powO(const boost::python::object& Line 2245  Data::powO(const boost::python::object&
2245  Data  Data
2246  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2247  {  {
2248    if (isLazy() || right.isLazy())    MAKELAZYBIN(right,POW)
   {  
     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);  
     return Data(c);  
   }  
2249    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2250  }  }
2251    
# Line 2175  Data::powD(const Data& right) const Line 2254  Data::powD(const Data& right) const
2254  Data  Data
2255  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2256  {  {
2257    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,ADD)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);  
     return Data(c);  
   }  
2258    return C_TensorBinaryOperation(left, right, plus<double>());    return C_TensorBinaryOperation(left, right, plus<double>());
2259  }  }
2260    
# Line 2188  escript::operator+(const Data& left, con Line 2263  escript::operator+(const Data& left, con
2263  Data  Data
2264  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2265  {  {
2266    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,SUB)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);  
     return Data(c);  
   }  
2267    return C_TensorBinaryOperation(left, right, minus<double>());    return C_TensorBinaryOperation(left, right, minus<double>());
2268  }  }
2269    
# Line 2201  escript::operator-(const Data& left, con Line 2272  escript::operator-(const Data& left, con
2272  Data  Data
2273  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2274  {  {
2275    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,MUL)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);  
     return Data(c);  
   }  
2276    return C_TensorBinaryOperation(left, right, multiplies<double>());    return C_TensorBinaryOperation(left, right, multiplies<double>());
2277  }  }
2278    
# Line 2214  escript::operator*(const Data& left, con Line 2281  escript::operator*(const Data& left, con
2281  Data  Data
2282  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2283  {  {
2284    if (left.isLazy() || right.isLazy())    MAKELAZYBIN2(left,right,DIV)
   {  
     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);  
     return Data(c);  
   }  
2285    return C_TensorBinaryOperation(left, right, divides<double>());    return C_TensorBinaryOperation(left, right, divides<double>());
2286  }  }
2287    
# Line 2227  escript::operator/(const Data& left, con Line 2290  escript::operator/(const Data& left, con
2290  Data  Data
2291  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2292  {  {
2293    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2294    {    MAKELAZYBIN2(left,tmp,ADD)
2295      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);  
2296  }  }
2297    
2298  //  //
# Line 2240  escript::operator+(const Data& left, con Line 2300  escript::operator+(const Data& left, con
2300  Data  Data
2301  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2302  {  {
2303    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2304    {    MAKELAZYBIN2(left,tmp,SUB)
2305      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);  
2306  }  }
2307    
2308  //  //
# Line 2253  escript::operator-(const Data& left, con Line 2310  escript::operator-(const Data& left, con
2310  Data  Data
2311  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2312  {  {
2313    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2314    {    MAKELAZYBIN2(left,tmp,MUL)
2315      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);  
2316  }  }
2317    
2318  //  //
# Line 2266  escript::operator*(const Data& left, con Line 2320  escript::operator*(const Data& left, con
2320  Data  Data
2321  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2322  {  {
2323    if (left.isLazy())    Data tmp(right,left.getFunctionSpace(),false);
2324    {    MAKELAZYBIN2(left,tmp,DIV)
2325      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);  
2326  }  }
2327    
2328  //  //
# Line 2279  escript::operator/(const Data& left, con Line 2330  escript::operator/(const Data& left, con
2330  Data  Data
2331  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2332  {  {
2333    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2334    {    MAKELAZYBIN2(tmp,right,ADD)
2335      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;  
2336  }  }
2337    
2338  //  //
# Line 2292  escript::operator+(const boost::python:: Line 2340  escript::operator+(const boost::python::
2340  Data  Data
2341  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2342  {  {
2343    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2344    {    MAKELAZYBIN2(tmp,right,SUB)
2345      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;  
2346  }  }
2347    
2348  //  //
# Line 2305  escript::operator-(const boost::python:: Line 2350  escript::operator-(const boost::python::
2350  Data  Data
2351  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2352  {  {
2353    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2354    {    MAKELAZYBIN2(tmp,right,MUL)
2355      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;  
2356  }  }
2357    
2358  //  //
# Line 2318  escript::operator*(const boost::python:: Line 2360  escript::operator*(const boost::python::
2360  Data  Data
2361  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2362  {  {
2363    if (right.isLazy())    Data tmp(left,right.getFunctionSpace(),false);
2364    {    MAKELAZYBIN2(tmp,right,DIV)
2365      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;  
2366  }  }
2367    
2368    
# Line 2364  void Line 2403  void
2403  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2404                 const Data& value)                 const Data& value)
2405  {  {
 //  const DataArrayView& view=getPointDataView();  
   
2406    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2407    if (slice_region.size()!=getDataPointRank()) {    if (slice_region.size()!=getDataPointRank()) {
2408      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2409    }    }
2410      exclusiveWrite();
2411    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2412       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2413    } else {    } else {
# Line 2384  Data::setSlice(const Data& value, Line 2422  Data::setSlice(const Data& value,
2422    if (isProtected()) {    if (isProtected()) {
2423          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2424    }    }
2425    FORCERESOLVE;    forceResolve();
2426  /*  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.");  
   }*/  
2427    Data tempValue(value);    Data tempValue(value);
2428    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2429    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2431  Data::typeMatchRight(const Data& right) Line 2466  Data::typeMatchRight(const Data& right)
2466    }    }
2467  }  }
2468    
2469    // The normal TaggedValue adds the tag if it is not already present
2470    // This form does not. It throws instead.
2471    // This is because the names are maintained by the domain and cannot be added
2472    // without knowing the tag number to map it to.
2473  void  void
2474  Data::setTaggedValueByName(std::string name,  Data::setTaggedValueByName(std::string name,
2475                             const boost::python::object& value)                             const boost::python::object& value)
2476  {  {
2477       if (getFunctionSpace().getDomain()->isValidTagName(name)) {       if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2478      FORCERESOLVE;      forceResolve();
2479        exclusiveWrite();
2480          int tagKey=getFunctionSpace().getDomain()->getTag(name);          int tagKey=getFunctionSpace().getDomain()->getTag(name);
2481          setTaggedValue(tagKey,value);          setTaggedValue(tagKey,value);
2482       }       }
2483         else
2484         {                  // The
2485        throw DataException("Error - unknown tag in setTaggedValueByName.");
2486         }
2487  }  }
2488    
2489  void  void
2490  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2491                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2450  Data::setTaggedValue(int tagKey, Line 2495  Data::setTaggedValue(int tagKey,
2495    }    }
2496    //    //
2497    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2498    FORCERESOLVE;    forceResolve();
2499      exclusiveWrite();
2500    if (isConstant()) tag();    if (isConstant()) tag();
2501    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]));  
   }  
2502    
2503    DataVector temp_data2;    DataVector temp_data2;
2504    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromArray(w,1);
2505    
2506    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2507  }  }
2508    
2509    
# Line 2478  Data::setTaggedValueFromCPP(int tagKey, Line 2518  Data::setTaggedValueFromCPP(int tagKey,
2518    }    }
2519    //    //
2520    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2521    FORCERESOLVE;    forceResolve();
2522    if (isConstant()) tag();    if (isConstant()) tag();
2523      exclusiveWrite();
2524    //    //
2525    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2526    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2553  escript::C_GeneralTensorProduct(Data& ar
2553    // 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().
2554    
2555    // deal with any lazy data    // deal with any lazy data
2556    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2557    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2558      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2559      {
2560        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2561        return Data(c);
2562      }
2563    
2564    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2565    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2636  escript::C_GeneralTensorProduct(Data& ar
2636       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
2637    }    }
2638    
2639      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2640      {
2641         ostringstream os;
2642         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2643         throw DataException(os.str());
2644      }
2645    
2646    // Declare output Data object    // Declare output Data object
2647    Data res;    Data res;
2648    
2649    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2650      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2651      double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));      const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2652      double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2653      double *ptr_2 = &(res.getDataAtOffset(0));      double *ptr_2 = &(res.getDataAtOffsetRW(0));
2654      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);
2655    }    }
2656    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 2671  escript::C_GeneralTensorProduct(Data& ar
2671    
2672      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2673      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2674      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));  
2675    
2676        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2677        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2678    
2679      // Compute an MVP for the default      // Compute an MVP for the default
2680      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 2683  escript::C_GeneralTensorProduct(Data& ar
2683      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
2684      for (i=lookup_1.begin();i!=lookup_1.end();i++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2685        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);  
2686    
2687        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2688        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2689            
2690        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);
2691      }      }
# Line 2667  escript::C_GeneralTensorProduct(Data& ar Line 2709  escript::C_GeneralTensorProduct(Data& ar
2709        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2710          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2711          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2712          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2713          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2714          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2715          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);
2716        }        }
2717      }      }
# Line 2693  escript::C_GeneralTensorProduct(Data& ar Line 2735  escript::C_GeneralTensorProduct(Data& ar
2735    
2736      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2737      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2738      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2739      // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2740  //     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));  
2741    
2742      // Compute an MVP for the default      // Compute an MVP for the default
2743      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 2745  escript::C_GeneralTensorProduct(Data& ar
2745      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2746      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
2747      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);  
2748    
2749        tmp_2->addTag(i->first);        tmp_2->addTag(i->first);
2750        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2751        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2752        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);
2753      }      }
2754    
# Line 2739  escript::C_GeneralTensorProduct(Data& ar Line 2769  escript::C_GeneralTensorProduct(Data& ar
2769      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2770      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2771    
2772  //     // Get the views      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2773  //     DataArrayView view_0 = tmp_0->getDefaultValue();      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2774  //     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));  
   
2775    
2776      // Compute an MVP for the default      // Compute an MVP for the default
2777      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 2788  escript::C_GeneralTensorProduct(Data& ar
2788      // Compute an MVP for each tag      // Compute an MVP for each tag
2789      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2790      for (i=lookup_2.begin();i!=lookup_2.end();i++) {      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2791  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2792  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2793  //       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));  
2794    
2795        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);
2796      }      }
# Line 2799  escript::C_GeneralTensorProduct(Data& ar Line 2812  escript::C_GeneralTensorProduct(Data& ar
2812      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2813      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2814        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
2815        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2816        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2817          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2818          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2819          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2820          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2821          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);
2822        }        }
2823      }      }
# Line 2828  escript::C_GeneralTensorProduct(Data& ar Line 2841  escript::C_GeneralTensorProduct(Data& ar
2841        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2842          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2843          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2844          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2845          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2846          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2847          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);
2848        }        }
2849      }      }
# Line 2853  escript::C_GeneralTensorProduct(Data& ar Line 2866  escript::C_GeneralTensorProduct(Data& ar
2866      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2867      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2868        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2869        double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2870        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2871          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2872          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2873          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2874          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2875          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);
2876        }        }
2877      }      }
# Line 2883  escript::C_GeneralTensorProduct(Data& ar Line 2896  escript::C_GeneralTensorProduct(Data& ar
2896          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2897          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2898          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2899          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));          const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2900          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2901          double *ptr_2 = &(res.getDataAtOffset(offset_2));          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2902          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);
2903        }        }
2904      }      }
# Line 2923  Data::borrowReadyPtr() const Line 2936  Data::borrowReadyPtr() const
2936  std::string  std::string
2937  Data::toString() const  Data::toString() const
2938  {  {
2939      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2940      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2941        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2942      {      {
2943      stringstream temp;      stringstream temp;
2944      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 2948  Data::toString() const
2948  }  }
2949    
2950    
2951    // This method is not thread-safe
2952    DataTypes::ValueType::reference
2953    Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
2954    {
2955        checkExclusiveWrite();
2956        return getReady()->getDataAtOffsetRW(i);
2957    }
2958    
2959    // This method is not thread-safe
2960  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2961  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const  Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
2962  {  {
2963      if (isLazy())      forceResolve();
2964      {      return getReady()->getDataAtOffsetRO(i);
     throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");  
     }  
     return getReady()->getDataAtOffset(i);  
2965  }  }
2966    
2967    
2968  DataTypes::ValueType::reference  // DataTypes::ValueType::const_reference
2969  Data::getDataAtOffset(DataTypes::ValueType::size_type i)  // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
2970  {  // {
2971  //     if (isLazy())  //     if (isLazy())
2972  //     {  //     {
2973  //  throw DataException("getDataAtOffset not permitted on lazy data.");  //  throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
2974  //     }  //     }
2975      FORCERESOLVE;  //     return getReady()->getDataAtOffsetRO(i);
2976      return getReady()->getDataAtOffset(i);  // }
2977  }  
2978    
2979  DataTypes::ValueType::const_reference  DataTypes::ValueType::const_reference
2980  Data::getDataPoint(int sampleNo, int dataPointNo) const  Data::getDataPointRO(int sampleNo, int dataPointNo)
2981  {  {
2982      forceResolve();
2983    if (!isReady())    if (!isReady())
2984    {    {
2985      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.");
2986    }    }
2987    else    else
2988    {    {
2989      const DataReady* dr=getReady();      const DataReady* dr=getReady();
2990      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));      return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
2991    }    }
2992  }  }
2993    
2994    
2995    
2996    
2997  DataTypes::ValueType::reference  DataTypes::ValueType::reference
2998  Data::getDataPoint(int sampleNo, int dataPointNo)  Data::getDataPointRW(int sampleNo, int dataPointNo)
2999  {  {
3000    FORCERESOLVE;    checkExclusiveWrite();
3001    if (!isReady())    DataReady* dr=getReady();
3002    {    return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
3003      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");  }
3004    }  
3005    else  BufferGroup*
3006    {  Data::allocSampleBuffer() const
3007      DataReady* dr=getReady();  {
3008      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));       if (isLazy())
3009    }       {
3010        #ifdef _OPENMP
3011        int tnum=omp_get_max_threads();
3012        #else
3013        int tnum=1;
3014        #endif
3015        return new BufferGroup(getSampleBufferSize(),tnum);
3016         }
3017         else
3018         {
3019        return NULL;
3020         }
3021    }
3022    
3023    void
3024    Data::freeSampleBuffer(BufferGroup* bufferg)
3025    {
3026         if (bufferg!=0)
3027         {
3028        delete bufferg;
3029         }
3030    }
3031    
3032    
3033    Data
3034    Data::interpolateFromTable(boost::python::object table, double Amin, double Astep,
3035             double undef, Data& B, double Bmin, double Bstep)
3036    {
3037        WrappedArray t(table);
3038        return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep);
3039    }
3040    
3041            
3042    Data
3043    Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3044                           double undef, Data& B, double Bmin, double Bstep)
3045    {
3046        int error=0;
3047        if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3048        {
3049            throw DataException("Inputs to 2D interpolation must be scalar");
3050        }
3051        if (table.getRank()!=2)
3052        {
3053        throw DataException("Table for 2D interpolation must be 2D");
3054        }
3055        if (getFunctionSpace()!=B.getFunctionSpace())
3056        {
3057        Data n=B.interpolate(getFunctionSpace());
3058        return interpolateFromTable2D(table, Amin, Astep, undef,
3059            n , Bmin, Bstep);
3060        }
3061        if (!isExpanded())
3062        {
3063        expand();
3064        }
3065        if (!B.isExpanded())
3066        {
3067        B.expand();
3068        }
3069        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3070        do                                   // to make breaks useful
3071        {
3072        try
3073        {
3074            int numpts=getNumDataPoints();
3075            const DataVector& adat=getReady()->getVectorRO();
3076            const DataVector& bdat=B.getReady()->getVectorRO();
3077            DataVector& rdat=res.getReady()->getVectorRW();
3078            const DataTypes::ShapeType& ts=table.getShape();
3079            for (int l=0; l<numpts; ++l)
3080            {
3081            double a=adat[l];
3082            double b=bdat[l];
3083            int x=static_cast<int>((a-Amin)/Astep);
3084            int y=static_cast<int>((b-Bmin)/Bstep);
3085            if ((a<Amin) || (b<Bmin))
3086            {
3087                error=1;
3088                break;  
3089            }
3090            if ((x>=(ts[0]-1)) || (y>=(ts[1]-1)))
3091            {
3092                error=1;
3093                break;
3094            }
3095            else        // x and y are in bounds
3096            {
3097                double sw=table.getElt(x,y);
3098                double nw=table.getElt(x,y+1);
3099                double se=table.getElt(x+1,y);
3100                double ne=table.getElt(x+1,y+1);
3101                if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3102                {
3103                error=2;
3104                break;
3105                }
3106                // map x*Astep <= a << (x+1)*Astep to [-1,1]
3107                // same with b
3108                double la = 2.0*(a-(x*Astep))/Astep-1;
3109                double lb = 2.0*(b-(y*Bstep))/Bstep-1;
3110                rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3111                     (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3112            }
3113            }
3114        } catch (DataException d)
3115        {
3116            error=3;
3117            break;
3118        }
3119        } while (false);
3120    #ifdef PASO_MPI
3121        int rerror=0;
3122        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3123        error=rerror;
3124    #endif
3125        if (error)
3126        {
3127        switch (error)
3128        {
3129        case 1: throw DataException("Point out of bounds");
3130        case 2: throw DataException("Interpolated value too large");
3131        default:
3132            throw DataException("Unknown error in interpolation");      
3133        }
3134        }
3135        return res;
3136  }  }
3137    
3138    
# Line 3000  Data::print() Line 3148  Data::print()
3148    {    {
3149      printf( "[%6d]", i );      printf( "[%6d]", i );
3150      for( j=0; j<getNumDataPointsPerSample(); j++ )      for( j=0; j<getNumDataPointsPerSample(); j++ )
3151        printf( "\t%10.7g", (getSampleData(i))[j] );        printf( "\t%10.7g", (getSampleDataRW(i))[j] );    // doesn't really need RW access
3152      printf( "\n" );      printf( "\n" );
3153    }    }
3154  }  }
# Line 3020  Data::dump(const std::string fileName) c Line 3168  Data::dump(const std::string fileName) c
3168            return m_data->dump(fileName);            return m_data->dump(fileName);
3169      }      }
3170    }    }
3171    catch (exception& e)    catch (std::exception& e)
3172    {    {
3173          cout << e.what() << endl;          cout << e.what() << endl;
3174    }    }
# Line 3062  Data::get_MPIComm() const Line 3210  Data::get_MPIComm() const
3210  #endif  #endif
3211  }  }
3212    
   

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

  ViewVC Help
Powered by ViewVC 1.1.26