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