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