/[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 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC revision 2724 by jfenwick, Sun Oct 18 23:50:18 2009 UTC
# Line 1  Line 1 
 // $Id$  
1    
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Licensed under the Open Software License version 3.0
10   *                                                          *  * http://www.opensource.org/licenses/osl-3.0.php
11   ************************************************************  *
12  */  *******************************************************/
13    
14    
15  #include "Data.h"  #include "Data.h"
16    
17  #include "DataExpanded.h"  #include "DataExpanded.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
19  #include "DataTagged.h"  #include "DataTagged.h"
20  #include "DataEmpty.h"  #include "DataEmpty.h"
21  #include "DataArray.h"  #include "DataLazy.h"
 #include "DataArrayView.h"  
 #include "DataProf.h"  
22  #include "FunctionSpaceFactory.h"  #include "FunctionSpaceFactory.h"
23  #include "AbstractContinuousDomain.h"  #include "AbstractContinuousDomain.h"
24  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
25    #include "FunctionSpaceException.h"
26    #include "EscriptParams.h"
27    
28    extern "C" {
29    #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;
45  using namespace boost;  using namespace boost;
46  using namespace escript;  using namespace escript;
47    
48  #if defined DOPROF  // ensure the current object is not a DataLazy
49  //  // The idea was that we could add an optional warning whenever a resolve is forced
50  // global table of profiling data for all Data objects  // #define forceResolve() if (isLazy()) {#resolve();}
51  DataProf dataProfTable;  
52  #endif  #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    #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
91    
92    namespace
93    {
94    
95    template <class ARR>
96    inline
97    boost::python::tuple
98    pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
99    {
100        using namespace boost::python;
101        using boost::python::tuple;
102        using boost::python::list;
103    
104        list l;
105        unsigned int dim=shape[0];
106        for (size_t i=0;i<dim;++i)
107        {
108        l.append(v[i+offset]);
109        }
110        return tuple(l);
111    }
112    
113    template <class ARR>
114    inline
115    boost::python::tuple
116    pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
117    {
118        using namespace boost::python;
119        using boost::python::tuple;
120        using boost::python::list;
121    
122        unsigned int shape0=shape[0];
123        unsigned int shape1=shape[1];
124        list lj;
125        for (size_t j=0;j<shape0;++j)
126        {
127            list li;
128        for (size_t i=0;i<shape1;++i)
129        {
130            li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
131        }
132        lj.append(tuple(li));
133        }
134        return tuple(lj);
135    }
136    
137    template <class ARR>
138    inline
139    boost::python::tuple
140    pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
141    {
142        using namespace boost::python;
143        using boost::python::tuple;
144        using boost::python::list;
145    
146        unsigned int shape0=shape[0];
147        unsigned int shape1=shape[1];
148        unsigned int shape2=shape[2];
149    
150        list lk;
151        for (size_t k=0;k<shape0;++k)
152        {
153            list lj;
154        for (size_t j=0;j<shape1;++j)
155        {
156            list li;
157            for (size_t i=0;i<shape2;++i)
158            {
159                    li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
160                }
161            lj.append(tuple(li));
162            }
163            lk.append(tuple(lj));
164        }
165        return tuple(lk);
166    }
167    
168    template <class ARR>
169    inline
170    boost::python::tuple
171    pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
172    {
173        using namespace boost::python;
174        using boost::python::tuple;
175        using boost::python::list;
176    
177        unsigned int shape0=shape[0];
178        unsigned int shape1=shape[1];
179        unsigned int shape2=shape[2];
180        unsigned int shape3=shape[3];
181    
182        list ll;
183        for (size_t l=0;l<shape0;++l)
184        {
185            list lk;
186        for (size_t k=0;k<shape1;++k)
187        {
188                list lj;
189                for (size_t j=0;j<shape2;++j)
190                {
191                    list li;
192                    for (size_t i=0;i<shape3;++i)
193                    {
194                        li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
195                    }
196                    lj.append(tuple(li));
197                }
198                lk.append(tuple(lj));
199        }
200            ll.append(tuple(lk));
201        }
202        return tuple(ll);
203    }
204    
205    
206    // This should be safer once the DataC RO changes have been brought in
207    template <class ARR>
208    boost::python::tuple
209    pointToTuple( const DataTypes::ShapeType& shape,ARR v)
210    {
211       using namespace boost::python;
212       using boost::python::list;
213       int rank=shape.size();
214       if (rank==0)
215       {
216        return make_tuple(v[0]);
217       }
218       else if (rank==1)
219       {
220            return pointToTuple1(shape,v,0);
221       }
222       else if (rank==2)
223       {
224        return pointToTuple2(shape,v,0);
225       }
226       else if (rank==3)
227       {
228        return pointToTuple3(shape,v,0);
229       }
230       else if (rank==4)
231       {
232        return pointToTuple4(shape,v,0);
233       }
234       else
235         throw DataException("Unknown rank in pointToTuple.");
236    }
237    
238    }  // anonymous namespace
239    
240  Data::Data()  Data::Data()
241        : m_shared(false), m_lazy(false)
242  {  {
243    //    //
244    // Default data is type DataEmpty    // Default data is type DataEmpty
245    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
246    shared_ptr<DataAbstract> temp_data(temp);  //   m_data=temp->getPtr();
247    m_data=temp_data;    set_m_data(temp->getPtr());
248    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
249  }  }
250    
251  Data::Data(double value,  Data::Data(double value,
252             const tuple& shape,             const tuple& shape,
253             const FunctionSpace& what,             const FunctionSpace& what,
254             bool expanded)             bool expanded)
255        : m_shared(false), m_lazy(false)
256  {  {
257    DataArrayView::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
258    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
259      dataPointShape.push_back(extract<const int>(shape[i]));      dataPointShape.push_back(extract<const int>(shape[i]));
260    }    }
261    DataArray temp(dataPointShape,value);  
262    initialise(temp.getView(),what,expanded);    int len = DataTypes::noValues(dataPointShape);
263      DataVector temp_data(len,value,len);
264      initialise(temp_data, dataPointShape, what, expanded);
265    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
266  }  }
267    
268  Data::Data(double value,  Data::Data(double value,
269         const DataArrayView::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
270         const FunctionSpace& what,         const FunctionSpace& what,
271             bool expanded)             bool expanded)
272        : m_shared(false), m_lazy(false)
273  {  {
274    DataArray temp(dataPointShape,value);    int len = DataTypes::noValues(dataPointShape);
275    pair<int,int> dataShape=what.getDataShape();    DataVector temp_data(len,value,len);
276    initialise(temp.getView(),what,expanded);    initialise(temp_data, dataPointShape, what, expanded);
277    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
278  }  }
279    
280  Data::Data(const Data& inData)  Data::Data(const Data& inData)
281        : m_shared(false), m_lazy(false)
282  {  {
283    m_data=inData.m_data;    set_m_data(inData.m_data);
284    m_protected=inData.isProtected();    m_protected=inData.isProtected();
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
285  }  }
286    
287    
288  Data::Data(const Data& inData,  Data::Data(const Data& inData,
289             const DataArrayView::RegionType& region)             const DataTypes::RegionType& region)
290        : m_shared(false), m_lazy(false)
291  {  {
292      DataAbstract_ptr dat=inData.m_data;
293      if (inData.isLazy())
294      {
295        dat=inData.m_data->resolve();
296      }
297      else
298      {
299        dat=inData.m_data;
300      }
301    //    //
302    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
303    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
304    shared_ptr<DataAbstract> temp_data(tmp);    set_m_data(DataAbstract_ptr(tmp));
   m_data=temp_data;  
305    m_protected=false;    m_protected=false;
306  #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
307  }  }
308    
309  Data::Data(const Data& inData,  Data::Data(const Data& inData,
310             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
311        : m_shared(false), m_lazy(false)
312  {  {
313  #if defined DOPROF    if (inData.isEmpty())
314    // create entry in global profiling table for this object    {
315    profData = dataProfTable.newData();      throw DataException("Error - will not interpolate for instances of DataEmpty.");
316  #endif    }
317    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
318      m_data=inData.m_data;      set_m_data(inData.m_data);
319    } else {    }
320      #if defined DOPROF    else
321      profData->interpolate++;    {
322      #endif  
323      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space
324      // Note: Must use a reference or pointer to a derived object        if (!inData.probeInterpolation(functionspace))
325      // in order to get polymorphic behaviour. Shouldn't really        {           // Even though this is constant, we still need to check whether interpolation is allowed
326      // be able to create an instance of AbstractDomain but that was done      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
327      // as a boost:python work around which may no longer be required.        }
328      const AbstractDomain& inDataDomain=inData.getDomain();        // if the data is not lazy, this will just be a cast to DataReady
329      if  (inDataDomain==functionspace.getDomain()) {        DataReady_ptr dr=inData.m_data->resolve();
330        inDataDomain.interpolateOnDomain(tmp,inData);        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
331    //       m_data=DataAbstract_ptr(dc);
332          set_m_data(DataAbstract_ptr(dc));
333      } else {      } else {
334        inDataDomain.interpolateACross(tmp,inData);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
335          // Note: Must use a reference or pointer to a derived object
336          // in order to get polymorphic behaviour. Shouldn't really
337          // be able to create an instance of AbstractDomain but that was done
338          // as a boost:python work around which may no longer be required.
339          /*const AbstractDomain& inDataDomain=inData.getDomain();*/
340          const_Domain_ptr inDataDomain=inData.getDomain();
341          if  (inDataDomain==functionspace.getDomain()) {
342            inDataDomain->interpolateOnDomain(tmp,inData);
343          } else {
344            inDataDomain->interpolateACross(tmp,inData);
345          }
346    //       m_data=tmp.m_data;
347          set_m_data(tmp.m_data);
348      }      }
     m_data=tmp.m_data;  
349    }    }
350    m_protected=false;    m_protected=false;
351  }  }
352    
353  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(DataAbstract* underlyingdata)
354             const DataTagged::ValueListType & values,      : m_shared(false), m_lazy(false)
            const DataArrayView& defaultValue,  
            const FunctionSpace& what,  
            bool expanded)  
355  {  {
356    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);      set_m_data(underlyingdata->getPtr());
357    shared_ptr<DataAbstract> temp_data(temp);      m_protected=false;
   m_data=temp_data;  
   m_protected=false;  
   if (expanded) {  
     expand();  
   }  
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
358  }  }
359    
360  Data::Data(const numeric::array& value,  Data::Data(DataAbstract_ptr underlyingdata)
361         const FunctionSpace& what,      : m_shared(false), m_lazy(false)
            bool expanded)  
