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