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