362  {  {
363    initialise(value,what,expanded);      set_m_data(underlyingdata);
364    m_protected=false;      m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
365  }  }
366    
367  Data::Data(const DataArrayView& value,  Data::Data(const DataTypes::ValueType& value,
368         const FunctionSpace& what,           const DataTypes::ShapeType& shape,
369             bool expanded)                   const FunctionSpace& what,
370                     bool expanded)
371        : m_shared(false), m_lazy(false)
372  {  {
373    initialise(value,what,expanded);     initialise(value,shape,what,expanded);
374    m_protected=false;     m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
375  }  }
376    
377    
378  Data::Data(const object& value,  Data::Data(const object& value,
379         const FunctionSpace& what,         const FunctionSpace& what,
380             bool expanded)             bool expanded)
381        : m_shared(false), m_lazy(false)
382  {  {
383    numeric::array asNumArray(value);    WrappedArray w(value);
384    initialise(asNumArray,what,expanded);    initialise(w,what,expanded);
385    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
386  }  }
387    
388    
389  Data::Data(const object& value,  Data::Data(const object& value,
390             const Data& other)             const Data& other)
391        : m_shared(false), m_lazy(false)
392  {  {
393    //    WrappedArray w(value);
394    // Create DataConstant using the given value and all other parameters  
395    // copied from other. If value is a rank 0 object this Data    // extract the shape of the array
396    // will assume the point data shape of other.    const DataTypes::ShapeType& tempShape=w.getShape();
397    DataArray temp(value);    if (w.getRank()==0) {
398    if (temp.getView().getRank()==0) {  
399      //  
400      // Create a DataArray with the scalar value for all elements      // get the space for the data vector
401      DataArray temp2(other.getPointDataView().getShape(),temp.getView()());      int len1 = DataTypes::noValues(tempShape);
402      initialise(temp2.getView(),other.getFunctionSpace(),false);      DataVector temp_data(len1, 0.0, len1);
403        temp_data.copyFromArray(w,1);
404    
405        int len = DataTypes::noValues(other.getDataPointShape());
406    
407        DataVector temp2_data(len, temp_data[0], len);
408        DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
409    //     m_data=DataAbstract_ptr(t);
410        set_m_data(DataAbstract_ptr(t));
411    
412    } else {    } else {
413      //      //
414      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
415      initialise(temp.getView(),other.getFunctionSpace(),false);      DataConstant* t=new DataConstant(w,other.getFunctionSpace());
416    //     m_data=DataAbstract_ptr(t);
417        set_m_data(DataAbstract_ptr(t));
418    }    }
419    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
420  }  }
421    
422  Data::~Data()  Data::~Data()
423  {  {
424      set_m_data(DataAbstract_ptr());
425    }
426    
427    
428    // only call in thread safe contexts.
429    // This method should be atomic
430    void Data::set_m_data(DataAbstract_ptr p)
431    {
432      if (m_data.get()!=0)  // release old ownership
433      {
434        m_data->removeOwner(this);
435      }
436      if (p.get()!=0)
437      {
438        m_data=p;
439        m_data->addOwner(this);
440        m_shared=m_data->isShared();
441        m_lazy=m_data->isLazy();
442      }
443    }
444    
445    void Data::initialise(const WrappedArray& value,
446                     const FunctionSpace& what,
447                     bool expanded)
448    {
449      //
450      // Construct a Data object of the appropriate type.
451      // Construct the object first as there seems to be a bug which causes
452      // undefined behaviour if an exception is thrown during construction
453      // within the shared_ptr constructor.
454      if (expanded) {
455        DataAbstract* temp=new DataExpanded(value, what);
456    //     m_data=temp->getPtr();
457        set_m_data(temp->getPtr());
458      } else {
459        DataAbstract* temp=new DataConstant(value, what);
460    //     m_data=temp->getPtr();
461        set_m_data(temp->getPtr());
462      }
463    }
464    
465    
466    void
467    Data::initialise(const DataTypes::ValueType& value,
468             const DataTypes::ShapeType& shape,
469                     const FunctionSpace& what,
470                     bool expanded)
471    {
472      //
473      // Construct a Data object of the appropriate type.
474      // Construct the object first as there seems to be a bug which causes
475      // undefined behaviour if an exception is thrown during construction
476      // within the shared_ptr constructor.
477      if (expanded) {
478        DataAbstract* temp=new DataExpanded(what, shape, value);
479    //     m_data=temp->getPtr();
480        set_m_data(temp->getPtr());
481      } else {
482        DataAbstract* temp=new DataConstant(what, shape, value);
483    //     m_data=temp->getPtr();
484        set_m_data(temp->getPtr());
485      }
486  }  }
487    
488    
489  escriptDataC  escriptDataC
490  Data::getDataC()  Data::getDataC()
491  {  {
# Line 247  Data::getDataC() const Line 502  Data::getDataC() const
502    return temp;    return temp;
503  }  }
504    
505    size_t
506    Data::getSampleBufferSize() const
507    {
508      return m_data->getSampleBufferSize();
509    }
510    
511    
512  const boost::python::tuple  const boost::python::tuple
513  Data::getShapeTuple() const  Data::getShapeTuple() const
514  {  {
515    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataTypes::ShapeType& shape=getDataPointShape();
516    switch(getDataPointRank()) {    switch(getDataPointRank()) {
517       case 0:       case 0:
518          return make_tuple();          return make_tuple();
# Line 267  Data::getShapeTuple() const Line 529  Data::getShapeTuple() const
529    }    }
530  }  }
531    
532    
533    // The different name is needed because boost has trouble with overloaded functions.
534    // It can't work out what type the function is based soley on its name.
535    // There are ways to fix this involving creating function pointer variables for each form
536    // but there doesn't seem to be a need given that the methods have the same name from the python point of view
537    Data
538    Data::copySelf()
539    {
540       DataAbstract* temp=m_data->deepCopy();
541       return Data(temp);
542    }
543    
544  void  void
545  Data::copy(const Data& other)  Data::copy(const Data& other)
546  {  {
547    //    DataAbstract* temp=other.m_data->deepCopy();
548    // Perform a deep copy    DataAbstract_ptr p=temp->getPtr();
549    //   m_data=p;
550      set_m_data(p);
551    }
552    
553    
554    Data
555    Data::delay()
556    {
557      if (!isLazy())
558    {    {
559      DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());        DataLazy* dl=new DataLazy(m_data);
560      if (temp!=0) {        return Data(dl);
561        //    }
562        // Construct a DataExpanded copy    return *this;
563        DataAbstract* newData=new DataExpanded(*temp);  }
564        shared_ptr<DataAbstract> temp_data(newData);  
565        m_data=temp_data;  void
566        return;  Data::delaySelf()
567      }  {
568      if (!isLazy())
569      {
570    //  m_data=(new DataLazy(m_data))->getPtr();
571        set_m_data((new DataLazy(m_data))->getPtr());
572    }    }
573    }
574    
575    
576    // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
577    // to zero, all the tags from all the DataTags would be in the result.
578    // However since they all have the same value (0) whether they are there or not should not matter.
579    // So I have decided that for all types this method will create a constant 0.
580    // It can be promoted up as required.
581    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
582    // but we can deal with that if it arises.
583    //
584    void
585    Data::setToZero()
586    {
587      if (isEmpty())
588    {    {
589      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
     if (temp!=0) {  
       //  
       // Construct a DataTagged copy  
       DataAbstract* newData=new DataTagged(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
590    }    }
591      if (isLazy())
592    {    {
593      DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());       DataTypes::ValueType v(getNoValues(),0);
594      if (temp!=0) {       DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
595        //       DataLazy* dl=new DataLazy(dc->getPtr());
596        // Construct a DataConstant copy       set_m_data(dl->getPtr());
       DataAbstract* newData=new DataConstant(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
597    }    }
598      else
599    {    {
600      DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());       exclusiveWrite();
601      if (temp!=0) {       m_data->setToZero();
       //  
       // Construct a DataEmpty copy  
       DataAbstract* newData=new DataEmpty();  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
602    }    }
   throw DataException("Error - Copy not implemented for this Data type.");  
603  }  }
604    
605    
606  void  void
607  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
608                     const Data& mask)                     const Data& mask)
609  {  {
610    Data mask1;    // 1. Interpolate if required so all Datas use the same FS as this
611    Data mask2;    // 2. Tag or Expand so that all Data's are the same type
612      // 3. Iterate over the data vectors copying values where mask is >0
613      if (other.isEmpty() || mask.isEmpty())
614      {
615        throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
616      }
617      Data other2(other);
618      Data mask2(mask);
619      other2.resolve();
620      mask2.resolve();
621      this->resolve();
622      FunctionSpace myFS=getFunctionSpace();
623      FunctionSpace oFS=other2.getFunctionSpace();
624      FunctionSpace mFS=mask2.getFunctionSpace();
625      if (oFS!=myFS)
626      {
627         if (other2.probeInterpolation(myFS))
628         {
629        other2=other2.interpolate(myFS);
630         }
631         else
632         {
633        throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
634         }
635      }
636      if (mFS!=myFS)
637      {
638         if (mask2.probeInterpolation(myFS))
639         {
640        mask2=mask2.interpolate(myFS);
641         }
642         else
643         {
644        throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
645         }
646      }
647                // Ensure that all args have the same type
648      if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
649      {
650        this->expand();
651        other2.expand();
652        mask2.expand();
653      }
654      else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
655      {
656        this->tag();
657        other2.tag();
658        mask2.tag();
659      }
660      else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
661      {
662      }
663      else
664      {
665        throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
666      }
667      unsigned int selfrank=getDataPointRank();
668      unsigned int otherrank=other2.getDataPointRank();
669      unsigned int maskrank=mask2.getDataPointRank();
670      if ((selfrank==0) && (otherrank>0 || maskrank>0))
671      {
672        // to get here we must be copying from a large object into a scalar
673        // I am not allowing this.
674        // If you are calling copyWithMask then you are considering keeping some existing values
675        // and so I'm going to assume that you don't want your data objects getting a new shape.
676        throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
677      }
678      exclusiveWrite();
679      // Now we iterate over the elements
680      DataVector& self=getReady()->getVectorRW();;
681      const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
682      const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
683    
684    mask1 = mask.wherePositive();    if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
685    mask2.copy(mask1);    {
686        // Not allowing this combination.
687        // it is not clear what the rank of the target should be.
688        // Should it be filled with the scalar (rank stays the same);
689        // or should the target object be reshaped to be a scalar as well.
690        throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
691      }
692      if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
693      {
694        if (mvec[0]>0)      // copy whole object if scalar is >0
695        {
696            copy(other);
697        }
698        return;
699      }
700      if (isTagged())       // so all objects involved will also be tagged
701      {
702        // note the !
703        if (!((getDataPointShape()==mask2.getDataPointShape()) &&
704            ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
705        {
706            throw DataException("copyWithMask, shape mismatch.");
707        }
708    
709    mask1 *= other;      // We need to consider the possibility that tags are missing or in the wrong order
710    mask2 *= *this;      // My guiding assumption here is: All tagged Datas are assumed to have the default value for
711    mask2 = *this - mask2;      // all tags which are not explicitly defined
712    
713        const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
714        const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
715        DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
716    
717        // first, add any tags missing from other or mask
718        const DataTagged::DataMapType& olookup=optr->getTagLookup();
719            const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
720        const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
721        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
722        for (i=olookup.begin();i!=olookup.end();i++)
723        {
724               tptr->addTag(i->first);
725            }
726            for (i=mlookup.begin();i!=mlookup.end();i++) {
727               tptr->addTag(i->first);
728            }
729        // now we know that *this has all the required tags but they aren't guaranteed to be in
730        // the same order
731    
732    *this = mask1 + mask2;      // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
733        if ((selfrank==otherrank) && (otherrank==maskrank))
734        {
735            for (i=tlookup.begin();i!=tlookup.end();i++)
736            {
737                // get the target offset
738                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
739                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
740                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
741                for (int j=0;j<getDataPointSize();++j)
742                {
743                    if (mvec[j+moff]>0)
744                    {
745                        self[j+toff]=ovec[j+ooff];
746                    }
747                }
748                }
749            // now for the default value
750            for (int j=0;j<getDataPointSize();++j)
751            {
752                if (mvec[j+mptr->getDefaultOffset()]>0)
753                {
754                    self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
755                }
756            }
757        }
758        else    // other is a scalar
759        {
760            for (i=tlookup.begin();i!=tlookup.end();i++)
761            {
762                // get the target offset
763                DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
764                    DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
765                DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
766                for (int j=0;j<getDataPointSize();++j)
767                {
768                    if (mvec[j+moff]>0)
769                    {
770                        self[j+toff]=ovec[ooff];
771                    }
772                }
773                }
774            // now for the default value
775            for (int j=0;j<getDataPointSize();++j)
776            {
777                if (mvec[j+mptr->getDefaultOffset()]>0)
778                {
779                    self[j+tptr->getDefaultOffset()]=ovec[0];
780                }
781            }
782        }
783    
784        return;         // ugly
785      }
786      // mixed scalar and non-scalar operation
787      if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
788      {
789            size_t num_points=self.size();
790        // OPENMP 3.0 allows unsigned loop vars.
791        #if defined(_OPENMP) && (_OPENMP < 200805)
792        long i;
793        #else
794        size_t i;
795        #endif  
796        size_t psize=getDataPointSize();    
797        #pragma omp parallel for private(i) schedule(static)
798        for (i=0;i<num_points;++i)
799        {
800            if (mvec[i]>0)
801            {
802                self[i]=ovec[i/psize];      // since this is expanded there is one scalar
803            }                   // dest point
804        }
805        return;         // ugly!
806      }
807      // tagged data is already taken care of so we only need to worry about shapes
808      // special cases with scalars are already dealt with so all we need to worry about is shape
809      if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
810      {
811        ostringstream oss;
812        oss <<"Error - size mismatch in arguments to copyWithMask.";
813        oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
814        oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
815        oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
816        throw DataException(oss.str());
817      }
818      size_t num_points=self.size();
819    
820      // OPENMP 3.0 allows unsigned loop vars.
821    #if defined(_OPENMP) && (_OPENMP < 200805)
822      long i;
823    #else
824      size_t i;
825    #endif
826      #pragma omp parallel for private(i) schedule(static)
827      for (i=0;i<num_points;++i)
828      {
829        if (mvec[i]>0)
830        {
831           self[i]=ovec[i];
832        }
833      }
834  }  }
835    
836  bool  bool
# Line 344  Data::isExpanded() const Line 841  Data::isExpanded() const
841  }  }
842    
843  bool  bool
844    Data::actsExpanded() const
845    {
846      return m_data->actsExpanded();
847    
848    }
849    
850    
851    bool
852  Data::isTagged() const  Data::isTagged() const
853  {  {
854    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());    DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
855    return (temp!=0);    return (temp!=0);
856  }  }
857    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
858  bool  bool
859  Data::isEmpty() const  Data::isEmpty() const
860  {  {
# Line 366  Data::isConstant() const Line 869  Data::isConstant() const
869    return (temp!=0);    return (temp!=0);
870  }  }
871    
872    bool
873    Data::isLazy() const
874    {
875      return m_lazy;    // not asking m_data because we need to be able to ask this while m_data is changing
876    }
877    
878    // at the moment this is synonymous with !isLazy() but that could change
879    bool
880    Data::isReady() const
881    {
882      return (dynamic_cast<DataReady*>(m_data.get())!=0);
883    }
884    
885    
886  void  void
887  Data::setProtection()  Data::setProtection()
888  {  {
889     m_protected=true;     m_protected=true;
890  }  }
891    
892  bool  bool
893  Data::isProtected() const  Data::isProtected() const
894  {  {
895     return m_protected;     return m_protected;
896  }  }
897    
# Line 386  Data::expand() Line 903  Data::expand()
903    if (isConstant()) {    if (isConstant()) {
904      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
905      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
906      shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
907      m_data=temp_data;      set_m_data(temp->getPtr());
908    } else if (isTagged()) {    } else if (isTagged()) {
909      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
910      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
911      shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
912      m_data=temp_data;      set_m_data(temp->getPtr());
913    } else if (isExpanded()) {    } else if (isExpanded()) {
914      //      //
915      // do nothing      // do nothing
916    } else if (isEmpty()) {    } else if (isEmpty()) {
917      throw DataException("Error - Expansion of DataEmpty not possible.");      throw DataException("Error - Expansion of DataEmpty not possible.");
918      } else if (isLazy()) {
919        resolve();
920        expand();       // resolve might not give us expanded data
921    } else {    } else {
922      throw DataException("Error - Expansion not implemented for this Data type.");      throw DataException("Error - Expansion not implemented for this Data type.");
923    }    }
# Line 409  Data::tag() Line 929  Data::tag()
929    if (isConstant()) {    if (isConstant()) {
930      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
931      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
932      shared_ptr<DataAbstract> temp_data(temp);  //     m_data=temp->getPtr();
933      m_data=temp_data;      set_m_data(temp->getPtr());
934    } else if (isTagged()) {    } else if (isTagged()) {
935      // do nothing      // do nothing
936    } else if (isExpanded()) {    } else if (isExpanded()) {
937      throw DataException("Error - Creating tag data from DataExpanded not possible.");      throw DataException("Error - Creating tag data from DataExpanded not possible.");
938    } else if (isEmpty()) {    } else if (isEmpty()) {
939      throw DataException("Error - Creating tag data from DataEmpty not possible.");      throw DataException("Error - Creating tag data from DataEmpty not possible.");
940      } else if (isLazy()) {
941         DataAbstract_ptr res=m_data->resolve();
942         if (m_data->isExpanded())
943         {
944        throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
945         }
946    //      m_data=res;
947         set_m_data(res);
948         tag();
949    } else {    } else {
950      throw DataException("Error - Tagging not implemented for this Data type.");      throw DataException("Error - Tagging not implemented for this Data type.");
951    }    }
952  }  }
953    
954  void  void
955  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::resolve()
956    {
957      if (isLazy())
958      {
959    //      m_data=m_data->resolve();
960        set_m_data(m_data->resolve());
961      }
962    }
963    
964    void
965    Data::requireWrite()
966  {  {
967    m_data->reshapeDataPoint(shape);    resolve();
968      exclusiveWrite();
969    }
970    
971    Data
972    Data::oneOver() const
973    {
974      MAKELAZYOP(RECIP)
975      return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
976  }  }
977    
978  Data  Data
979  Data::wherePositive() const  Data::wherePositive() const
980  {  {
981  #if defined DOPROF    MAKELAZYOP(GZ)
982    profData->where++;    return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
 #endif  
   return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));  
983  }  }
984    
985  Data  Data
986  Data::whereNegative() const  Data::whereNegative() const
987  {  {
988  #if defined DOPROF    MAKELAZYOP(LZ)
989    profData->where++;    return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
 #endif  
   return escript::unaryOp(*this,bind2nd(less<double>(),0.0));  
990  }  }
991    
992  Data  Data
993  Data::whereNonNegative() const  Data::whereNonNegative() const
994  {  {
995  #if defined DOPROF    MAKELAZYOP(GEZ)
996    profData->where++;    return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
 #endif  
   return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));  
997  }  }
998    
999  Data  Data
1000  Data::whereNonPositive() const  Data::whereNonPositive() const
1001  {  {
1002  #if defined DOPROF    MAKELAZYOP(LEZ)
1003    profData->where++;    return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
 #endif  
   return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));  
1004  }  }
1005    
1006  Data  Data
1007  Data::whereZero(double tol) const  Data::whereZero(double tol) const
1008  {  {
1009  #if defined DOPROF  //   Data dataAbs=abs();
1010    profData->where++;  //   return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1011  #endif     MAKELAZYOPOFF(EZ,tol)
1012    Data dataAbs=abs();     return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1013    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));  
1014  }  }
1015    
1016  Data  Data
1017  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
1018  {  {
1019  #if defined DOPROF  //   Data dataAbs=abs();
1020    profData->where++;  //   return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1021  #endif    MAKELAZYOPOFF(NEZ,tol)
1022    Data dataAbs=abs();    return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1023    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));  
1024  }  }
1025    
1026  Data  Data
1027  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
1028  {  {
 #if defined DOPROF  
   profData->interpolate++;  
 #endif  
1029    return Data(*this,functionspace);    return Data(*this,functionspace);
1030  }  }
1031    
1032  bool  bool
1033  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
1034  {  {
1035    if (getFunctionSpace()==functionspace) {    return getFunctionSpace().probeInterpolation(functionspace);
     return true;  
   } else {  
     const AbstractDomain& domain=getDomain();  
     if  (domain==functionspace.getDomain()) {  
       return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());  
     } else {  
       return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());  
     }  
   }  
1036  }  }
1037    
1038  Data  Data
1039  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
1040  {  {
1041  #if defined DOPROF    if (isEmpty())
1042    profData->grad++;    {
1043  #endif      throw DataException("Error - operation not permitted on instances of DataEmpty.");
1044      }
1045      double blocktimer_start = blocktimer_time();
1046    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
1047      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
1048    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataTypes::ShapeType grad_shape=getDataPointShape();
1049    grad_shape.push_back(functionspace.getDim());    grad_shape.push_back(functionspace.getDim());
1050    Data out(0.0,grad_shape,functionspace,true);    Data out(0.0,grad_shape,functionspace,true);
1051    getDomain().setToGradient(out,*this);    getDomain()->setToGradient(out,*this);
1052      blocktimer_increment("grad()", blocktimer_start);
1053    return out;    return out;
1054  }  }
1055    
1056  Data  Data
1057  Data::grad() const  Data::grad() const
1058  {  {
1059    return gradOn(escript::function(getDomain()));    if (isEmpty())
1060      {
1061        throw DataException("Error - operation not permitted on instances of DataEmpty.");
1062      }
1063      return gradOn(escript::function(*getDomain()));
1064  }  }
1065    
1066  int  int
1067  Data::getDataPointSize() const  Data::getDataPointSize() const
1068  {  {
1069    return getPointDataView().noValues();    return m_data->getNoValues();
1070  }  }
1071    
1072  DataArrayView::ValueType::size_type  
1073    DataTypes::ValueType::size_type
1074  Data::getLength() const  Data::getLength() const
1075  {  {
1076    return m_data->getLength();    return m_data->getLength();
1077  }  }
1078    
1079  const DataArrayView::ShapeType&  
1080  Data::getDataPointShape() const  // There is no parallelism here ... elements need to be added in the correct order.
1081    //   If we could presize the list and then fill in the elements it might work
1082    //   This would need setting elements to be threadsafe.
1083    //   Having mulitple C threads calling into one interpreter is aparently a no-no.
1084    const boost::python::object
1085    Data::toListOfTuples(bool scalarastuple)
1086    {
1087        using namespace boost::python;
1088        using boost::python::list;
1089        if (get_MPISize()>1)
1090        {
1091            throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1092        }
1093        unsigned int rank=getDataPointRank();
1094        unsigned int size=getDataPointSize();
1095    
1096        int npoints=getNumDataPoints();
1097        expand();           // This will also resolve if required
1098        const DataTypes::ValueType& vec=getReady()->getVectorRO();
1099        boost::python::list temp;
1100        temp.append(object());
1101        boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints"  trick
1102        if (rank==0)
1103        {
1104            long count;
1105            if (scalarastuple)
1106            {
1107                for (count=0;count<npoints;++count)
1108                {
1109            res[count]=make_tuple(vec[count]);
1110                }
1111            }
1112            else
1113            {
1114                for (count=0;count<npoints;++count)
1115                {
1116                    res[count]=vec[count];
1117                }
1118            }
1119        }
1120        else if (rank==1)
1121        {
1122            size_t count;
1123            size_t offset=0;
1124            for (count=0;count<npoints;++count,offset+=size)
1125            {
1126                res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1127            }
1128        }
1129        else if (rank==2)
1130        {
1131            size_t count;
1132            size_t offset=0;
1133            for (count=0;count<npoints;++count,offset+=size)
1134            {
1135            res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1136            }
1137        }
1138        else if (rank==3)
1139        {
1140            size_t count;
1141            size_t offset=0;
1142            for (count=0;count<npoints;++count,offset+=size)
1143            {
1144                res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1145            }
1146        }
1147        else if (rank==4)
1148        {
1149            size_t count;
1150            size_t offset=0;
1151            for (count=0;count<npoints;++count,offset+=size)
1152            {
1153                res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1154            }
1155        }
1156        else
1157        {
1158            throw DataException("Unknown rank in ::toListOfTuples()");
1159        }
1160        return res;
1161    }
1162    
1163    const boost::python::object
1164    Data::getValueOfDataPointAsTuple(int dataPointNo)
1165    {
1166       forceResolve();
1167       if (getNumDataPointsPerSample()>0) {
1168           int sampleNo = dataPointNo/getNumDataPointsPerSample();
1169           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1170           //
1171           // Check a valid sample number has been supplied
1172           if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1173               throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1174           }
1175    
1176           //
1177           // Check a valid data point number has been supplied
1178           if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1179               throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1180           }
1181           // TODO: global error handling
1182           DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1183           return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1184      }
1185      else
1186      {
1187        // The pre-numpy method would return an empty array of the given shape
1188        // I'm going to throw an exception because if we have zero points per sample we have problems
1189        throw DataException("Error - need at least 1 datapoint per sample.");
1190      }
1191    
1192    }
1193    
1194    
1195    void
1196    Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1197  {  {
1198    return getPointDataView().getShape();      // this will throw if the value cannot be represented
1199        setValueOfDataPointToArray(dataPointNo,py_object);
1200  }  }
1201    
1202  void  void
1203  Data::fillFromNumArray(const boost::python::numeric::array num_array)  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1204  {  {
1205    if (isProtected()) {    if (isProtected()) {
1206          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1207    }    }
1208    
1209      WrappedArray w(obj);
1210    //    //
1211    // check rank    // check rank
1212    if (num_array.getrank()<getDataPointRank())    if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1213        throw DataException("Rank of numarray does not match Data object rank");        throw DataException("Rank of array does not match Data object rank");
1214    
1215    //    //
1216    // check shape of num_array    // check shape of array
1217    for (int i=0; i<getDataPointRank(); i++) {    for (unsigned int i=0; i<getDataPointRank(); i++) {
1218      if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])      if (w.getShape()[i]!=getDataPointShape()[i])
1219         throw DataException("Shape of numarray does not match Data object rank");         throw DataException("Shape of array does not match Data object rank");
1220    }    }
1221    
1222      exclusiveWrite();
1223    
1224    //    //
1225    // make sure data is expanded:    // make sure data is expanded:
1226      //
1227    if (!isExpanded()) {    if (!isExpanded()) {
1228      expand();      expand();
1229    }    }
1230      if (getNumDataPointsPerSample()>0) {
1231    //         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1232    // and copy over         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1233    m_data->copyAll(num_array);         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1234      } else {
1235           m_data->copyToDataPoint(-1, 0,w);
1236      }
1237  }  }
1238    
1239  const  void
1240  boost::python::numeric::array  Data::setValueOfDataPoint(int dataPointNo, const double value)
 Data::convertToNumArray()  
1241  {  {
1242    //    if (isProtected()) {
1243    // determine the total number of data points          throw DataException("Error - attempt to update protected Data object.");
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
   
   //  
   // create the numeric array to be returned  
   boost::python::numeric::array numArray(0.0);  
   
   //  
   // the rank of the returned numeric array will be the rank of  
   // the data points, plus one. Where the rank of the array is n,  
   // the last n-1 dimensions will be equal to the shape of the  
   // data points, whilst the first dimension will be equal to the  
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPoints);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
1244    }    }
   
1245    //    //
1246    // resize the numeric array to the shape just calculated    // make sure data is expanded:
1247    if (arrayRank==1) {    exclusiveWrite();
1248      numArray.resize(arrayShape[0]);    if (!isExpanded()) {
1249    }      expand();
   if (arrayRank==2) {  
     numArray.resize(arrayShape[0],arrayShape[1]);  
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
   }  
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
1250    }    }
1251      if (getNumDataPointsPerSample()>0) {
1252    //         int sampleNo = dataPointNo/getNumDataPointsPerSample();
1253    // loop through each data point in turn, loading the values for that data point         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1254    // into the numeric array.         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1255    int dataPoint = 0;    } else {
1256    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         m_data->copyToDataPoint(-1, 0,value);
     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {  
       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
       if (dataPointRank==0) {  
         numArray[dataPoint]=dataPointView();  
       }  
       if (dataPointRank==1) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           numArray[dataPoint][i]=dataPointView(i);  
         }  
       }  
       if (dataPointRank==2) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             numArray[dataPoint][i][j] = dataPointView(i,j);  
           }  
         }  
       }  
       if (dataPointRank==3) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             for (int k=0; k<dataPointShape[2]; k++) {  
               numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
             }  
           }  
         }  
       }  
       if (dataPointRank==4) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             for (int k=0; k<dataPointShape[2]; k++) {  
               for (int l=0; l<dataPointShape[3]; l++) {  
                 numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
               }  
             }  
           }  
         }  
       }  
       dataPoint++;  
     }  
1257    }    }
   
   //  
   // return the loaded array  
   return numArray;  
1258  }  }
1259    
1260  const  const
1261  boost::python::numeric::array  boost::python::object
1262  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1263  {  {
1264    //    // This could be lazier than it is now
1265    // Check a valid sample number has been supplied    forceResolve();
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
1266    
1267    //    // copy datapoint into a buffer
1268    // determine the number of data points per sample    // broadcast buffer to all nodes
1269    int numDataPointsPerSample = getNumDataPointsPerSample();    // convert buffer to tuple
1270      // return tuple
1271    //  
1272    // determine the rank and shape of each data point    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1273    int dataPointRank = getDataPointRank();    size_t length=DataTypes::noValues(dataPointShape);
1274    DataArrayView::ShapeType dataPointShape = getDataPointShape();  
1275      // added for the MPI communication
1276    //    double *tmpData = new double[length];
1277    // create the numeric array to be returned  
1278    boost::python::numeric::array numArray(0.0);    // updated for the MPI case
1279      if( get_MPIRank()==procNo ){
1280    //        if (getNumDataPointsPerSample()>0) {
1281    // the rank of the returned numeric array will be the rank of      int sampleNo = dataPointNo/getNumDataPointsPerSample();
1282    // the data points, plus one. Where the rank of the array is n,      int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1283    // the last n-1 dimensions will be equal to the shape of the      //
1284    // data points, whilst the first dimension will be equal to the      // Check a valid sample number has been supplied
1285    // total number of data points. Thus the array will consist of      if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1286    // a serial vector of the data points.          throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
1287    int arrayRank = dataPointRank + 1;      }
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
1288    
1289    //      //
1290    // resize the numeric array to the shape just calculated      // Check a valid data point number has been supplied
1291    if (arrayRank==1) {      if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1292      numArray.resize(arrayShape[0]);          throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
1293    }      }
1294    if (arrayRank==2) {      // TODO: global error handling
1295      numArray.resize(arrayShape[0],arrayShape[1]);      DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
   }  
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
1296    
1297    //      memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1298    // loop through each data point in turn, loading the values for that data point       }
   // into the numeric array.  
   for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {  
     DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);  
     if (dataPointRank==0) {  
       numArray[dataPoint]=dataPointView();  
     }  
     if (dataPointRank==1) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         numArray[dataPoint][i]=dataPointView(i);  
       }  
     }  
     if (dataPointRank==2) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           numArray[dataPoint][i][j] = dataPointView(i,j);  
         }  
       }  
     }  
     if (dataPointRank==3) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
           }  
         }  
       }  
     }  
     if (dataPointRank==4) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             for (int l=0; l<dataPointShape[3]; l++) {  
               numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
             }  
           }  
         }  
       }  
     }  
1299    }    }
1300    #ifdef PASO_MPI
1301      // broadcast the data to all other processes
1302      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1303    #endif
1304    
1305      boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1306      delete [] tmpData;
1307    //    //
1308    // return the loaded array    // return the loaded array
1309    return numArray;    return t;
1310    
1311  }  }
1312    
 const  
 boost::python::numeric::array  
 Data::convertToNumArrayFromDPNo(int procNo,  
                                 int sampleNo,  
                                                                 int dataPointNo)  
1313    
1314    boost::python::object
1315    Data::integrateToTuple_const() const
1316  {  {
1317      size_t length=0;    if (isLazy())
1318      int i, j, k, l, pos;    {
1319        throw DataException("Error - cannot integrate for constant lazy data.");
   //  
   // Check a valid sample number has been supplied  
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // Check a valid data point number has been supplied  
   if (dataPointNo >= getNumDataPointsPerSample()) {  
     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");  
   }  
   
   //  
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
   
   //  
   // create the numeric array to be returned  
   boost::python::numeric::array numArray(0.0);  
   
   //  
   // the shape of the returned numeric array will be the same  
   // as that of the data point  
   int arrayRank = dataPointRank;  
   DataArrayView::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]);  
1320    }    }
1321      return integrateWorker();
1322    }
1323    
1324      // added for the MPI communication  boost::python::object
1325      length=1;  Data::integrateToTuple()
1326      for( i=0; i<arrayRank; i++ )  {
1327          length *= arrayShape[i];    if (isLazy())
1328      double *tmpData = new double[length];    {
1329        expand();   // Can't do a non-resolving version of this without changing the domain object
1330    //    }         // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1331    // load the values for the data point into the numeric array.    return integrateWorker();
   
     // updated for the MPI case  
     if( get_MPIRank()==procNo ){  
         // create a view of the data if it is stored locally  
         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
           
         // pack the data from the view into tmpData for MPI communication  
         pos=0;  
         switch( dataPointRank ){  
             case 0 :  
                 tmpData[0] = dataPointView();  
                 break;  
             case 1 :          
                 for( i=0; i<dataPointShape[0]; i++ )  
                     tmpData[i]=dataPointView(i);  
                 break;  
             case 2 :          
                 for( i=0; i<dataPointShape[0]; i++ )  
                     for( j=0; j<dataPointShape[1]; j++, pos++ )  
                         tmpData[pos]=dataPointView(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]=dataPointView(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]=dataPointView(i,j,k,l);  
                 break;  
         }  
     }  
 #ifdef PASO_MPI  
         // broadcast the data to all other processes  
         MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );  
 #endif  
   
     // unpack the data  
     switch( dataPointRank ){  
         case 0 :  
             numArray[i]=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++ )  
                     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++ )  
                         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++ )  
                             tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];  
             break;  
     }  
   
     delete [] tmpData;    
 /*  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
 */  
1332    
   //  
   // return the loaded array  
   return numArray;  
1333  }  }
1334    
1335  boost::python::numeric::array  boost::python::object
1336  Data::integrate() const  Data::integrateWorker() const
1337  {  {
1338    int index;    DataTypes::ShapeType shape = getDataPointShape();
1339    int rank = getDataPointRank();    int dataPointSize = getDataPointSize();
   DataArrayView::ShapeType shape = getDataPointShape();  
   
 #if defined DOPROF  
   profData->integrate++;  
 #endif  
1340    
1341    //    //
1342    // calculate the integral values    // calculate the integral values
1343    vector<double> integrals(getDataPointSize());    vector<double> integrals(dataPointSize);
1344    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);    vector<double> integrals_local(dataPointSize);
1345      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1346    //    if (dom==0)
1347    // create the numeric array to be returned    {            
1348    // and load the array with the integral values      throw DataException("Can not integrate over non-continuous domains.");
   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];  
           }  
         }  
       }  
     }  
1349    }    }
1350    #ifdef PASO_MPI
1351      dom->setToIntegrals(integrals_local,*this);
1352      // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1353      double *tmp = new double[dataPointSize];
1354      double *tmp_local = new double[dataPointSize];
1355      for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1356      MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1357      for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1358      tuple result=pointToTuple(shape,tmp);
1359      delete[] tmp;
1360      delete[] tmp_local;
1361    #else
1362      dom->setToIntegrals(integrals,*this);
1363    /*  double *tmp = new double[dataPointSize];
1364      for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1365      tuple result=pointToTuple(shape,integrals);
1366    //   delete tmp;
1367    #endif
1368    
1369    //  
1370    // return the loaded array    return result;
   return bp_array;  
1371  }  }
1372    
1373  Data  Data
1374  Data::sin() const  Data::sin() const
1375  {  {
1376  #if defined DOPROF    MAKELAZYOP(SIN)
1377    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);  
1378  }  }
1379    
1380  Data  Data
1381  Data::cos() const  Data::cos() const
1382  {  {
1383  #if defined DOPROF    MAKELAZYOP(COS)
1384    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);  
1385  }  }
1386    
1387  Data  Data
1388  Data::tan() const  Data::tan() const
1389  {  {
1390  #if defined DOPROF    MAKELAZYOP(TAN)
1391    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);  
1392  }  }
1393    
1394  Data  Data
1395  Data::asin() const  Data::asin() const
1396  {  {
1397  #if defined DOPROF    MAKELAZYOP(ASIN)
1398    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);  
1399  }  }
1400    
1401  Data  Data
1402  Data::acos() const  Data::acos() const
1403  {  {
1404  #if defined DOPROF    MAKELAZYOP(ACOS)
1405    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);  
1406  }  }
1407    
1408    
1409  Data  Data
1410  Data::atan() const  Data::atan() const
1411  {  {
1412  #if defined DOPROF    MAKELAZYOP(ATAN)
1413    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);  
1414  }  }
1415    
1416  Data  Data
1417  Data::sinh() const  Data::sinh() const
1418  {  {
1419  #if defined DOPROF    MAKELAZYOP(SINH)
1420    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);  
1421  }  }
1422    
1423  Data  Data
1424  Data::cosh() const  Data::cosh() const
1425  {  {
1426  #if defined DOPROF    MAKELAZYOP(COSH)
1427    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);  
1428  }  }
1429    
1430  Data  Data
1431  Data::tanh() const  Data::tanh() const
1432  {  {
1433  #if defined DOPROF    MAKELAZYOP(TANH)
1434    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1435    }
1436    
1437    
1438    Data
1439    Data::erf() const
1440    {
1441    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1442      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1443    #else
1444      MAKELAZYOP(ERF)
1445      return C_TensorUnaryOperation(*this, ::erf);
1446  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);  
1447  }  }
1448    
1449  Data  Data
1450  Data::asinh() const  Data::asinh() const
1451  {  {
1452  #if defined DOPROF    MAKELAZYOP(ASINH)
1453    profData->unary++;  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1454      return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1455    #else
1456      return C_TensorUnaryOperation(*this, ::asinh);
1457  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);  
1458  }  }
1459    
1460  Data  Data
1461  Data::acosh() const  Data::acosh() const
1462  {  {
1463  #if defined DOPROF    MAKELAZYOP(ACOSH)
1464    profData->unary++;  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1465      return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1466    #else
1467      return C_TensorUnaryOperation(*this, ::acosh);
1468  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);  
1469  }  }
1470    
1471  Data  Data
1472  Data::atanh() const  Data::atanh() const
1473  {  {
1474  #if defined DOPROF    MAKELAZYOP(ATANH)
1475    profData->unary++;  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1476      return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1477    #else
1478      return C_TensorUnaryOperation(*this, ::atanh);
1479  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);  
1480  }  }
1481    
1482  Data  Data
1483  Data::log10() const  Data::log10() const
1484  {  {
1485  #if defined DOPROF    MAKELAZYOP(LOG10)
1486    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);  
1487  }  }
1488    
1489  Data  Data
1490  Data::log() const  Data::log() const
1491  {  {
1492  #if defined DOPROF    MAKELAZYOP(LOG)
1493    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);  
1494  }  }
1495    
1496  Data  Data
1497  Data::sign() const  Data::sign() const
1498  {  {
1499  #if defined DOPROF    MAKELAZYOP(SIGN)
1500    profData->unary++;    return C_TensorUnaryOperation(*this, escript::fsign);
 #endif  
   return escript::unaryOp(*this,escript::fsign);  
1501  }  }
1502    
1503  Data  Data
1504  Data::abs() const  Data::abs() const
1505  {  {
1506  #if defined DOPROF    MAKELAZYOP(ABS)
1507    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);  
1508  }  }
1509    
1510  Data  Data
1511  Data::neg() const  Data::neg() const
1512  {  {
1513  #if defined DOPROF    MAKELAZYOP(NEG)
1514    profData->unary++;    return C_TensorUnaryOperation(*this, negate<double>());
 #endif  
   return escript::unaryOp(*this,negate<double>());  
1515  }  }
1516    
1517  Data  Data
1518  Data::pos() const  Data::pos() const
1519  {  {
1520  #if defined DOPROF      // not doing lazy check here is deliberate.
1521    profData->unary++;      // since a deep copy of lazy data should be cheap, I'll just let it happen now
 #endif  
1522    Data result;    Data result;
1523    // perform a deep copy    // perform a deep copy
1524    result.copy(*this);    result.copy(*this);
# Line 1196  Data::pos() const Line 1528  Data::pos() const
1528  Data  Data
1529  Data::exp() const  Data::exp() const
1530  {  {
1531  #if defined DOPROF    MAKELAZYOP(EXP)
1532    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
 #endif  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);  
1533  }  }
1534    
1535  Data  Data
1536  Data::sqrt() const  Data::sqrt() const
1537  {  {
1538  #if defined DOPROF    MAKELAZYOP(SQRT)
1539    profData->unary++;    return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1540    }
1541    
1542    double
1543    Data::Lsup_const() const
1544    {
1545       if (isLazy())
1546       {
1547        throw DataException("Error - cannot compute Lsup for constant lazy data.");
1548       }
1549       return LsupWorker();
1550    }
1551    
1552    double
1553    Data::Lsup()
1554    {
1555       if (isLazy())
1556       {
1557        if (!actsExpanded() || CHECK_DO_CRES)
1558        {
1559            resolve();
1560        }
1561        else
1562        {
1563    #ifdef PASO_MPI
1564            return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1565    #else
1566            return lazyAlgWorker<AbsMax>(0);
1567  #endif  #endif
1568    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);      }
1569       }
1570       return LsupWorker();
1571    }
1572    
1573    double
1574    Data::sup_const() const
1575    {
1576       if (isLazy())
1577       {
1578        throw DataException("Error - cannot compute sup for constant lazy data.");
1579       }
1580       return supWorker();
1581  }  }
1582    
1583  double  double
1584  Data::Lsup() const  Data::sup()
1585  {  {
1586    double localValue, globalValue;     if (isLazy())
1587  #if defined DOPROF     {
1588    profData->reduction1++;      if (!actsExpanded() || CHECK_DO_CRES)
1589        {
1590            resolve();
1591        }
1592        else
1593        {
1594    #ifdef PASO_MPI
1595            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1596    #else
1597            return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1);
1598  #endif  #endif
1599    //      }
1600    // set the initial absolute maximum value to zero     }
1601       return supWorker();
1602    }
1603    
1604    AbsMax abs_max_func;  double
1605    localValue = algorithm(abs_max_func,0);  Data::inf_const() const
1606    {
1607       if (isLazy())
1608       {
1609        throw DataException("Error - cannot compute inf for constant lazy data.");
1610       }
1611       return infWorker();
1612    }
1613    
1614    double
1615    Data::inf()
1616    {
1617       if (isLazy())
1618       {
1619        if (!actsExpanded() || CHECK_DO_CRES)
1620        {
1621            resolve();
1622        }
1623        else
1624        {
1625  #ifdef PASO_MPI  #ifdef PASO_MPI
1626    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );          return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
   return globalValue;  
1627  #else  #else
1628    return localValue;          return lazyAlgWorker<FMin>(numeric_limits<double>::max());
1629  #endif  #endif
1630        }
1631       }
1632       return infWorker();
1633  }  }
1634    
1635    template <class BinaryOp>
1636  double  double
1637  Data::Linf() const  #ifdef PASO_MPI
1638    Data::lazyAlgWorker(double init, MPI_Op mpiop_type)
1639    #else
1640    Data::lazyAlgWorker(double init)
1641    #endif
1642  {  {
1643    double localValue, globalValue;     if (!isLazy() || !m_data->actsExpanded())
1644  #if defined DOPROF     {
1645    profData->reduction1++;      throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1646       }
1647       DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1648       EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1649       BufferGroup* bg=allocSampleBuffer();
1650       double val=init;
1651       int i=0;
1652       const size_t numsamples=getNumSamples();
1653       const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1654       BinaryOp operation;
1655       #pragma omp parallel private(i)
1656       {
1657        double localtot=init;
1658        #pragma omp for schedule(static) private(i)
1659        for (i=0;i<numsamples;++i)
1660        {
1661            size_t roffset=0;
1662            const DataTypes::ValueType* v=dl->resolveSample(*bg, i, roffset);
1663            // Now we have the sample, run operation on all points
1664            for (size_t j=0;j<samplesize;++j)
1665            {
1666            localtot=operation(localtot,(*v)[j+roffset]);
1667            }
1668        }
1669        #pragma omp critical
1670        val=operation(val,localtot);
1671       }
1672       freeSampleBuffer(bg);
1673    #ifdef PASO_MPI
1674       double globalValue;
1675       MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1676       return globalValue;
1677    #else
1678       return val;
1679  #endif  #endif
1680    }
1681    
1682    double
1683    Data::LsupWorker() const
1684    {
1685      double localValue;
1686    //    //
1687    // set the initial absolute minimum value to max double    // set the initial absolute maximum value to zero
   AbsMin abs_min_func;  
   localValue = algorithm(abs_min_func,numeric_limits<double>::max());  
1688    
1689      AbsMax abs_max_func;
1690      localValue = algorithm(abs_max_func,0);
1691  #ifdef PASO_MPI  #ifdef PASO_MPI
1692    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );    double globalValue;
1693      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1694    return globalValue;    return globalValue;
1695  #else  #else
1696    return localValue;    return localValue;
# Line 1252  Data::Linf() const Line 1698  Data::Linf() const
1698  }  }
1699    
1700  double  double
1701  Data::sup() const  Data::supWorker() const
1702  {  {
1703    double localValue, globalValue;    double localValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1704    //    //
1705    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1706    FMax fmax_func;    FMax fmax_func;
1707    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1708  #ifdef PASO_MPI  #ifdef PASO_MPI
1709      double globalValue;
1710    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1711    return globalValue;    return globalValue;
1712  #else  #else
# Line 1271  Data::sup() const Line 1715  Data::sup() const
1715  }  }
1716    
1717  double  double
1718  Data::inf() const  Data::infWorker() const
1719  {  {
1720    double localValue, globalValue;    double localValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1721    //    //
1722    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1723    FMin fmin_func;    FMin fmin_func;
1724    localValue = algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1725  #ifdef PASO_MPI  #ifdef PASO_MPI
1726      double globalValue;
1727    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1728    return globalValue;    return globalValue;
1729  #else  #else
# Line 1294  Data::inf() const Line 1736  Data::inf() const
1736  Data  Data
1737  Data::maxval() const  Data::maxval() const
1738  {  {
1739  #if defined DOPROF     MAKELAZYOP(MAXVAL)
   profData->reduction2++;  
 #endif  
1740    //    //
1741    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1742    FMax fmax_func;    FMax fmax_func;
# Line 1306  Data::maxval() const Line 1746  Data::maxval() const
1746  Data  Data
1747  Data::minval() const  Data::minval() const
1748  {  {
1749  #if defined DOPROF    MAKELAZYOP(MINVAL)
   profData->reduction2++;  
 #endif  
1750    //    //
1751    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1752    FMin fmin_func;    FMin fmin_func;
# Line 1316  Data::minval() const Line 1754  Data::minval() const
1754  }  }
1755    
1756  Data  Data
1757  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1758  {  {
1759  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1760    profData->reduction2++;       DataTypes::ShapeType s=getDataPointShape();
1761  #endif       DataTypes::ShapeType ev_shape;
1762    Trace trace_func;       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1763    return dp_algorithm(trace_func,0);       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1764         int rank=getDataPointRank();
1765         if (rank<2) {
1766            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1767         }
1768         if (axis0<0 || axis0>rank-1) {
1769            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1770         }
1771         if (axis1<0 || axis1>rank-1) {
1772             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1773         }
1774         if (axis0 == axis1) {
1775             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1776         }
1777         MAKELAZYOP2(SWAP,axis0,axis1)
1778         if (axis0 > axis1)
1779         {
1780        axis0_tmp=axis1;
1781        axis1_tmp=axis0;
1782         }
1783         else
1784         {
1785        axis0_tmp=axis0;
1786        axis1_tmp=axis1;
1787         }
1788         for (int i=0; i<rank; i++)
1789         {
1790        if (i == axis0_tmp)
1791        {
1792            ev_shape.push_back(s[axis1_tmp]);
1793        }
1794        else if (i == axis1_tmp)
1795        {
1796            ev_shape.push_back(s[axis0_tmp]);
1797        }
1798        else
1799        {
1800            ev_shape.push_back(s[i]);
1801        }
1802         }
1803         Data ev(0.,ev_shape,getFunctionSpace());
1804         ev.typeMatchRight(*this);
1805         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1806         return ev;
1807  }  }
1808    
1809  Data  Data
1810  Data::symmetric() const  Data::symmetric() const
1811  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1812       // check input       // check input
1813       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1814       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1815          if(s[0] != s[1])          if(s[0] != s[1])
1816             throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");             throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1817       }       }
1818       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
# Line 1344  Data::symmetric() const Line 1822  Data::symmetric() const
1822       else {       else {
1823          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.");
1824       }       }
1825         MAKELAZYOP(SYM)
1826       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1827       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1828       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1353  Data::symmetric() const Line 1832  Data::symmetric() const
1832  Data  Data
1833  Data::nonsymmetric() const  Data::nonsymmetric() const
1834  {  {
1835       #if defined DOPROF       MAKELAZYOP(NSYM)
         profData->unary++;  
      #endif  
1836       // check input       // check input
1837       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1838       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1839          if(s[0] != s[1])          if(s[0] != s[1])
1840             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1841          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1842          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
1843          ev_shape.push_back(s[1]);          ev_shape.push_back(s[1]);
1844          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
# Line 1372  Data::nonsymmetric() const Line 1849  Data::nonsymmetric() const
1849       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
1850          if(!(s[0] == s[2] && s[1] == s[3]))          if(!(s[0] == s[2] && s[1] == s[3]))
1851             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1852          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1853          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
1854          ev_shape.push_back(s[1]);          ev_shape.push_back(s[1]);
1855          ev_shape.push_back(s[2]);          ev_shape.push_back(s[2]);
# Line 1388  Data::nonsymmetric() const Line 1865  Data::nonsymmetric() const
1865  }  }
1866    
1867  Data  Data
1868  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1869  {  {    
1870       #if defined DOPROF       MAKELAZYOPOFF(TRACE,axis_offset)
1871          profData->unary++;       if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1872       #endif       {
1873       DataArrayView::ShapeType s=getDataPointShape();      throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
1874         }
1875         DataTypes::ShapeType s=getDataPointShape();
1876       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1877          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1878          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1879          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1880          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1881          return ev;          return ev;
1882       }       }
1883       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
1884          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1885          if (axis_offset==0) {          if (axis_offset==0) {
1886            int s2=s[2];            int s2=s[2];
1887            ev_shape.push_back(s2);            ev_shape.push_back(s2);
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1892  Data::matrixtrace(int axis_offset) const
1892          }          }
1893          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1894          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1895          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1896          return ev;          return ev;
1897       }       }
1898       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
1899          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1900          if (axis_offset==0) {          if (axis_offset==0) {
1901            ev_shape.push_back(s[2]);            ev_shape.push_back(s[2]);
1902            ev_shape.push_back(s[3]);            ev_shape.push_back(s[3]);
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1911  Data::matrixtrace(int axis_offset) const
1911      }      }
1912          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1913          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1914      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1915          return ev;          return ev;
1916       }       }
1917       else {       else {
1918          throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");          throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1919       }       }
1920  }  }
1921    
1922  Data  Data
1923  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1924  {  {    
1925  #if defined DOPROF       MAKELAZYOPOFF(TRANS,axis_offset)
1926       profData->reduction2++;       DataTypes::ShapeType s=getDataPointShape();
1927  #endif       DataTypes::ShapeType ev_shape;
      DataArrayView::ShapeType s=getDataPointShape();  
      DataArrayView::ShapeType ev_shape;  
1928       // 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]
1929       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1930       int rank=getDataPointRank();       int rank=getDataPointRank();
# Line 1455  Data::transpose(int axis_offset) const Line 1932  Data::transpose(int axis_offset) const
1932          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);          throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1933       }       }
1934       for (int i=0; i<rank; i++) {       for (int i=0; i<rank; i++) {
1935    
1936         int index = (axis_offset+i)%rank;         int index = (axis_offset+i)%rank;
1937         ev_shape.push_back(s[index]); // Append to new shape         ev_shape.push_back(s[index]); // Append to new shape
1938       }       }
# Line 1467  Data::transpose(int axis_offset) const Line 1945  Data::transpose(int axis_offset) const
1945  Data  Data
1946  Data::eigenvalues() const  Data::eigenvalues() const
1947  {  {
1948       #if defined DOPROF       if (isLazy())
1949          profData->unary++;       {
1950       #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1951        temp.resolve();
1952        return temp.eigenvalues();
1953         }
1954       // check input       // check input
1955       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1956       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1957          throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");          throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1958       if(s[0] != s[1])       if(s[0] != s[1])
1959          throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");          throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1960       // create return       // create return
1961       DataArrayView::ShapeType ev_shape(1,s[0]);       DataTypes::ShapeType ev_shape(1,s[0]);
1962       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1963       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1964       m_data->eigenvalues(ev.m_data.get());       m_data->eigenvalues(ev.m_data.get());
# Line 1487  Data::eigenvalues() const Line 1968  Data::eigenvalues() const
1968  const boost::python::tuple  const boost::python::tuple
1969  Data::eigenvalues_and_eigenvectors(const double tol) const  Data::eigenvalues_and_eigenvectors(const double tol) const
1970  {  {
1971       #if defined DOPROF       if (isLazy())
1972          profData->unary++;       {
1973       #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1974       DataArrayView::ShapeType s=getDataPointShape();      temp.resolve();
1975       if (getDataPointRank()!=2)      return temp.eigenvalues_and_eigenvectors(tol);
1976         }
1977         DataTypes::ShapeType s=getDataPointShape();
1978         if (getDataPointRank()!=2)
1979          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1980       if(s[0] != s[1])       if(s[0] != s[1])
1981          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1982       // create return       // create return
1983       DataArrayView::ShapeType ev_shape(1,s[0]);       DataTypes::ShapeType ev_shape(1,s[0]);
1984       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1985       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1986       DataArrayView::ShapeType V_shape(2,s[0]);       DataTypes::ShapeType V_shape(2,s[0]);
1987       Data V(0.,V_shape,getFunctionSpace());       Data V(0.,V_shape,getFunctionSpace());
1988       V.typeMatchRight(*this);       V.typeMatchRight(*this);
1989       m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);       m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
# Line 1507  Data::eigenvalues_and_eigenvectors(const Line 1991  Data::eigenvalues_and_eigenvectors(const
1991  }  }
1992    
1993  const boost::python::tuple  const boost::python::tuple
1994  Data::mindp() const  Data::minGlobalDataPoint() const
1995  {  {
1996    // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an    // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1997    // abort (for unknown reasons) if there are openmp directives with it in the    // abort (for unknown reasons) if there are openmp directives with it in the
1998    // surrounding function    // surrounding function
1999    
   int SampleNo;  
2000    int DataPointNo;    int DataPointNo;
2001      int ProcNo;    int ProcNo;
2002    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
2003    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
2004  }  }
2005    
2006  void  void
2007  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
2008                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
2009  {  {
2010      if (isLazy())
2011      {
2012        Data temp(*this);   // to get around the fact that you can't resolve a const Data
2013        temp.resolve();
2014        return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2015      }
2016    int i,j;    int i,j;
2017    int lowi=0,lowj=0;    int lowi=0,lowj=0;
2018    double min=numeric_limits<double>::max();    double min=numeric_limits<double>::max();
# Line 1535  Data::calc_mindp(  int& ProcNo, Line 2023  Data::calc_mindp(  int& ProcNo,
2023    int numDPPSample=temp.getNumDataPointsPerSample();    int numDPPSample=temp.getNumDataPointsPerSample();
2024    
2025    double next,local_min;    double next,local_min;
2026    int local_lowi,local_lowj;    int local_lowi=0,local_lowj=0;    
2027    
2028    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
2029    {    {
2030      local_min=min;      local_min=min;
2031      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
2032      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
2033        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
2034          next=temp.getDataPoint(i,j)();          next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2035          if (next<local_min) {          if (next<local_min) {
2036            local_min=next;            local_min=next;
2037            local_lowi=i;            local_lowi=i;
# Line 1552  Data::calc_mindp(  int& ProcNo, Line 2040  Data::calc_mindp(  int& ProcNo,
2040        }        }
2041      }      }
2042      #pragma omp critical      #pragma omp critical
2043      if (local_min<min) {      if (local_min<min) {    // If we found a smaller value than our sentinel
2044        min=local_min;        min=local_min;
2045        lowi=local_lowi;        lowi=local_lowi;
2046        lowj=local_lowj;        lowj=local_lowj;
# Line 1560  Data::calc_mindp(  int& ProcNo, Line 2048  Data::calc_mindp(  int& ProcNo,
2048    }    }
2049    
2050  #ifdef PASO_MPI  #ifdef PASO_MPI
2051      // determine the processor on which the minimum occurs    // determine the processor on which the minimum occurs
2052      next = temp.getDataPoint(lowi,lowj)();    next = temp.getDataPointRO(lowi,lowj);
2053      int lowProc = 0;    int lowProc = 0;
2054      double *globalMins = new double[get_MPISize()+1];    double *globalMins = new double[get_MPISize()+1];
2055      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );    int error;
2056          error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
2057      if( get_MPIRank()==0 ){  
2058          next = globalMins[lowProc];    if( get_MPIRank()==0 ){
2059          for( i=1; i<get_MPISize(); i++ )      next = globalMins[lowProc];
2060              if( next>globalMins[i] ){      for( i=1; i<get_MPISize(); i++ )
2061                  lowProc = i;          if( next>globalMins[i] ){
2062                  next = globalMins[i];              lowProc = i;
2063              }              next = globalMins[i];
2064            }
2065      }
2066      MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2067      DataPointNo = lowj + lowi * numDPPSample;
2068      MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2069      delete [] globalMins;
2070      ProcNo = lowProc;
2071    #else
2072      ProcNo = 0;
2073      DataPointNo = lowj + lowi * numDPPSample;
2074    #endif
2075    }
2076    
2077    
2078    const boost::python::tuple
2079    Data::maxGlobalDataPoint() const
2080    {
2081      int DataPointNo;
2082      int ProcNo;
2083      calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2084      return make_tuple(ProcNo,DataPointNo);
2085    }
2086    
2087    void
2088    Data::calc_maxGlobalDataPoint(int& ProcNo,
2089                            int& DataPointNo) const
2090    {
2091      if (isLazy())
2092      {
2093        Data temp(*this);   // to get around the fact that you can't resolve a const Data
2094        temp.resolve();
2095        return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2096      }
2097      int i,j;
2098      int highi=0,highj=0;
2099    //-------------
2100      double max= -numeric_limits<double>::max();
2101    
2102      Data temp=maxval();
2103    
2104      int numSamples=temp.getNumSamples();
2105      int numDPPSample=temp.getNumDataPointsPerSample();
2106    
2107      double next,local_max;
2108      int local_highi=0,local_highj=0;  
2109    
2110      #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2111      {
2112        local_max=max;
2113        #pragma omp for private(i,j) schedule(static)
2114        for (i=0; i<numSamples; i++) {
2115          for (j=0; j<numDPPSample; j++) {
2116            next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2117            if (next>local_max) {
2118              local_max=next;
2119              local_highi=i;
2120              local_highj=j;
2121            }
2122          }
2123        }
2124        #pragma omp critical
2125        if (local_max>max) {    // If we found a larger value than our sentinel
2126          max=local_max;
2127          highi=local_highi;
2128          highj=local_highj;
2129        }
2130      }
2131    #ifdef PASO_MPI
2132      // determine the processor on which the maximum occurs
2133      next = temp.getDataPointRO(highi,highj);
2134      int highProc = 0;
2135      double *globalMaxs = new double[get_MPISize()+1];
2136      int error;
2137      error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2138      if( get_MPIRank()==0 ){
2139        next = globalMaxs[highProc];
2140        for( i=1; i<get_MPISize(); i++ )
2141        {
2142        if( next<globalMaxs[i] )
2143        {
2144            highProc = i;
2145            next = globalMaxs[i];
2146      }      }
2147      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );      }
2148      }
2149      MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2150      DataPointNo = highj + highi * numDPPSample;  
2151      MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2152    
2153      delete [] globalMins;    delete [] globalMaxs;
2154      ProcNo = lowProc;    ProcNo = highProc;
2155  #else  #else
2156      ProcNo = 0;    ProcNo = 0;
2157      DataPointNo = highj + highi * numDPPSample;
2158  #endif  #endif
   SampleNo = lowi;  
   DataPointNo = lowj;  
2159  }  }
2160    
2161  void  void
2162  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
2163  {  {
2164      if (isEmpty())
2165      {
2166        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2167      }
2168      if (isLazy())
2169      {
2170         Data temp(*this);  // to get around the fact that you can't resolve a const Data
2171         temp.resolve();
2172         temp.saveDX(fileName);
2173         return;
2174      }
2175    boost::python::dict args;    boost::python::dict args;
2176    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2177    getDomain().saveDX(fileName,args);    getDomain()->saveDX(fileName,args);
2178    return;    return;
2179  }  }
2180    
2181  void  void
2182  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
2183  {  {
2184      if (isEmpty())
2185      {
2186        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2187      }
2188      if (isLazy())
2189      {
2190         Data temp(*this);  // to get around the fact that you can't resolve a const Data
2191         temp.resolve();
2192         temp.saveVTK(fileName);
2193         return;
2194      }
2195    boost::python::dict args;    boost::python::dict args;
2196    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
2197    getDomain().saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args,"","");
2198    return;    return;
2199  }  }
2200    
2201    
2202    
2203  Data&  Data&
2204  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
2205  {  {
2206    if (isProtected()) {    if (isProtected()) {
2207          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2208    }    }
2209  #if defined DOPROF    MAKELAZYBINSELF(right,ADD)    // for lazy + is equivalent to +=
2210    profData->binary++;    exclusiveWrite();         // Since Lazy data does not modify its leaves we only need to worry here
 #endif  
2211    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
2212    return (*this);    return (*this);
2213  }  }
# Line 1622  Data::operator+=(const boost::python::ob Line 2218  Data::operator+=(const boost::python::ob
2218    if (isProtected()) {    if (isProtected()) {
2219          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2220    }    }
2221  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2222    profData->binary++;    (*this)+=tmp;
2223  #endif    return *this;
2224    binaryOp(right,plus<double>());  }
2225    
2226    // Hmmm, operator= makes a deep copy but the copy constructor does not?
2227    Data&
2228    Data::operator=(const Data& other)
2229    {
2230      m_protected=false;        // since any changes should be caught by exclusiveWrite();
2231    //   m_data=other.m_data;
2232      set_m_data(other.m_data);
2233    return (*this);    return (*this);
2234  }  }
2235    
# Line 1635  Data::operator-=(const Data& right) Line 2239  Data::operator-=(const Data& right)
2239    if (isProtected()) {    if (isProtected()) {
2240          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2241    }    }
2242  #if defined DOPROF    MAKELAZYBINSELF(right,SUB)
2243    profData->binary++;    exclusiveWrite();
 #endif  
2244    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
2245    return (*this);    return (*this);
2246  }  }
# Line 1648  Data::operator-=(const boost::python::ob Line 2251  Data::operator-=(const boost::python::ob
2251    if (isProtected()) {    if (isProtected()) {
2252          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2253    }    }
2254  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2255    profData->binary++;    (*this)-=tmp;
 #endif  
   binaryOp(right,minus<double>());  
2256    return (*this);    return (*this);
2257  }  }
2258    
# Line 1661  Data::operator*=(const Data& right) Line 2262  Data::operator*=(const Data& right)
2262    if (isProtected()) {    if (isProtected()) {
2263          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2264    }    }
2265  #if defined DOPROF    MAKELAZYBINSELF(right,MUL)
2266    profData->binary++;    exclusiveWrite();
 #endif  
2267    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
2268    return (*this);    return (*this);
2269  }  }
2270    
2271  Data&  Data&
2272  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
2273  {  {  
2274    if (isProtected()) {    if (isProtected()) {
2275          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2276    }    }
2277  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2278    profData->binary++;    (*this)*=tmp;
 #endif  
   binaryOp(right,multiplies<double>());  
2279    return (*this);    return (*this);
2280  }  }
2281    
# Line 1687  Data::operator/=(const Data& right) Line 2285  Data::operator/=(const Data& right)
2285    if (isProtected()) {    if (isProtected()) {
2286          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2287    }    }
2288  #if defined DOPROF    MAKELAZYBINSELF(right,DIV)
2289    profData->binary++;    exclusiveWrite();
 #endif  
2290    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
2291    return (*this);    return (*this);
2292  }  }
# Line 1700  Data::operator/=(const boost::python::ob Line 2297  Data::operator/=(const boost::python::ob
2297    if (isProtected()) {    if (isProtected()) {
2298          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2299    }    }
2300  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2301    profData->binary++;    (*this)/=tmp;
 #endif  
   binaryOp(right,divides<double>());  
2302    return (*this);    return (*this);
2303  }  }
2304    
2305  Data  Data
2306  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
2307  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
2308    Data left_d(left,*this);    Data left_d(left,*this);
2309    return left_d.powD(*this);    return left_d.powD(*this);
2310  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 2312  Data::rpowO(const boost::python::object&
2312  Data  Data
2313  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
2314  {  {
2315  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2316    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
2317  }  }
2318    
2319  Data  Data
2320  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2321  {  {
2322  #if defined DOPROF    MAKELAZYBIN(right,POW)
2323    profData->binary++;    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
2324  }  }
2325    
   
2326  //  //
2327  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
2328  Data  Data
2329  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2330  {  {
2331    Data result;    MAKELAZYBIN2(left,right,ADD)
2332    //    return C_TensorBinaryOperation(left, right, plus<double>());
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
2333  }  }
2334    
2335  //  //
# Line 1763  escript::operator+(const Data& left, con Line 2337  escript::operator+(const Data& left, con
2337  Data  Data
2338  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
2339  {  {
2340    Data result;    MAKELAZYBIN2(left,right,SUB)
2341    //    return C_TensorBinaryOperation(left, right, minus<double>());
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
2342  }  }
2343    
2344  //  //
# Line 1776  escript::operator-(const Data& left, con Line 2346  escript::operator-(const Data& left, con
2346  Data  Data
2347  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2348  {  {
2349    Data result;    MAKELAZYBIN2(left,right,MUL)
2350    //    return C_TensorBinaryOperation(left, right, multiplies<double>());
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
2351  }  }
2352    
2353  //  //
# Line 1789  escript::operator*(const Data& left, con Line 2355  escript::operator*(const Data& left, con
2355  Data  Data
2356  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2357  {  {
2358    Data result;    MAKELAZYBIN2(left,right,DIV)
2359    //    return C_TensorBinaryOperation(left, right, divides<double>());
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
2360  }  }
2361    
2362  //  //
# Line 1802  escript::operator/(const Data& left, con Line 2364  escript::operator/(const Data& left, con
2364  Data  Data
2365  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
2366  {  {
2367    //    Data tmp(right,left.getFunctionSpace(),false);
2368    // Convert to DataArray format if possible    MAKELAZYBIN2(left,tmp,ADD)
2369    DataArray temp(right);    return left+tmp;
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
2370  }  }
2371    
2372  //  //
# Line 1818  escript::operator+(const Data& left, con Line 2374  escript::operator+(const Data& left, con
2374  Data  Data
2375  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2376  {  {
2377    //    Data tmp(right,left.getFunctionSpace(),false);
2378    // Convert to DataArray format if possible    MAKELAZYBIN2(left,tmp,SUB)
2379    DataArray temp(right);    return left-tmp;
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
2380  }  }
2381    
2382  //  //
# Line 1834  escript::operator-(const Data& left, con Line 2384  escript::operator-(const Data& left, con
2384  Data  Data
2385  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2386  {  {
2387    //    Data tmp(right,left.getFunctionSpace(),false);
2388    // Convert to DataArray format if possible    MAKELAZYBIN2(left,tmp,MUL)
2389    DataArray temp(right);    return left*tmp;
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
2390  }  }
2391    
2392  //  //
# Line 1850  escript::operator*(const Data& left, con Line 2394  escript::operator*(const Data& left, con
2394  Data  Data
2395  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2396  {  {
2397    //    Data tmp(right,left.getFunctionSpace(),false);
2398    // Convert to DataArray format if possible    MAKELAZYBIN2(left,tmp,DIV)
2399    DataArray temp(right);    return left/tmp;
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
2400  }  }
2401    
2402  //  //
# Line 1866  escript::operator/(const Data& left, con Line 2404  escript::operator/(const Data& left, con
2404  Data  Data
2405  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2406  {  {
2407    //    Data tmp(left,right.getFunctionSpace(),false);
2408    // Construct the result using the given value and the other parameters    MAKELAZYBIN2(tmp,right,ADD)
2409    // from right    return tmp+right;
   Data result(left,right);  
   result+=right;  
   return result;  
2410  }  }
2411    
2412  //  //
# Line 1879  escript::operator+(const boost::python:: Line 2414  escript::operator+(const boost::python::
2414  Data  Data
2415  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2416  {  {
2417    //    Data tmp(left,right.getFunctionSpace(),false);
2418    // Construct the result using the given value and the other parameters    MAKELAZYBIN2(tmp,right,SUB)
2419    // from right    return tmp-right;
   Data result(left,right);  
   result-=right;  
   return result;  
2420  }  }
2421    
2422  //  //
# Line 1892  escript::operator-(const boost::python:: Line 2424  escript::operator-(const boost::python::
2424  Data  Data
2425  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2426  {  {
2427    //    Data tmp(left,right.getFunctionSpace(),false);
2428    // Construct the result using the given value and the other parameters    MAKELAZYBIN2(tmp,right,MUL)
2429    // from right    return tmp*right;
   Data result(left,right);  
   result*=right;  
   return result;  
2430  }  }
2431    
2432  //  //
# Line 1905  escript::operator*(const boost::python:: Line 2434  escript::operator*(const boost::python::
2434  Data  Data
2435  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2436  {  {
2437    //    Data tmp(left,right.getFunctionSpace(),false);
2438    // Construct the result using the given value and the other parameters    MAKELAZYBIN2(tmp,right,DIV)
2439    // from right    return tmp/right;
   Data result(left,right);  
   result/=right;  
   return result;  
2440  }  }
2441    
 //  
 //bool escript::operator==(const Data& left, const Data& right)  
 //{  
 //  /*  
 //  NB: this operator does very little at this point, and isn't to  
 //  be relied on. Requires further implementation.  
 //  */  
 //  
 //  bool ret;  
 //  
 //  if (left.isEmpty()) {  
 //    if(!right.isEmpty()) {  
 //      ret = false;  
 //    } else {  
 //      ret = true;  
 //    }  
 //  }  
 //  
 //  if (left.isConstant()) {  
 //    if(!right.isConstant()) {  
 //      ret = false;  
 //    } else {  
 //      ret = true;  
 //    }  
 // }  
 //  
 //  if (left.isTagged()) {  
 //   if(!right.isTagged()) {  
 //      ret = false;  
 //    } else {  
 //      ret = true;  
 //    }  
 //  }  
 //  
 //  if (left.isExpanded()) {  
 //    if(!right.isExpanded()) {  
 //      ret = false;  
 //    } else {  
 //      ret = true;  
 //    }  
 //  }  
 //  
 //  return ret;  
 //}  
2442    
2443  /* TODO */  /* TODO */
2444  /* global reduction */  /* global reduction */
2445  Data  Data
2446  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
2447  {  {
   const DataArrayView& view=getPointDataView();  
2448    
2449    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2450    
2451    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=getDataPointRank()) {
2452      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2453    }    }
2454    
# Line 1977  Data::getItem(const boost::python::objec Line 2458  Data::getItem(const boost::python::objec
2458  /* TODO */  /* TODO */
2459  /* global reduction */  /* global reduction */
2460  Data  Data
2461  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataTypes::RegionType& region) const
2462  {  {
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
2463    return Data(*this,region);    return Data(*this,region);
2464  }  }
2465    
# Line 1995  Data::setItemO(const boost::python::obje Line 2473  Data::setItemO(const boost::python::obje
2473    setItemD(key,tempData);    setItemD(key,tempData);
2474  }  }
2475    
 /* TODO */  
 /* global reduction */  
2476  void  void
2477  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2478                 const Data& value)                 const Data& value)
2479  {  {
2480    const DataArrayView& view=getPointDataView();    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2481      if (slice_region.size()!=getDataPointRank()) {
   DataArrayView::RegionType slice_region=view.getSliceRegion(key);  
   if (slice_region.size()!=view.getRank()) {  
2482      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2483    }    }
2484      exclusiveWrite();
2485    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
2486       setSlice(Data(value,getFunctionSpace()),slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
2487    } else {    } else {
# Line 2014  Data::setItemD(const boost::python::obje Line 2489  Data::setItemD(const boost::python::obje
2489    }    }
2490  }  }
2491    
 /* TODO */  
 /* global reduction */  
2492  void  void
2493  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2494                 const DataArrayView::RegionType& region)                 const DataTypes::RegionType& region)
2495  {  {
2496    if (isProtected()) {    if (isProtected()) {
2497          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2498    }    }
2499  #if defined DOPROF    forceResolve();
2500    profData->slicing++;    exclusiveWrite();     // In case someone finds a way to call this without going through setItemD
 #endif  
2501    Data tempValue(value);    Data tempValue(value);
2502    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2503    typeMatchRight(tempValue);    typeMatchRight(tempValue);
2504    m_data->setSlice(tempValue.m_data.get(),region);    getReady()->setSlice(tempValue.m_data.get(),region);
2505  }  }
2506    
2507  void  void
2508  Data::typeMatchLeft(Data& right) const  Data::typeMatchLeft(Data& right) const
2509  {  {
2510      if (right.isLazy() && !isLazy())
2511      {
2512        right.resolve();
2513      }
2514    if (isExpanded()){    if (isExpanded()){
2515      right.expand();      right.expand();
2516    } else if (isTagged()) {    } else if (isTagged()) {
# Line 2047  Data::typeMatchLeft(Data& right) const Line 2523  Data::typeMatchLeft(Data& right) const
2523  void  void
2524  Data::typeMatchRight(const Data& right)  Data::typeMatchRight(const Data& right)
2525  {  {
2526      if (isLazy() && !right.isLazy())
2527      {
2528        resolve();
2529      }
2530    if (isTagged()) {    if (isTagged()) {
2531      if (right.isExpanded()) {      if (right.isExpanded()) {
2532        expand();        expand();
# Line 2060  Data::typeMatchRight(const Data& right) Line 2540  Data::typeMatchRight(const Data& right)
2540    }    }
2541  }  }
2542    
2543  /* TODO */  // The normal TaggedValue adds the tag if it is not already present
2544  /* global reduction */  // This form does not. It throws instead.
2545    // This is because the names are maintained by the domain and cannot be added
2546    // without knowing the tag number to map it to.
2547    void
2548    Data::setTaggedValueByName(std::string name,
2549                               const boost::python::object& value)
2550    {
2551         if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2552        forceResolve();
2553        exclusiveWrite();
2554            int tagKey=getFunctionSpace().getDomain()->getTag(name);
2555            setTaggedValue(tagKey,value);
2556         }
2557         else
2558         {                  // The
2559        throw DataException("Error - unknown tag in setTaggedValueByName.");
2560         }
2561    }
2562    
2563  void  void
2564  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2565                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2071  Data::setTaggedValue(int tagKey, Line 2569  Data::setTaggedValue(int tagKey,
2569    }    }
2570    //    //
2571    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2572    tag();    forceResolve();
2573      exclusiveWrite();
2574      if (isConstant()) tag();
2575      WrappedArray w(value);
2576    
2577    if (!isTagged()) {    DataVector temp_data2;
2578      throw DataException("Error - DataTagged conversion failed!!");    temp_data2.copyFromArray(w,1);
   }  
2579    
2580    //    m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setTaggedValue  
   m_data->setTaggedValue(tagKey,valueDataArray.getView());  
2581  }  }
2582    
2583  /* TODO */  
 /* global reduction */  
2584  void  void
2585  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2586                              const DataArrayView& value)                  const DataTypes::ShapeType& pointshape,
2587                                const DataTypes::ValueType& value,
2588                    int dataOffset)
2589  {  {
2590    if (isProtected()) {    if (isProtected()) {
2591          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2592    }    }
2593    //    //
2594    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2595    tag();    forceResolve();
2596      if (isConstant()) tag();
2597    if (!isTagged()) {    exclusiveWrite();
     throw DataException("Error - DataTagged conversion failed!!");  
   }  
                                                                                                                 
2598    //    //
2599    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2600    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2601  }  }
2602    
 /* TODO */  
 /* global reduction */  
2603  int  int
2604  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2605  {  {
2606    return m_data->getTagNumber(dpno);    if (isEmpty())
2607  }    {
2608        throw DataException("Error - operation not permitted on instances of DataEmpty.");
 /* TODO */  
 /* global reduction */  
 void  
 Data::setRefValue(int ref,  
                   const boost::python::numeric::array& value)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
2609    }    }
2610    //    return getFunctionSpace().getTagFromDataPointNo(dpno);
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
2611  }  }
2612    
 /* TODO */  
 /* global reduction */  
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
2613    
2614    //  ostream& escript::operator<<(ostream& o, const Data& data)
2615    // Load DataArray with values from data-points specified by ref  {
2616    m_data->getRefValue(ref,valueDataArray);    o << data.toString();
2617      return o;
2618    //  }
   // Load values from valueDataArray into return numarray  
2619    
2620    // extract the shape of the numarray  Data
2621    int rank = value.getrank();  escript::C_GeneralTensorProduct(Data& arg_0,
2622    DataArrayView::ShapeType shape;                       Data& arg_1,
2623    for (int i=0; i < rank; i++) {                       int axis_offset,
2624      shape.push_back(extract<int>(value.getshape()[i]));                       int transpose)
2625    {
2626      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2627      // SM is the product of the last axis_offset entries in arg_0.getShape().
2628    
2629      // deal with any lazy data
2630    //   if (arg_0.isLazy()) {arg_0.resolve();}
2631    //   if (arg_1.isLazy()) {arg_1.resolve();}
2632      if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2633      {
2634        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2635        return Data(c);
2636    }    }
2637    
2638    // and load the numarray with the data from the DataArray    // Interpolate if necessary and find an appropriate function space
2639    DataArrayView valueView = valueDataArray.getView();    Data arg_0_Z, arg_1_Z;
2640      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2641    if (rank==0) {      if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2642        boost::python::numeric::array temp_numArray(valueView());        arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2643        value = temp_numArray;        arg_1_Z = Data(arg_1);
2644    }      }
2645    if (rank==1) {      else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2646      for (int i=0; i < shape[0]; i++) {        arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2647        value[i] = valueView(i);        arg_0_Z =Data(arg_0);
2648      }      }
2649    }      else {
2650    if (rank==2) {        throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         value[i][j] = valueView(i,j);  
       }  
2651      }      }
2652      } else {
2653          arg_0_Z = Data(arg_0);
2654          arg_1_Z = Data(arg_1);
2655    }    }
2656    if (rank==3) {    // Get rank and shape of inputs
2657      for (int i=0; i < shape[0]; i++) {    int rank0 = arg_0_Z.getDataPointRank();
2658        for (int j=0; j < shape[1]; j++) {    int rank1 = arg_1_Z.getDataPointRank();
2659          for (int k=0; k < shape[2]; k++) {    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2660            value[i][j][k] = valueView(i,j,k);    const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2661          }  
2662        }    // Prepare for the loops of the product and verify compatibility of shapes
2663      int start0=0, start1=0;
2664      if (transpose == 0)       {}
2665      else if (transpose == 1)  { start0 = axis_offset; }
2666      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2667      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2668    
2669    
2670      // Adjust the shapes for transpose
2671      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
2672      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
2673      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2674      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2675    
2676    #if 0
2677      // For debugging: show shape after transpose
2678      char tmp[100];
2679      std::string shapeStr;
2680      shapeStr = "(";
2681      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2682      shapeStr += ")";
2683      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2684      shapeStr = "(";
2685      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2686      shapeStr += ")";
2687      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2688    #endif
2689    
2690      // Prepare for the loops of the product
2691      int SL=1, SM=1, SR=1;
2692      for (int i=0; i<rank0-axis_offset; i++)   {
2693        SL *= tmpShape0[i];
2694      }
2695      for (int i=rank0-axis_offset; i<rank0; i++)   {
2696        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2697          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2698        }
2699        SM *= tmpShape0[i];
2700      }
2701      for (int i=axis_offset; i<rank1; i++)     {
2702        SR *= tmpShape1[i];
2703      }
2704    
2705      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2706      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
2707      {         // block to limit the scope of out_index
2708         int out_index=0;
2709         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2710         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2711      }
2712    
2713      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2714      {
2715         ostringstream os;
2716         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2717         throw DataException(os.str());
2718      }
2719    
2720      // Declare output Data object
2721      Data res;
2722    
2723      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2724        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2725        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2726        const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2727        double *ptr_2 = &(res.getDataAtOffsetRW(0));
2728        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2729      }
2730      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2731    
2732        // Prepare the DataConstant input
2733        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2734        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2735    
2736        // Borrow DataTagged input from Data object
2737        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2738        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2739    
2740        // Prepare a DataTagged output 2
2741        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2742        res.tag();
2743        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2744        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2745    
2746        // Prepare offset into DataConstant
2747        int offset_0 = tmp_0->getPointOffset(0,0);
2748        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2749    
2750        const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2751        double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2752    
2753        // Compute an MVP for the default
2754        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2755        // Compute an MVP for each tag
2756        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2757        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2758        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2759          tmp_2->addTag(i->first);
2760    
2761          const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2762          double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2763        
2764          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2765      }      }
2766    
2767    }    }
2768    if (rank==4) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2769      for (int i=0; i < shape[0]; i++) {  
2770        for (int j=0; j < shape[1]; j++) {      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2771          for (int k=0; k < shape[2]; k++) {      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2772            for (int l=0; l < shape[3]; l++) {      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2773              value[i][j][k][l] = valueView(i,j,k,l);      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2774            }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2775          }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2776        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2777        int sampleNo_1,dataPointNo_1;
2778        int numSamples_1 = arg_1_Z.getNumSamples();
2779        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2780        int offset_0 = tmp_0->getPointOffset(0,0);
2781        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2782        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2783          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2784            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2785            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2786            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2787            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2788            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2789            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2790        }        }
2791      }      }
2792    
2793    }    }
2794      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2795    
2796  }      // Borrow DataTagged input from Data object
2797        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2798        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2799    
2800  void      // Prepare the DataConstant input
2801  Data::archiveData(const std::string fileName)      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2802  {      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
   cout << "Archiving Data object to: " << fileName << endl;  
2803    
2804    //      // Prepare a DataTagged output 2
2805    // Determine type of this Data object      res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2806    int dataType = -1;      res.tag();
2807        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2808        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2809    
2810    if (isEmpty()) {      // Prepare offset into DataConstant
2811      dataType = 0;      int offset_1 = tmp_1->getPointOffset(0,0);
2812      cout << "\tdataType: DataEmpty" << endl;      const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2813    }      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2814    if (isConstant()) {      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
     dataType = 1;  
     cout << "\tdataType: DataConstant" << endl;  
   }  
   if (isTagged()) {  
     dataType = 2;  
     cout << "\tdataType: DataTagged" << endl;  
   }  
   if (isExpanded()) {  
     dataType = 3;  
     cout << "\tdataType: DataExpanded" << endl;  
   }  
2815    
2816    if (dataType == -1) {      // Compute an MVP for the default
2817      throw DataException("archiveData Error: undefined dataType");      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2818    }      // Compute an MVP for each tag
2819        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2820        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2821        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2822    
2823    //        tmp_2->addTag(i->first);
2824    // Collect data items common to all Data types        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2825    int noSamples = getNumSamples();        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2826    int noDPPSample = getNumDataPointsPerSample();        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
   int functionSpaceType = getFunctionSpace().getTypeCode();  
   int dataPointRank = getDataPointRank();  
   int dataPointSize = getDataPointSize();  
   int dataLength = getLength();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
   vector<int> referenceNumbers(noSamples);  
   for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
     referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);  
   }  
   vector<int> tagNumbers(noSamples);  
   if (isTagged()) {  
     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
       tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);  
2827      }      }
2828    
2829    }    }
2830      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2831    
2832    cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;      // Borrow DataTagged input from Data object
2833    cout << "\tfunctionSpaceType: " << functionSpaceType << endl;      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2834    cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2835    
2836    //      // Borrow DataTagged input from Data object
2837    // Flatten Shape to an array of integers suitable for writing to file      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2838    int flatShape[4] = {0,0,0,0};      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
   cout << "\tshape: < ";  
   for (int dim=0; dim<dataPointRank; dim++) {  
     flatShape[dim] = dataPointShape[dim];  
     cout << dataPointShape[dim] << " ";  
   }  
   cout << ">" << endl;  
2839    
2840    //      // Prepare a DataTagged output 2
2841    // Open archive file      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2842    ofstream archiveFile;      res.tag();  // DataTagged output
2843    archiveFile.open(fileName.data(), ios::out);      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2844        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2845    
2846    if (!archiveFile.good()) {      const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2847      throw DataException("archiveData Error: problem opening archive file");      const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2848    }      double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2849    
2850    //      // Compute an MVP for the default
2851    // Write common data items to archive file      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2852    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));      // Merge the tags
2853    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2854    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2855    archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2856    archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2857    archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));        tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2858    archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));      }
2859    for (int dim = 0; dim < 4; dim++) {      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2860      archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));        tmp_2->addTag(i->first);
2861    }      }
2862    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      // Compute an MVP for each tag
2863      archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2864    }      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2865    if (isTagged()) {        const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2866      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {        const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2867        archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2868    
2869          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2870      }      }
   }  
2871    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem writing to archive file");  
2872    }    }
2873      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2874    
2875    //      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2876    // Archive underlying data values for each Data type      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2877    int noValues;      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2878    switch (dataType) {      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2879      case 0:      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2880        // DataEmpty      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2881        noValues = 0;      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2882        archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2883        cout << "\tnoValues: " << noValues << endl;      int sampleNo_0,dataPointNo_0;
2884        break;      int numSamples_0 = arg_0_Z.getNumSamples();
2885      case 1:      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2886        // DataConstant      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2887        noValues = m_data->getLength();      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2888        archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2889        cout << "\tnoValues: " << noValues << endl;        const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2890        if (m_data->archiveData(archiveFile,noValues)) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2891          throw DataException("archiveData Error: problem writing data to archive file");          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2892        }          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2893        break;          const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2894      case 2:          double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2895        // DataTagged          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
       noValues = m_data->getLength();  
       archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));  
       cout << "\tnoValues: " << noValues << endl;  
       if (m_data->archiveData(archiveFile,noValues)) {  
         throw DataException("archiveData Error: problem writing data to archive file");  
2896        }        }
2897        break;      }
2898      case 3:  
2899        // DataExpanded    }
2900        noValues = m_data->getLength();    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2901        archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));  
2902        cout << "\tnoValues: " << noValues << endl;      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2903        if (m_data->archiveData(archiveFile,noValues)) {      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2904          throw DataException("archiveData Error: problem writing data to archive file");      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2905        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2906        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2907        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2908        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2909        int sampleNo_0,dataPointNo_0;
2910        int numSamples_0 = arg_0_Z.getNumSamples();
2911        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2912        int offset_1 = tmp_1->getPointOffset(0,0);
2913        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2914        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2915          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2916            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2917            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2918            const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2919            const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2920            double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2921            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2922        }        }
2923        break;      }
2924    
2925    
2926    }    }
2927      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2928    
2929        // After finding a common function space above the two inputs have the same numSamples and num DPPS
<