/[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 2048 by phornby, Mon Nov 17 08:46:00 2008 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    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
25  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
26    #endif
27    #include "FunctionSpaceException.h"
28    #include "EscriptParams.h"
29    
30    extern "C" {
31    #include "escript/blocktimer.h"
32    }
33    
34  #include <fstream>  #include <fstream>
35  #include <algorithm>  #include <algorithm>
# Line 38  using namespace boost::python; Line 45  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();}
 DataProf dataProfTable;  
 #endif  
51    
52  Data::Data()  Data::Data()
53  {  {
54    //    //
55    // Default data is type DataEmpty    // Default data is type DataEmpty
56    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
57    shared_ptr<DataAbstract> temp_data(temp);    m_data=temp->getPtr();
   m_data=temp_data;  
58    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
59  }  }
60    
61  Data::Data(double value,  Data::Data(double value,
# Line 63  Data::Data(double value, Line 63  Data::Data(double value,
63             const FunctionSpace& what,             const FunctionSpace& what,
64             bool expanded)             bool expanded)
65  {  {
66    DataArrayView::ShapeType dataPointShape;    DataTypes::ShapeType dataPointShape;
67    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
68      dataPointShape.push_back(extract<const int>(shape[i]));      dataPointShape.push_back(extract<const int>(shape[i]));
69    }    }
70    DataArray temp(dataPointShape,value);  
71    initialise(temp.getView(),what,expanded);    int len = DataTypes::noValues(dataPointShape);
72      DataVector temp_data(len,value,len);
73      initialise(temp_data, dataPointShape, what, expanded);
74    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
75  }  }
76    
77  Data::Data(double value,  Data::Data(double value,
78         const DataArrayView::ShapeType& dataPointShape,         const DataTypes::ShapeType& dataPointShape,
79         const FunctionSpace& what,         const FunctionSpace& what,
80             bool expanded)             bool expanded)
81  {  {
82    DataArray temp(dataPointShape,value);    int len = DataTypes::noValues(dataPointShape);
83    pair<int,int> dataShape=what.getDataShape();  
84    initialise(temp.getView(),what,expanded);    DataVector temp_data(len,value,len);
85    //   DataArrayView temp_dataView(temp_data, dataPointShape);
86    
87    //   initialise(temp_dataView, what, expanded);
88      initialise(temp_data, dataPointShape, what, expanded);
89    
90    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
91  }  }
92    
93  Data::Data(const Data& inData)  Data::Data(const Data& inData)
94  {  {
95    m_data=inData.m_data;    m_data=inData.m_data;
96    m_protected=inData.isProtected();    m_protected=inData.isProtected();
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
97  }  }
98    
99    
100  Data::Data(const Data& inData,  Data::Data(const Data& inData,
101             const DataArrayView::RegionType& region)             const DataTypes::RegionType& region)
102  {  {
103      DataAbstract_ptr dat=inData.m_data;
104      if (inData.isLazy())
105      {
106        dat=inData.m_data->resolve();
107      }
108      else
109      {
110        dat=inData.m_data;
111      }
112    //    //
113    // Create Data which is a slice of another Data    // Create Data which is a slice of another Data
114    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = dat->getSlice(region);
115    shared_ptr<DataAbstract> temp_data(tmp);    m_data=DataAbstract_ptr(tmp);
   m_data=temp_data;  
116    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
117  }  }
118    
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122  #if defined DOPROF    if (inData.isEmpty())
123    // create entry in global profiling table for this object    {
124    profData = dataProfTable.newData();      throw DataException("Error - will not interpolate for instances of DataEmpty.");
125  #endif    }
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    }
129      #if defined DOPROF    else
130      profData->interpolate++;    {
131      #endif  
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space
133      // Note: Must use a reference or pointer to a derived object        if (!inData.probeInterpolation(functionspace))
134      // in order to get polymorphic behaviour. Shouldn't really        {           // Even though this is constant, we still need to check whether interpolation is allowed
135      // be able to create an instance of AbstractDomain but that was done      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");
136      // as a boost:python work around which may no longer be required.        }
137      const AbstractDomain& inDataDomain=inData.getDomain();        // if the data is not lazy, this will just be a cast to DataReady
138      if  (inDataDomain==functionspace.getDomain()) {        DataReady_ptr dr=inData.m_data->resolve();
139        inDataDomain.interpolateOnDomain(tmp,inData);        DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVector());  
140          m_data=DataAbstract_ptr(dc);
141      } else {      } else {
142        inDataDomain.interpolateACross(tmp,inData);        Data tmp(0,inData.getDataPointShape(),functionspace,true);
143          // Note: Must use a reference or pointer to a derived object
144          // in order to get polymorphic behaviour. Shouldn't really
145          // be able to create an instance of AbstractDomain but that was done
146          // as a boost:python work around which may no longer be required.
147          /*const AbstractDomain& inDataDomain=inData.getDomain();*/
148          const_Domain_ptr inDataDomain=inData.getDomain();
149          if  (inDataDomain==functionspace.getDomain()) {
150            inDataDomain->interpolateOnDomain(tmp,inData);
151          } else {
152            inDataDomain->interpolateACross(tmp,inData);
153          }
154          m_data=tmp.m_data;
155      }      }
     m_data=tmp.m_data;  
156    }    }
157    m_protected=false;    m_protected=false;
158  }  }
159    
160  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(DataAbstract* underlyingdata)
            const DataTagged::ValueListType & values,  
            const DataArrayView& defaultValue,  
            const FunctionSpace& what,  
            bool expanded)  
161  {  {
162    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);  //  m_data=shared_ptr<DataAbstract>(underlyingdata);
163    shared_ptr<DataAbstract> temp_data(temp);      m_data=underlyingdata->getPtr();
164    m_data=temp_data;      m_protected=false;
   m_protected=false;  
   if (expanded) {  
     expand();  
   }  
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
165  }  }
166    
167    Data::Data(DataAbstract_ptr underlyingdata)
168    {
169        m_data=underlyingdata;
170        m_protected=false;
171    }
172    
173    
174  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
175         const FunctionSpace& what,         const FunctionSpace& what,
176             bool expanded)             bool expanded)
177  {  {
178    initialise(value,what,expanded);    initialise(value,what,expanded);
179    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
180  }  }
181    /*
182  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
183         const FunctionSpace& what,         const FunctionSpace& what,
184             bool expanded)             bool expanded)
185  {  {
186    initialise(value,what,expanded);    initialise(value,what,expanded);
187    m_protected=false;    m_protected=false;
188  #if defined DOPROF  }*/
189    // create entry in global profiling table for this object  
190    profData = dataProfTable.newData();  Data::Data(const DataTypes::ValueType& value,
191  #endif           const DataTypes::ShapeType& shape,
192                     const FunctionSpace& what,
193                     bool expanded)
194    {
195       initialise(value,shape,what,expanded);
196       m_protected=false;
197  }  }
198    
199    
200  Data::Data(const object& value,  Data::Data(const object& value,
201         const FunctionSpace& what,         const FunctionSpace& what,
202             bool expanded)             bool expanded)
# Line 195  Data::Data(const object& value, Line 204  Data::Data(const object& value,
204    numeric::array asNumArray(value);    numeric::array asNumArray(value);
205    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
206    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
207  }  }
208    
209    
210  Data::Data(const object& value,  Data::Data(const object& value,
211             const Data& other)             const Data& other)
212  {  {
213      numeric::array asNumArray(value);
214    
215      // extract the shape of the numarray
216      DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);
217    // /*  for (int i=0; i < asNumArray.getrank(); i++) {
218    //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
219    //   }*/
220    //   // get the space for the data vector
221    //   int len = DataTypes::noValues(tempShape);
222    //   DataVector temp_data(len, 0.0, len);
223    // /*  DataArrayView temp_dataView(temp_data, tempShape);
224    //   temp_dataView.copy(asNumArray);*/
225    //   temp_data.copyFromNumArray(asNumArray);
226    
227    //    //
228    // Create DataConstant using the given value and all other parameters    // Create DataConstant using the given value and all other parameters
229    // copied from other. If value is a rank 0 object this Data    // copied from other. If value is a rank 0 object this Data
230    // will assume the point data shape of other.    // will assume the point data shape of other.
231    DataArray temp(value);  
232    if (temp.getView().getRank()==0) {    if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {
233      //  
234      // Create a DataArray with the scalar value for all elements  
235      DataArray temp2(other.getPointDataView().getShape(),temp.getView()());      // get the space for the data vector
236      initialise(temp2.getView(),other.getFunctionSpace(),false);      int len1 = DataTypes::noValues(tempShape);
237        DataVector temp_data(len1, 0.0, len1);
238        temp_data.copyFromNumArray(asNumArray);
239    
240        int len = DataTypes::noValues(other.getDataPointShape());
241    
242        DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);
243        //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());
244    //     initialise(temp2_dataView, other.getFunctionSpace(), false);
245    
246        DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
247    //     boost::shared_ptr<DataAbstract> sp(t);
248    //     m_data=sp;
249        m_data=DataAbstract_ptr(t);
250    
251    } else {    } else {
252      //      //
253      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
254      initialise(temp.getView(),other.getFunctionSpace(),false);  //     initialise(temp_dataView, other.getFunctionSpace(), false);
255        DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());
256    //     boost::shared_ptr<DataAbstract> sp(t);
257    //     m_data=sp;
258        m_data=DataAbstract_ptr(t);
259    }    }
260    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
261  }  }
262    
263  Data::~Data()  Data::~Data()
# Line 231  Data::~Data() Line 265  Data::~Data()
265    
266  }  }
267    
268    
269    
270    void
271    Data::initialise(const boost::python::numeric::array& value,
272                     const FunctionSpace& what,
273                     bool expanded)
274    {
275      //
276      // Construct a Data object of the appropriate type.
277      // Construct the object first as there seems to be a bug which causes
278      // undefined behaviour if an exception is thrown during construction
279      // within the shared_ptr constructor.
280      if (expanded) {
281        DataAbstract* temp=new DataExpanded(value, what);
282    //     boost::shared_ptr<DataAbstract> temp_data(temp);
283    //     m_data=temp_data;
284        m_data=temp->getPtr();
285      } else {
286        DataAbstract* temp=new DataConstant(value, what);
287    //     boost::shared_ptr<DataAbstract> temp_data(temp);
288    //     m_data=temp_data;
289        m_data=temp->getPtr();
290      }
291    }
292    
293    
294    void
295    Data::initialise(const DataTypes::ValueType& value,
296             const DataTypes::ShapeType& shape,
297                     const FunctionSpace& what,
298                     bool expanded)
299    {
300      //
301      // Construct a Data object of the appropriate type.
302      // Construct the object first as there seems to be a bug which causes
303      // undefined behaviour if an exception is thrown during construction
304      // within the shared_ptr constructor.
305      if (expanded) {
306        DataAbstract* temp=new DataExpanded(what, shape, value);
307    //     boost::shared_ptr<DataAbstract> temp_data(temp);
308    //     m_data=temp_data;
309        m_data=temp->getPtr();
310      } else {
311        DataAbstract* temp=new DataConstant(what, shape, value);
312    //     boost::shared_ptr<DataAbstract> temp_data(temp);
313    //     m_data=temp_data;
314        m_data=temp->getPtr();
315      }
316    }
317    
318    
319    // void
320    // Data::CompareDebug(const Data& rd)
321    // {
322    //  using namespace std;
323    //  bool mismatch=false;
324    //  std::cout << "Comparing left and right" << endl;
325    //  const DataTagged* left=dynamic_cast<DataTagged*>(m_data.get());
326    //  const DataTagged* right=dynamic_cast<DataTagged*>(rd.m_data.get());
327    //  
328    //  if (left==0)
329    //  {
330    //      cout << "left arg is not a DataTagged\n";
331    //      return;
332    //  }
333    //  
334    //  if (right==0)
335    //  {
336    //      cout << "right arg is not a DataTagged\n";
337    //      return;
338    //  }
339    //  cout << "Num elements=" << left->getVector().size() << ":" << right->getVector().size() << std::endl;
340    //  cout << "Shapes ";
341    //  if (left->getShape()==right->getShape())
342    //  {
343    //      cout << "ok\n";
344    //  }
345    //  else
346    //  {
347    //      cout << "Problem: shapes do not match\n";
348    //      mismatch=true;
349    //  }
350    //  int lim=left->getVector().size();
351    //  if (right->getVector().size()) lim=right->getVector().size();
352    //  for (int i=0;i<lim;++i)
353    //  {
354    //      if (left->getVector()[i]!=right->getVector()[i])
355    //      {
356    //          cout << "[" << i << "] value mismatch " << left->getVector()[i] << ":" << right->getVector()[i] << endl;
357    //          mismatch=true;
358    //      }
359    //  }
360    //
361    //  // still need to check the tag map
362    //  // also need to watch what is happening to function spaces, are they copied or what?
363    //
364    //  const DataTagged::DataMapType& mapleft=left->getTagLookup();
365    //  const DataTagged::DataMapType& mapright=right->getTagLookup();
366    //
367    //  if (mapleft.size()!=mapright.size())
368    //  {
369    //      cout << "Maps are different sizes " << mapleft.size() << ":" << mapright.size() << endl;
370    //      mismatch=true;
371    //      cout << "Left map\n";
372    //      DataTagged::DataMapType::const_iterator i,j;
373    //      for (i=mapleft.begin();i!=mapleft.end();++i) {
374    //          cout << "(" << i->first << "=>" << i->second << ")\n";
375    //      }
376    //      cout << "Right map\n";
377    //      for (i=mapright.begin();i!=mapright.end();++i) {
378    //          cout << "(" << i->first << "=>" << i->second << ")\n";
379    //      }
380    //      cout << "End map\n";
381    //
382    //  }
383    //
384    //  DataTagged::DataMapType::const_iterator i,j;
385    //  for (i=mapleft.begin(),j=mapright.begin();i!=mapleft.end() && j!=mapright.end();++i,++j) {
386    //     if ((i->first!=j->first) || (i->second!=j->second))
387    //     {
388    //      cout << "(" << i->first << "=>" << i->second << ")";
389    //      cout << ":(" << j->first << "=>" << j->second << ") ";
390    //      mismatch=true;
391    //            }
392    //  }
393    //  if (mismatch)
394    //  {
395    //      cout << "#Mismatch\n";
396    //  }
397    // }
398    
399  escriptDataC  escriptDataC
400  Data::getDataC()  Data::getDataC()
401  {  {
# Line 250  Data::getDataC() const Line 415  Data::getDataC() const
415  const boost::python::tuple  const boost::python::tuple
416  Data::getShapeTuple() const  Data::getShapeTuple() const
417  {  {
418    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataTypes::ShapeType& shape=getDataPointShape();
419    switch(getDataPointRank()) {    switch(getDataPointRank()) {
420       case 0:       case 0:
421          return make_tuple();          return make_tuple();
# Line 267  Data::getShapeTuple() const Line 432  Data::getShapeTuple() const
432    }    }
433  }  }
434    
435    
436    // The different name is needed because boost has trouble with overloaded functions.
437    // It can't work out what type the function is based soley on its name.
438    // There are ways to fix this involving creating function pointer variables for each form
439    // but there doesn't seem to be a need given that the methods have the same name from the python point of view
440    Data*
441    Data::copySelf()
442    {
443       DataAbstract* temp=m_data->deepCopy();
444       return new Data(temp);
445    }
446    
447  void  void
448  Data::copy(const Data& other)  Data::copy(const Data& other)
449  {  {
450    //    DataAbstract* temp=other.m_data->deepCopy();
451    // Perform a deep copy    DataAbstract_ptr p=temp->getPtr();
452    {    m_data=p;
453      DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());  }
454      if (temp!=0) {  
455        //  
456        // Construct a DataExpanded copy  Data
457        DataAbstract* newData=new DataExpanded(*temp);  Data::delay()
458        shared_ptr<DataAbstract> temp_data(newData);  {
459        m_data=temp_data;    DataLazy* dl=new DataLazy(m_data);
460        return;    return Data(dl);
461      }  }
462    }  
463    {  void
464      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());  Data::delaySelf()
465      if (temp!=0) {  {
466        //    if (!isLazy())
       // Construct a DataTagged copy  
       DataAbstract* newData=new DataTagged(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
   }  
467    {    {
468      DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());      m_data=(new DataLazy(m_data))->getPtr();
     if (temp!=0) {  
       //  
       // Construct a DataConstant copy  
       DataAbstract* newData=new DataConstant(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
469    }    }
470    }
471    
472    void
473    Data::setToZero()
474    {
475      if (isEmpty())
476    {    {
477      DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());       throw DataException("Error - Operations not permitted on instances of DataEmpty.");
     if (temp!=0) {  
       //  
       // Construct a DataEmpty copy  
       DataAbstract* newData=new DataEmpty();  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
478    }    }
479    throw DataException("Error - Copy not implemented for this Data type.");    m_data->setToZero();
480  }  }
481    
482  void  void
483  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
484                     const Data& mask)                     const Data& mask)
485  {  {
486    Data mask1;    // 1. Interpolate if required so all Datas use the same FS as this
487    Data mask2;    // 2. Tag or Expand so that all Data's are the same type
488      // 3. Iterate over the data vectors copying values where mask is >0
489      if (other.isEmpty() || mask.isEmpty())
490      {
491        throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
492      }
493      Data other2(other);
494      Data mask2(mask);
495      other2.resolve();
496      mask2.resolve();
497      this->resolve();
498      FunctionSpace myFS=getFunctionSpace();
499      FunctionSpace oFS=other2.getFunctionSpace();
500      FunctionSpace mFS=mask2.getFunctionSpace();
501      if (oFS!=myFS)
502      {
503         if (other2.probeInterpolation(myFS))
504         {
505        other2=other2.interpolate(myFS);
506         }
507         else
508         {
509        throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
510         }
511      }
512      if (mFS!=myFS)
513      {
514         if (mask2.probeInterpolation(myFS))
515         {
516        mask2=mask2.interpolate(myFS);
517         }
518         else
519         {
520        throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
521         }
522      }
523                // Ensure that all args have the same type
524      if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
525      {
526        this->expand();
527        other2.expand();
528        mask2.expand();
529      }
530      else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
531      {
532        this->tag();
533        other2.tag();
534        mask2.tag();
535      }
536      else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
537      {
538      }
539      else
540      {
541        throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
542      }
543      // Now we iterate over the elements
544      DataVector& self=getReadyPtr()->getVector();
545      const DataVector& ovec=other2.getReadyPtr()->getVector();
546      const DataVector& mvec=mask2.getReadyPtr()->getVector();
547      if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))
548      {
549        throw DataException("Error - size mismatch in arguments to copyWithMask.");
550      }
551      size_t num_points=self.size();
552    
553    mask1 = mask.wherePositive();    // OPENMP 3.0 allows unsigned loop vars.
554    mask2.copy(mask1);  #if defined(_OPENMP) && (_OPENMP < 200805)
555      long i;
556    #else
557      size_t i;
558    #endif
559      #pragma omp parallel for private(i) schedule(static)
560      for (i=0;i<num_points;++i)
561      {
562        if (mvec[i]>0)
563        {
564           self[i]=ovec[i];
565        }
566      }
567    }
568    
   mask1 *= other;  
   mask2 *= *this;  
   mask2 = *this - mask2;  
569    
   *this = mask1 + mask2;  
 }  
570    
571  bool  bool
572  Data::isExpanded() const  Data::isExpanded() const
# Line 350  Data::isTagged() const Line 582  Data::isTagged() const
582    return (temp!=0);    return (temp!=0);
583  }  }
584    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
585  bool  bool
586  Data::isEmpty() const  Data::isEmpty() const
587  {  {
# Line 366  Data::isConstant() const Line 596  Data::isConstant() const
596    return (temp!=0);    return (temp!=0);
597  }  }
598    
599    bool
600    Data::isLazy() const
601    {
602      return m_data->isLazy();
603    }
604    
605    // at the moment this is synonymous with !isLazy() but that could change
606    bool
607    Data::isReady() const
608    {
609      return (dynamic_cast<DataReady*>(m_data.get())!=0);
610    }
611    
612    
613  void  void
614  Data::setProtection()  Data::setProtection()
615  {  {
616     m_protected=true;     m_protected=true;
617  }  }
618    
619  bool  bool
620  Data::isProtected() const  Data::isProtected() const
621  {  {
622     return m_protected;     return m_protected;
623  }  }
624    
# Line 386  Data::expand() Line 630  Data::expand()
630    if (isConstant()) {    if (isConstant()) {
631      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
632      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
633      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
634      m_data=temp_data;  //     m_data=temp_data;
635        m_data=temp->getPtr();
636    } else if (isTagged()) {    } else if (isTagged()) {
637      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
638      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
639      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
640      m_data=temp_data;  //     m_data=temp_data;
641        m_data=temp->getPtr();
642    } else if (isExpanded()) {    } else if (isExpanded()) {
643      //      //
644      // do nothing      // do nothing
645    } else if (isEmpty()) {    } else if (isEmpty()) {
646      throw DataException("Error - Expansion of DataEmpty not possible.");      throw DataException("Error - Expansion of DataEmpty not possible.");
647      } else if (isLazy()) {
648        resolve();
649        expand();       // resolve might not give us expanded data
650    } else {    } else {
651      throw DataException("Error - Expansion not implemented for this Data type.");      throw DataException("Error - Expansion not implemented for this Data type.");
652    }    }
# Line 409  Data::tag() Line 658  Data::tag()
658    if (isConstant()) {    if (isConstant()) {
659      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
660      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
661      shared_ptr<DataAbstract> temp_data(temp);  //     shared_ptr<DataAbstract> temp_data(temp);
662      m_data=temp_data;  //     m_data=temp_data;
663        m_data=temp->getPtr();
664    } else if (isTagged()) {    } else if (isTagged()) {
665      // do nothing      // do nothing
666    } else if (isExpanded()) {    } else if (isExpanded()) {
667      throw DataException("Error - Creating tag data from DataExpanded not possible.");      throw DataException("Error - Creating tag data from DataExpanded not possible.");
668    } else if (isEmpty()) {    } else if (isEmpty()) {
669      throw DataException("Error - Creating tag data from DataEmpty not possible.");      throw DataException("Error - Creating tag data from DataEmpty not possible.");
670      } else if (isLazy()) {
671         DataAbstract_ptr res=m_data->resolve();
672         if (m_data->isExpanded())
673         {
674        throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
675         }
676         m_data=res;    
677         tag();
678    } else {    } else {
679      throw DataException("Error - Tagging not implemented for this Data type.");      throw DataException("Error - Tagging not implemented for this Data type.");
680    }    }
681  }  }
682    
683  void  void
684  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::resolve()
685  {  {
686    m_data->reshapeDataPoint(shape);    if (isLazy())
687      {
688         m_data=m_data->resolve();
689      }
690    }
691    
692    
693    Data
694    Data::oneOver() const
695    {
696      if (isLazy())
697      {
698        DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);
699        return Data(c);
700      }
701      return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
702  }  }
703    
704  Data  Data
705  Data::wherePositive() const  Data::wherePositive() const
706  {  {
707  #if defined DOPROF    if (isLazy())
708    profData->where++;    {
709  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),GZ);
710    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));      return Data(c);
711      }
712      return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
713  }  }
714    
715  Data  Data
716  Data::whereNegative() const  Data::whereNegative() const
717  {  {
718  #if defined DOPROF    if (isLazy())
719    profData->where++;    {
720  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),LZ);
721    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));      return Data(c);
722      }
723      return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
724  }  }
725    
726  Data  Data
727  Data::whereNonNegative() const  Data::whereNonNegative() const
728  {  {
729  #if defined DOPROF    if (isLazy())
730    profData->where++;    {
731  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);
732    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));      return Data(c);
733      }
734      return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
735  }  }
736    
737  Data  Data
738  Data::whereNonPositive() const  Data::whereNonPositive() const
739  {  {
740  #if defined DOPROF    if (isLazy())
741    profData->where++;    {
742  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);
743    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));      return Data(c);
744      }
745      return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
746  }  }
747    
748  Data  Data
749  Data::whereZero(double tol) const  Data::whereZero(double tol) const
750  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
751    Data dataAbs=abs();    Data dataAbs=abs();
752    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));    return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
753  }  }
754    
755  Data  Data
756  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
757  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
758    Data dataAbs=abs();    Data dataAbs=abs();
759    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));    return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
760  }  }
761    
762  Data  Data
763  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
764  {  {
 #if defined DOPROF  
   profData->interpolate++;  
 #endif  
765    return Data(*this,functionspace);    return Data(*this,functionspace);
766  }  }
767    
768  bool  bool
769  Data::probeInterpolation(const FunctionSpace& functionspace) const  Data::probeInterpolation(const FunctionSpace& functionspace) const
770  {  {
771    if (getFunctionSpace()==functionspace) {    return getFunctionSpace().probeInterpolation(functionspace);
772      return true;  //   if (getFunctionSpace()==functionspace) {
773    } else {  //     return true;
774      const AbstractDomain& domain=getDomain();  //   } else {
775      if  (domain==functionspace.getDomain()) {  //     const_Domain_ptr domain=getDomain();
776        return domain.probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());  //     if  (*domain==*functionspace.getDomain()) {
777      } else {  //       return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
778        return domain.probeInterpolationACross(getFunctionSpace().getTypeCode(),functionspace.getDomain(),functionspace.getTypeCode());  //     } else {
779      }  //       return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
780    }  //     }
781    //   }
782  }  }
783    
784  Data  Data
785  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
786  {  {
787  #if defined DOPROF    if (isEmpty())
788    profData->grad++;    {
789  #endif      throw DataException("Error - operation not permitted on instances of DataEmpty.");
790      }
791      double blocktimer_start = blocktimer_time();
792    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
793      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
794    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataTypes::ShapeType grad_shape=getDataPointShape();
795    grad_shape.push_back(functionspace.getDim());    grad_shape.push_back(functionspace.getDim());
796    Data out(0.0,grad_shape,functionspace,true);    Data out(0.0,grad_shape,functionspace,true);
797    getDomain().setToGradient(out,*this);    getDomain()->setToGradient(out,*this);
798      blocktimer_increment("grad()", blocktimer_start);
799    return out;    return out;
800  }  }
801    
802  Data  Data
803  Data::grad() const  Data::grad() const
804  {  {
805    return gradOn(escript::function(getDomain()));    if (isEmpty())
806      {
807        throw DataException("Error - operation not permitted on instances of DataEmpty.");
808      }
809      return gradOn(escript::function(*getDomain()));
810  }  }
811    
812  int  int
813  Data::getDataPointSize() const  Data::getDataPointSize() const
814  {  {
815    return getPointDataView().noValues();    return m_data->getNoValues();
816  }  }
817    
818  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
819  Data::getLength() const  Data::getLength() const
820  {  {
821    return m_data->getLength();    return m_data->getLength();
822  }  }
823    
 const DataArrayView::ShapeType&  
 Data::getDataPointShape() const  
 {  
   return getPointDataView().getShape();  
 }  
   
 void  
 Data::fillFromNumArray(const boost::python::numeric::array num_array)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
   //  
   // check rank  
   if (num_array.getrank()<getDataPointRank())  
       throw DataException("Rank of numarray does not match Data object rank");  
   
   //  
   // check shape of num_array  
   for (int i=0; i<getDataPointRank(); i++) {  
     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])  
        throw DataException("Shape of numarray does not match Data object rank");  
   }  
   
   //  
   // make sure data is expanded:  
   if (!isExpanded()) {  
     expand();  
   }  
   
   //  
   // and copy over  
   m_data->copyAll(num_array);  
 }  
   
824  const  const
825  boost::python::numeric::array  boost::python::numeric::array
826  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
827  {  {
828    //    int i, j, k, l;
829    // determine the total number of data points  
830    int numSamples = getNumSamples();    FORCERESOLVE;
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
831    
832    //    //
833    // determine the rank and shape of each data point    // determine the rank and shape of each data point
834    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
835    DataArrayView::ShapeType dataPointShape = getDataPointShape();    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
836    
837    //    //
838    // create the numeric array to be returned    // create the numeric array to be returned
839    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
840    
841    //    //
842    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
843    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
844    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
845    // data points, whilst the first dimension will be equal to the    const DataTypes::ShapeType& arrayShape = dataPointShape;
   // 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]);  
   }  
846    
847    //    //
848    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
849      if (arrayRank==0) {
850        numArray.resize(1);
851      }
852    if (arrayRank==1) {    if (arrayRank==1) {
853      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
854    }    }
# Line 623  Data::convertToNumArray() Line 861  Data::convertToNumArray()
861    if (arrayRank==4) {    if (arrayRank==4) {
862      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
863    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
864    
865    //    if (getNumDataPointsPerSample()>0) {
866    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
867    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
868    int dataPoint = 0;         //
869    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
870      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
871        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
872        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++;  
     }  
   }  
873    
874           //
875           // Check a valid data point number has been supplied
876           if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
877               throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
878           }
879           // TODO: global error handling
880           // create a view of the data if it is stored locally
881    //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
882           DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
883    
884    
885           switch( dataPointRank ){
886                case 0 :
887                    numArray[0] = getDataAtOffset(offset);
888                    break;
889                case 1 :
890                    for( i=0; i<dataPointShape[0]; i++ )
891                        numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
892                    break;
893                case 2 :
894                    for( i=0; i<dataPointShape[0]; i++ )
895                        for( j=0; j<dataPointShape[1]; j++)
896                            numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
897                    break;
898                case 3 :
899                    for( i=0; i<dataPointShape[0]; i++ )
900                        for( j=0; j<dataPointShape[1]; j++ )
901                            for( k=0; k<dataPointShape[2]; k++)
902                                numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
903                    break;
904                case 4 :
905                    for( i=0; i<dataPointShape[0]; i++ )
906                        for( j=0; j<dataPointShape[1]; j++ )
907                            for( k=0; k<dataPointShape[2]; k++ )
908                                for( l=0; l<dataPointShape[3]; l++)
909                                    numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
910                    break;
911        }
912      }
913    //    //
914    // return the loaded array    // return the array
915    return numArray;    return numArray;
916    
917  }  }
918    
919  const  void
920  boost::python::numeric::array  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
 Data::convertToNumArrayFromSampleNo(int sampleNo)  
921  {  {
922    //      // this will throw if the value cannot be represented
923    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
924    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
925      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();  
926    
927    void
928    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
929    {
930      if (isProtected()) {
931            throw DataException("Error - attempt to update protected Data object.");
932      }
933      FORCERESOLVE;
934    //    //
935    // create the numeric array to be returned    // check rank
936    boost::python::numeric::array numArray(0.0);    if (static_cast<unsigned int>(num_array.getrank())<getDataPointRank())
937          throw DataException("Rank of numarray does not match Data object rank");
938    
939    //    //
940    // the rank of the returned numeric array will be the rank of    // check shape of num_array
941    // the data points, plus one. Where the rank of the array is n,    for (unsigned int i=0; i<getDataPointRank(); i++) {
942    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
943    // data points, whilst the first dimension will be equal to the         throw DataException("Shape of numarray does not match Data object rank");
   // 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(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
944    }    }
   
945    //    //
946    // resize the numeric array to the shape just calculated    // make sure data is expanded:
947    if (arrayRank==1) {    //
948      numArray.resize(arrayShape[0]);    if (!isExpanded()) {
949    }      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]);  
950    }    }
951    if (arrayRank==5) {    if (getNumDataPointsPerSample()>0) {
952      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);         int sampleNo = dataPointNo/getNumDataPointsPerSample();
953           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
954           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
955      } else {
956           m_data->copyToDataPoint(-1, 0,num_array);
957    }    }
958    }
959    
960    //  void
961    // loop through each data point in turn, loading the values for that data point  Data::setValueOfDataPoint(int dataPointNo, const double value)
962    // into the numeric array.  {
963    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    if (isProtected()) {
964      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);          throw DataException("Error - attempt to update protected Data object.");
     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);  
             }  
           }  
         }  
       }  
     }  
965    }    }
   
966    //    //
967    // return the loaded array    // make sure data is expanded:
968    return numArray;    FORCERESOLVE;
969      if (!isExpanded()) {
970        expand();
971      }
972      if (getNumDataPointsPerSample()>0) {
973           int sampleNo = dataPointNo/getNumDataPointsPerSample();
974           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
975           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
976      } else {
977           m_data->copyToDataPoint(-1, 0,value);
978      }
979  }  }
980    
981  const  const
982  boost::python::numeric::array  boost::python::numeric::array
983  Data::convertToNumArrayFromDPNo(int procNo,  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
                                 int sampleNo,  
                                                                 int dataPointNo)  
   
984  {  {
985      size_t length=0;    size_t length=0;
986      int i, j, k, l, pos;    int i, j, k, l, pos;
987      FORCERESOLVE;
   //  
   // 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.");  
   }  
   
988    //    //
989    // determine the rank and shape of each data point    // determine the rank and shape of each data point
990    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
991    DataArrayView::ShapeType dataPointShape = getDataPointShape();    const DataTypes::ShapeType& dataPointShape = getDataPointShape();
992    
993    //    //
994    // create the numeric array to be returned    // create the numeric array to be returned
# Line 815  Data::convertToNumArrayFromDPNo(int proc Line 998  Data::convertToNumArrayFromDPNo(int proc
998    // the shape of the returned numeric array will be the same    // the shape of the returned numeric array will be the same
999    // as that of the data point    // as that of the data point
1000    int arrayRank = dataPointRank;    int arrayRank = dataPointRank;
1001    DataArrayView::ShapeType arrayShape = dataPointShape;    const DataTypes::ShapeType& arrayShape = dataPointShape;
1002    
1003    //    //
1004    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 1018  Data::convertToNumArrayFromDPNo(int proc
1018      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
1019    }    }
1020    
1021      // added for the MPI communication    // added for the MPI communication
1022      length=1;    length=1;
1023      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
1024          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
1025    
1026    //    //
1027    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
1028    
1029      // updated for the MPI case      // updated for the MPI case
1030      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
1031                 if (getNumDataPointsPerSample()>0) {
1032                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
1033                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1034                    //
1035                    // Check a valid sample number has been supplied
1036                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1037                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
1038                    }
1039    
1040                    //
1041                    // Check a valid data point number has been supplied
1042                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1043                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
1044                    }
1045                    // TODO: global error handling
1046          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
1047          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
1048                    DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1049    
1050          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
1051          pos=0;          pos=0;
1052          switch( dataPointRank ){          switch( dataPointRank ){
1053              case 0 :              case 0 :
1054                  tmpData[0] = dataPointView();                  tmpData[0] = getDataAtOffset(offset);
1055                  break;                  break;
1056              case 1 :                      case 1 :
1057                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
1058                      tmpData[i]=dataPointView(i);                      tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
1059                  break;                  break;
1060              case 2 :                      case 2 :
1061                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
1062                      for( j=0; j<dataPointShape[1]; j++, pos++ )                      for( j=0; j<dataPointShape[1]; j++, pos++ )
1063                          tmpData[pos]=dataPointView(i,j);                          tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
1064                  break;                  break;
1065              case 3 :                      case 3 :
1066                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
1067                      for( j=0; j<dataPointShape[1]; j++ )                      for( j=0; j<dataPointShape[1]; j++ )
1068                          for( k=0; k<dataPointShape[2]; k++, pos++ )                          for( k=0; k<dataPointShape[2]; k++, pos++ )
1069                              tmpData[pos]=dataPointView(i,j,k);                              tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
1070                  break;                  break;
1071              case 4 :              case 4 :
1072                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
1073                      for( j=0; j<dataPointShape[1]; j++ )                      for( j=0; j<dataPointShape[1]; j++ )
1074                          for( k=0; k<dataPointShape[2]; k++ )                          for( k=0; k<dataPointShape[2]; k++ )
1075                              for( l=0; l<dataPointShape[3]; l++, pos++ )                              for( l=0; l<dataPointShape[3]; l++, pos++ )
1076                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
1077                  break;                  break;
1078          }          }
1079                }
1080      }      }
1081  #ifdef PASO_MPI          #ifdef PASO_MPI
1082          // broadcast the data to all other processes          // broadcast the data to all other processes
1083          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1084  #endif          #endif
1085    
1086      // unpack the data      // unpack the data
1087      switch( dataPointRank ){      switch( dataPointRank ){
1088          case 0 :          case 0 :
1089              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
1090              break;              break;
1091          case 1 :                  case 1 :
1092              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
1093                  numArray[i]=tmpData[i];                  numArray[i]=tmpData[i];
1094              break;              break;
1095          case 2 :                  case 2 :
1096              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
1097                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
1098                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1099              break;              break;
1100          case 3 :                  case 3 :
1101              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
1102                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
1103                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
1104                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1105              break;              break;
1106          case 4 :          case 4 :
1107              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
1108                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
1109                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
1110                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
1111                              tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];                                  numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1112              break;              break;
1113      }      }
1114    
1115      delete [] tmpData;        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);  
           }  
         }  
       }  
     }  
   }  
 */  
   
1116    //    //
1117    // return the loaded array    // return the loaded array
1118    return numArray;    return numArray;
1119  }  }
1120    
1121    
1122    boost::python::numeric::array
1123    Data::integrate_const() const
1124    {
1125      if (isLazy())
1126      {
1127        throw DataException("Error - cannot integrate for constant lazy data.");
1128      }
1129      return integrateWorker();
1130    }
1131    
1132  boost::python::numeric::array  boost::python::numeric::array
1133  Data::integrate() const  Data::integrate()
1134    {
1135      if (isLazy())
1136      {
1137        expand();
1138      }
1139      return integrateWorker();
1140    }
1141    
1142    
1143    
1144    boost::python::numeric::array
1145    Data::integrateWorker() const
1146  {  {
1147    int index;    int index;
1148    int rank = getDataPointRank();    int rank = getDataPointRank();
1149    DataArrayView::ShapeType shape = getDataPointShape();    DataTypes::ShapeType shape = getDataPointShape();
1150      int dataPointSize = getDataPointSize();
 #if defined DOPROF  
   profData->integrate++;  
 #endif  
1151    
1152    //    //
1153    // calculate the integral values    // calculate the integral values
1154    vector<double> integrals(getDataPointSize());    vector<double> integrals(dataPointSize);
1155    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);    vector<double> integrals_local(dataPointSize);
1156    #ifdef PASO_MPI
1157      AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);
1158      // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1159      double *tmp = new double[dataPointSize];
1160      double *tmp_local = new double[dataPointSize];
1161      for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1162      MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1163      for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1164      delete[] tmp;
1165      delete[] tmp_local;
1166    #else
1167      AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1168    #endif
1169    
1170    //    //
1171    // create the numeric array to be returned    // create the numeric array to be returned
# Line 1031  Data::integrate() const Line 1225  Data::integrate() const
1225  Data  Data
1226  Data::sin() const  Data::sin() const
1227  {  {
1228  #if defined DOPROF    if (isLazy())
1229    profData->unary++;    {
1230  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),SIN);
1231    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);      return Data(c);
1232      }
1233      return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1234  }  }
1235    
1236  Data  Data
1237  Data::cos() const  Data::cos() const
1238  {  {
1239  #if defined DOPROF    if (isLazy())
1240    profData->unary++;    {
1241  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),COS);
1242    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);      return Data(c);
1243      }
1244      return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1245  }  }
1246    
1247  Data  Data
1248  Data::tan() const  Data::tan() const
1249  {  {
1250  #if defined DOPROF    if (isLazy())
1251    profData->unary++;    {
1252  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),TAN);
1253    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);      return Data(c);
1254      }
1255      return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1256  }  }
1257    
1258  Data  Data
1259  Data::asin() const  Data::asin() const
1260  {  {
1261  #if defined DOPROF    if (isLazy())
1262    profData->unary++;    {
1263  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);
1264    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);      return Data(c);
1265      }
1266      return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1267  }  }
1268    
1269  Data  Data
1270  Data::acos() const  Data::acos() const
1271  {  {
1272  #if defined DOPROF    if (isLazy())
1273    profData->unary++;    {
1274  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);
1275    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);      return Data(c);
1276      }
1277      return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1278  }  }
1279    
1280    
1281  Data  Data
1282  Data::atan() const  Data::atan() const
1283  {  {
1284  #if defined DOPROF    if (isLazy())
1285    profData->unary++;    {
1286  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);
1287    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);      return Data(c);
1288      }
1289      return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1290  }  }
1291    
1292  Data  Data
1293  Data::sinh() const  Data::sinh() const
1294  {  {
1295  #if defined DOPROF    if (isLazy())
1296    profData->unary++;    {
1297  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),SINH);
1298    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);      return Data(c);
1299      }
1300      return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1301  }  }
1302    
1303  Data  Data
1304  Data::cosh() const  Data::cosh() const
1305  {  {
1306  #if defined DOPROF    if (isLazy())
1307    profData->unary++;    {
1308  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),COSH);
1309    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);      return Data(c);
1310      }
1311      return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1312  }  }
1313    
1314  Data  Data
1315  Data::tanh() const  Data::tanh() const
1316  {  {
1317  #if defined DOPROF    if (isLazy())
1318    profData->unary++;    {
1319        DataLazy* c=new DataLazy(borrowDataPtr(),TANH);
1320        return Data(c);
1321      }
1322      return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1323    }
1324    
1325    
1326    Data
1327    Data::erf() const
1328    {
1329    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1330      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1331    #else
1332      if (isLazy())
1333      {
1334        DataLazy* c=new DataLazy(borrowDataPtr(),ERF);
1335        return Data(c);
1336      }
1337      return C_TensorUnaryOperation(*this, ::erf);
1338  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);  
1339  }  }
1340    
1341  Data  Data
1342  Data::asinh() const  Data::asinh() const
1343  {  {
1344  #if defined DOPROF    if (isLazy())
1345    profData->unary++;    {
1346        DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1347        return Data(c);
1348      }
1349    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1350      return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1351    #else
1352      return C_TensorUnaryOperation(*this, ::asinh);
1353  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);  
1354  }  }
1355    
1356  Data  Data
1357  Data::acosh() const  Data::acosh() const
1358  {  {
1359  #if defined DOPROF    if (isLazy())
1360    profData->unary++;    {
1361        DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1362        return Data(c);
1363      }
1364    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1365      return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1366    #else
1367      return C_TensorUnaryOperation(*this, ::acosh);
1368  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);  
1369  }  }
1370    
1371  Data  Data
1372  Data::atanh() const  Data::atanh() const
1373  {  {
1374  #if defined DOPROF    if (isLazy())
1375    profData->unary++;    {
1376        DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1377        return Data(c);
1378      }
1379    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1380      return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1381    #else
1382      return C_TensorUnaryOperation(*this, ::atanh);
1383  #endif  #endif
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);  
1384  }  }
1385    
1386  Data  Data
1387  Data::log10() const  Data::log10() const
1388  {  {  if (isLazy())
1389  #if defined DOPROF    {
1390    profData->unary++;      DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);
1391  #endif      return Data(c);
1392    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    }
1393      return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1394  }  }
1395    
1396  Data  Data
1397  Data::log() const  Data::log() const
1398  {  {
1399  #if defined DOPROF    if (isLazy())
1400    profData->unary++;    {
1401  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),LOG);
1402    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);      return Data(c);
1403      }
1404      return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1405  }  }
1406    
1407  Data  Data
1408  Data::sign() const  Data::sign() const
1409  {  {
1410  #if defined DOPROF    if (isLazy())
1411    profData->unary++;    {
1412  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);
1413    return escript::unaryOp(*this,escript::fsign);      return Data(c);
1414      }
1415      return C_TensorUnaryOperation(*this, escript::fsign);
1416  }  }
1417    
1418  Data  Data
1419  Data::abs() const  Data::abs() const
1420  {  {
1421  #if defined DOPROF    if (isLazy())
1422    profData->unary++;    {
1423  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),ABS);
1424    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);      return Data(c);
1425      }
1426      return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1427  }  }
1428    
1429  Data  Data
1430  Data::neg() const  Data::neg() const
1431  {  {
1432  #if defined DOPROF    if (isLazy())
1433    profData->unary++;    {
1434  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),NEG);
1435    return escript::unaryOp(*this,negate<double>());      return Data(c);
1436      }
1437      return C_TensorUnaryOperation(*this, negate<double>());
1438  }  }
1439    
1440  Data  Data
1441  Data::pos() const  Data::pos() const
1442  {  {
1443  #if defined DOPROF      // not doing lazy check here is deliberate.
1444    profData->unary++;      // since a deep copy of lazy data should be cheap, I'll just let it happen now
 #endif  
1445    Data result;    Data result;
1446    // perform a deep copy    // perform a deep copy
1447    result.copy(*this);    result.copy(*this);
# Line 1195  Data::pos() const Line 1450  Data::pos() const
1450    
1451  Data  Data
1452  Data::exp() const  Data::exp() const
1453  {  {  
1454  #if defined DOPROF    if (isLazy())
1455    profData->unary++;    {
1456  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),EXP);
1457    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);      return Data(c);
1458      }
1459      return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1460  }  }
1461    
1462  Data  Data
1463  Data::sqrt() const  Data::sqrt() const
1464  {  {
1465  #if defined DOPROF    if (isLazy())
1466    profData->unary++;    {
1467  #endif      DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);
1468    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);      return Data(c);
1469      }
1470      return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1471  }  }
1472    
1473  double  double
1474  Data::Lsup() const  Data::Lsup_const() const
1475  {  {
1476    double localValue, globalValue;     if (isLazy())
1477  #if defined DOPROF     {
1478    profData->reduction1++;      throw DataException("Error - cannot compute Lsup for constant lazy data.");
1479  #endif     }
1480    //     return LsupWorker();
1481    // set the initial absolute maximum value to zero  }
1482    
1483    AbsMax abs_max_func;  double
1484    localValue = algorithm(abs_max_func,0);  Data::Lsup()
1485  #ifdef PASO_MPI  {
1486    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );     if (isLazy())
1487    return globalValue;     {
1488  #else      expand();
1489    return localValue;     }
1490  #endif     return LsupWorker();
1491  }  }
1492    
1493  double  double
1494  Data::Linf() const  Data::sup_const() const
1495  {  {
1496    double localValue, globalValue;     if (isLazy())
1497  #if defined DOPROF     {
1498    profData->reduction1++;      throw DataException("Error - cannot compute sup for constant lazy data.");
1499  #endif     }
1500       return supWorker();
1501    }
1502    
1503    double
1504    Data::sup()
1505    {
1506       if (isLazy())
1507       {
1508        expand();
1509       }
1510       return supWorker();
1511    }
1512    
1513    double
1514    Data::inf_const() const
1515    {
1516       if (isLazy())
1517       {
1518        throw DataException("Error - cannot compute inf for constant lazy data.");
1519       }
1520       return infWorker();
1521    }
1522    
1523    double
1524    Data::inf()
1525    {
1526       if (isLazy())
1527       {
1528        expand();
1529       }
1530       return infWorker();
1531    }
1532    
1533    double
1534    Data::LsupWorker() const
1535    {
1536      double localValue;
1537    //    //
1538    // 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());  
1539    
1540      AbsMax abs_max_func;
1541      localValue = algorithm(abs_max_func,0);
1542  #ifdef PASO_MPI  #ifdef PASO_MPI
1543    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );    double globalValue;
1544      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1545    return globalValue;    return globalValue;
1546  #else  #else
1547    return localValue;    return localValue;
# Line 1252  Data::Linf() const Line 1549  Data::Linf() const
1549  }  }
1550    
1551  double  double
1552  Data::sup() const  Data::supWorker() const
1553  {  {
1554    double localValue, globalValue;    double localValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1555    //    //
1556    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1557    FMax fmax_func;    FMax fmax_func;
1558    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1559  #ifdef PASO_MPI  #ifdef PASO_MPI
1560      double globalValue;
1561    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1562    return globalValue;    return globalValue;
1563  #else  #else
# Line 1271  Data::sup() const Line 1566  Data::sup() const
1566  }  }
1567    
1568  double  double
1569  Data::inf() const  Data::infWorker() const
1570  {  {
1571    double localValue, globalValue;    double localValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1572    //    //
1573    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1574    FMin fmin_func;    FMin fmin_func;
1575    localValue = algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1576  #ifdef PASO_MPI  #ifdef PASO_MPI
1577      double globalValue;
1578    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );    MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1579    return globalValue;    return globalValue;
1580  #else  #else
# Line 1294  Data::inf() const Line 1587  Data::inf() const
1587  Data  Data
1588  Data::maxval() const  Data::maxval() const
1589  {  {
1590  #if defined DOPROF    if (isLazy())
1591    profData->reduction2++;    {
1592  #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1593        temp.resolve();
1594        return temp.maxval();
1595      }
1596    //    //
1597    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1598    FMax fmax_func;    FMax fmax_func;
# Line 1306  Data::maxval() const Line 1602  Data::maxval() const
1602  Data  Data
1603  Data::minval() const  Data::minval() const
1604  {  {
1605  #if defined DOPROF    if (isLazy())
1606    profData->reduction2++;    {
1607  #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1608        temp.resolve();
1609        return temp.minval();
1610      }
1611    //    //
1612    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1613    FMin fmin_func;    FMin fmin_func;
# Line 1316  Data::minval() const Line 1615  Data::minval() const
1615  }  }
1616    
1617  Data  Data
1618  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1619  {  {
1620  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1621    profData->reduction2++;       DataTypes::ShapeType s=getDataPointShape();
1622  #endif       DataTypes::ShapeType ev_shape;
1623    Trace trace_func;       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1624    return dp_algorithm(trace_func,0);       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1625         int rank=getDataPointRank();
1626         if (rank<2) {
1627            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1628         }
1629         if (axis0<0 || axis0>rank-1) {
1630            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1631         }
1632         if (axis1<0 || axis1>rank-1) {
1633             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1634         }
1635         if (axis0 == axis1) {
1636             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1637         }
1638         if (axis0 > axis1) {
1639             axis0_tmp=axis1;
1640             axis1_tmp=axis0;
1641         } else {
1642             axis0_tmp=axis0;
1643             axis1_tmp=axis1;
1644         }
1645         for (int i=0; i<rank; i++) {
1646           if (i == axis0_tmp) {
1647              ev_shape.push_back(s[axis1_tmp]);
1648           } else if (i == axis1_tmp) {
1649              ev_shape.push_back(s[axis0_tmp]);
1650           } else {
1651              ev_shape.push_back(s[i]);
1652           }
1653         }
1654         Data ev(0.,ev_shape,getFunctionSpace());
1655         ev.typeMatchRight(*this);
1656         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1657         return ev;
1658    
1659  }  }
1660    
1661  Data  Data
1662  Data::symmetric() const  Data::symmetric() const
1663  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1664       // check input       // check input
1665       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1666       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1667          if(s[0] != s[1])          if(s[0] != s[1])
1668             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.");
1669       }       }
1670       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
# Line 1344  Data::symmetric() const Line 1674  Data::symmetric() const
1674       else {       else {
1675          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.");
1676       }       }
1677         if (isLazy())
1678         {
1679        DataLazy* c=new DataLazy(borrowDataPtr(),SYM);
1680        return Data(c);
1681         }
1682       Data ev(0.,getDataPointShape(),getFunctionSpace());       Data ev(0.,getDataPointShape(),getFunctionSpace());
1683       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1684       m_data->symmetric(ev.m_data.get());       m_data->symmetric(ev.m_data.get());
# Line 1353  Data::symmetric() const Line 1688  Data::symmetric() const
1688  Data  Data
1689  Data::nonsymmetric() const  Data::nonsymmetric() const
1690  {  {
1691       #if defined DOPROF       if (isLazy())
1692          profData->unary++;       {
1693       #endif      DataLazy* c=new DataLazy(borrowDataPtr(),NSYM);
1694        return Data(c);
1695         }
1696       // check input       // check input
1697       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1698       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1699          if(s[0] != s[1])          if(s[0] != s[1])
1700             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.");
1701          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1702          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
1703          ev_shape.push_back(s[1]);          ev_shape.push_back(s[1]);
1704          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
# Line 1372  Data::nonsymmetric() const Line 1709  Data::nonsymmetric() const
1709       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
1710          if(!(s[0] == s[2] && s[1] == s[3]))          if(!(s[0] == s[2] && s[1] == s[3]))
1711             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.");
1712          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1713          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
1714          ev_shape.push_back(s[1]);          ev_shape.push_back(s[1]);
1715          ev_shape.push_back(s[2]);          ev_shape.push_back(s[2]);
# Line 1387  Data::nonsymmetric() const Line 1724  Data::nonsymmetric() const
1724       }       }
1725  }  }
1726    
1727  Data  
1728  Data::matrixtrace(int axis_offset) const  // Doing a lazy version of this would require some thought.
1729  {  // First it needs a parameter (which DataLazy doesn't support at the moment).
1730       #if defined DOPROF  // (secondly although it does not apply to trace) we can't handle operations which return
1731          profData->unary++;  // multiple results (like eigenvectors_values) or return values of different shapes to their input
1732       #endif  // (like eigenvalues).
1733       DataArrayView::ShapeType s=getDataPointShape();  Data
1734    Data::trace(int axis_offset) const
1735    {
1736         if (isLazy())
1737         {
1738        Data temp(*this);   // to get around the fact that you can't resolve a const Data
1739        temp.resolve();
1740        return temp.trace(axis_offset);
1741         }
1742         DataTypes::ShapeType s=getDataPointShape();
1743       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1744          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1745          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1746          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1747          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1748          return ev;          return ev;
1749       }       }
1750       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
1751          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1752          if (axis_offset==0) {          if (axis_offset==0) {
1753            int s2=s[2];            int s2=s[2];
1754            ev_shape.push_back(s2);            ev_shape.push_back(s2);
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1759  Data::matrixtrace(int axis_offset) const
1759          }          }
1760          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1761          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1762          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1763          return ev;          return ev;
1764       }       }
1765       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
1766          DataArrayView::ShapeType ev_shape;          DataTypes::ShapeType ev_shape;
1767          if (axis_offset==0) {          if (axis_offset==0) {
1768            ev_shape.push_back(s[2]);            ev_shape.push_back(s[2]);
1769            ev_shape.push_back(s[3]);            ev_shape.push_back(s[3]);
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1778  Data::matrixtrace(int axis_offset) const
1778      }      }
1779          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1780          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1781      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1782          return ev;          return ev;
1783       }       }
1784       else {       else {
1785          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.");
1786       }       }
1787  }  }
1788    
1789  Data  Data
1790  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1791  {  {    
1792  #if defined DOPROF       if (isLazy())
1793       profData->reduction2++;       {
1794  #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1795       DataArrayView::ShapeType s=getDataPointShape();      temp.resolve();
1796       DataArrayView::ShapeType ev_shape;      return temp.transpose(axis_offset);
1797         }
1798         DataTypes::ShapeType s=getDataPointShape();
1799         DataTypes::ShapeType ev_shape;
1800       // 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]
1801       // 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)
1802       int rank=getDataPointRank();       int rank=getDataPointRank();
# Line 1467  Data::transpose(int axis_offset) const Line 1816  Data::transpose(int axis_offset) const
1816  Data  Data
1817  Data::eigenvalues() const  Data::eigenvalues() const
1818  {  {
1819       #if defined DOPROF       if (isLazy())
1820          profData->unary++;       {
1821       #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1822        temp.resolve();
1823        return temp.eigenvalues();
1824         }
1825       // check input       // check input
1826       DataArrayView::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1827       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1828          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.");
1829       if(s[0] != s[1])       if(s[0] != s[1])
1830          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.");
1831       // create return       // create return
1832       DataArrayView::ShapeType ev_shape(1,s[0]);       DataTypes::ShapeType ev_shape(1,s[0]);
1833       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1834       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1835       m_data->eigenvalues(ev.m_data.get());       m_data->eigenvalues(ev.m_data.get());
# Line 1487  Data::eigenvalues() const Line 1839  Data::eigenvalues() const
1839  const boost::python::tuple  const boost::python::tuple
1840  Data::eigenvalues_and_eigenvectors(const double tol) const  Data::eigenvalues_and_eigenvectors(const double tol) const
1841  {  {
1842       #if defined DOPROF       if (isLazy())
1843          profData->unary++;       {
1844       #endif      Data temp(*this);   // to get around the fact that you can't resolve a const Data
1845       DataArrayView::ShapeType s=getDataPointShape();      temp.resolve();
1846       if (getDataPointRank()!=2)      return temp.eigenvalues_and_eigenvectors(tol);
1847         }
1848         DataTypes::ShapeType s=getDataPointShape();
1849         if (getDataPointRank()!=2)
1850          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.");
1851       if(s[0] != s[1])       if(s[0] != s[1])
1852          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.");
1853       // create return       // create return
1854       DataArrayView::ShapeType ev_shape(1,s[0]);       DataTypes::ShapeType ev_shape(1,s[0]);
1855       Data ev(0.,ev_shape,getFunctionSpace());       Data ev(0.,ev_shape,getFunctionSpace());
1856       ev.typeMatchRight(*this);       ev.typeMatchRight(*this);
1857       DataArrayView::ShapeType V_shape(2,s[0]);       DataTypes::ShapeType V_shape(2,s[0]);
1858       Data V(0.,V_shape,getFunctionSpace());       Data V(0.,V_shape,getFunctionSpace());
1859       V.typeMatchRight(*this);       V.typeMatchRight(*this);
1860       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 1862  Data::eigenvalues_and_eigenvectors(const
1862  }  }
1863    
1864  const boost::python::tuple  const boost::python::tuple
1865  Data::mindp() const  Data::minGlobalDataPoint() const
1866  {  {
1867    // 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
1868    // 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
1869    // surrounding function    // surrounding function
1870    
   int SampleNo;  
1871    int DataPointNo;    int DataPointNo;
1872      int ProcNo;    int ProcNo;
1873    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1874    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1875  }  }
1876    
1877  void  void
1878  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1879                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1880  {  {
1881      if (isLazy())
1882      {
1883        Data temp(*this);   // to get around the fact that you can't resolve a const Data
1884        temp.resolve();
1885        return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
1886      }
1887    int i,j;    int i,j;
1888    int lowi=0,lowj=0;    int lowi=0,lowj=0;
1889    double min=numeric_limits<double>::max();    double min=numeric_limits<double>::max();
# Line 1535  Data::calc_mindp(  int& ProcNo, Line 1894  Data::calc_mindp(  int& ProcNo,
1894    int numDPPSample=temp.getNumDataPointsPerSample();    int numDPPSample=temp.getNumDataPointsPerSample();
1895    
1896    double next,local_min;    double next,local_min;
1897    int local_lowi,local_lowj;    int local_lowi=0,local_lowj=0;    
1898    
1899    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1900    {    {
# Line 1543  Data::calc_mindp(  int& ProcNo, Line 1902  Data::calc_mindp(  int& ProcNo,
1902      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
1903      for (i=0; i<numSamples; i++) {      for (i=0; i<numSamples; i++) {
1904        for (j=0; j<numDPPSample; j++) {        for (j=0; j<numDPPSample; j++) {
1905          next=temp.getDataPoint(i,j)();          next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1906          if (next<local_min) {          if (next<local_min) {
1907            local_min=next;            local_min=next;
1908            local_lowi=i;            local_lowi=i;
# Line 1561  Data::calc_mindp(  int& ProcNo, Line 1920  Data::calc_mindp(  int& ProcNo,
1920    
1921  #ifdef PASO_MPI  #ifdef PASO_MPI
1922      // determine the processor on which the minimum occurs      // determine the processor on which the minimum occurs
1923      next = temp.getDataPoint(lowi,lowj)();      next = temp.getDataPoint(lowi,lowj);
1924      int lowProc = 0;      int lowProc = 0;
1925      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1926      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1927        
1928      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1929          next = globalMins[lowProc];          next = globalMins[lowProc];
1930          for( i=1; i<get_MPISize(); i++ )          for( i=1; i<get_MPISize(); i++ )
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1940  Data::calc_mindp(  int& ProcNo,
1940  #else  #else
1941      ProcNo = 0;      ProcNo = 0;
1942  #endif  #endif
1943    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1944  }  }
1945    
1946  void  void
1947  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1948  {  {
1949      if (isEmpty())
1950      {
1951        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1952      }
1953      if (isLazy())
1954      {
1955         Data temp(*this);  // to get around the fact that you can't resolve a const Data
1956         temp.resolve();
1957         temp.saveDX(fileName);
1958         return;
1959      }
1960    boost::python::dict args;    boost::python::dict args;
1961    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
1962    getDomain().saveDX(fileName,args);    getDomain()->saveDX(fileName,args);
1963    return;    return;
1964  }  }
1965    
1966  void  void
1967  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1968  {  {
1969      if (isEmpty())
1970      {
1971        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1972      }
1973      if (isLazy())
1974      {
1975         Data temp(*this);  // to get around the fact that you can't resolve a const Data
1976         temp.resolve();
1977         temp.saveVTK(fileName);
1978         return;
1979      }
1980    boost::python::dict args;    boost::python::dict args;
1981    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
1982    getDomain().saveVTK(fileName,args);    getDomain()->saveVTK(fileName,args);
1983    return;    return;
1984  }  }
1985    
# Line 1609  Data::operator+=(const Data& right) Line 1989  Data::operator+=(const Data& right)
1989    if (isProtected()) {    if (isProtected()) {
1990          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1991    }    }
1992  #if defined DOPROF    if (isLazy() || right.isLazy())
1993    profData->binary++;    {
1994  #endif      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
1995    binaryOp(right,plus<double>());          m_data=c->getPtr();
1996    return (*this);      return (*this);
1997      }
1998      else
1999      {
2000        binaryOp(right,plus<double>());
2001        return (*this);
2002      }
2003  }  }
2004    
2005  Data&  Data&
# Line 1622  Data::operator+=(const boost::python::ob Line 2008  Data::operator+=(const boost::python::ob
2008    if (isProtected()) {    if (isProtected()) {
2009          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2010    }    }
2011  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2012    profData->binary++;    if (isLazy())
2013  #endif    {
2014    binaryOp(right,plus<double>());      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD);   // for lazy + is equivalent to +=
2015            m_data=c->getPtr();
2016        return (*this);
2017      }
2018      else
2019      {
2020        binaryOp(tmp,plus<double>());
2021        return (*this);
2022      }
2023    }
2024    
2025    // Hmmm, operator= makes a deep copy but the copy constructor does not?
2026    Data&
2027    Data::operator=(const Data& other)
2028    {
2029      copy(other);
2030    return (*this);    return (*this);
2031  }  }
2032    
# Line 1635  Data::operator-=(const Data& right) Line 2036  Data::operator-=(const Data& right)
2036    if (isProtected()) {    if (isProtected()) {
2037          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2038    }    }
2039  #if defined DOPROF    if (isLazy() || right.isLazy())
2040    profData->binary++;    {
2041  #endif      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2042    binaryOp(right,minus<double>());          m_data=c->getPtr();
2043    return (*this);      return (*this);
2044      }
2045      else
2046      {
2047        binaryOp(right,minus<double>());
2048        return (*this);
2049      }
2050  }  }
2051    
2052  Data&  Data&
# Line 1648  Data::operator-=(const boost::python::ob Line 2055  Data::operator-=(const boost::python::ob
2055    if (isProtected()) {    if (isProtected()) {
2056          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2057    }    }
2058  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2059    profData->binary++;    if (isLazy())
2060  #endif    {
2061    binaryOp(right,minus<double>());      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB);   // for lazy - is equivalent to -=
2062    return (*this);          m_data=c->getPtr();
2063        return (*this);
2064      }
2065      else
2066      {
2067        binaryOp(tmp,minus<double>());
2068        return (*this);
2069      }
2070  }  }
2071    
2072  Data&  Data&
# Line 1661  Data::operator*=(const Data& right) Line 2075  Data::operator*=(const Data& right)
2075    if (isProtected()) {    if (isProtected()) {
2076          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2077    }    }
2078  #if defined DOPROF    if (isLazy() || right.isLazy())
2079    profData->binary++;    {
2080  #endif      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2081    binaryOp(right,multiplies<double>());          m_data=c->getPtr();
2082    return (*this);      return (*this);
2083      }
2084      else
2085      {
2086        binaryOp(right,multiplies<double>());
2087        return (*this);
2088      }
2089  }  }
2090    
2091  Data&  Data&
2092  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
2093  {  {  
2094    if (isProtected()) {    if (isProtected()) {
2095          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2096    }    }
2097  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2098    profData->binary++;    if (isLazy())
2099  #endif    {
2100    binaryOp(right,multiplies<double>());      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL);   // for lazy * is equivalent to *=
2101    return (*this);          m_data=c->getPtr();
2102        return (*this);
2103      }
2104      else
2105      {
2106        binaryOp(tmp,multiplies<double>());
2107        return (*this);
2108      }
2109  }  }
2110    
2111  Data&  Data&
# Line 1687  Data::operator/=(const Data& right) Line 2114  Data::operator/=(const Data& right)
2114    if (isProtected()) {    if (isProtected()) {
2115          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2116    }    }
2117  #if defined DOPROF    if (isLazy() || right.isLazy())
2118    profData->binary++;    {
2119  #endif      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2120    binaryOp(right,divides<double>());          m_data=c->getPtr();
2121    return (*this);      return (*this);
2122      }
2123      else
2124      {
2125        binaryOp(right,divides<double>());
2126        return (*this);
2127      }
2128  }  }
2129    
2130  Data&  Data&
# Line 1700  Data::operator/=(const boost::python::ob Line 2133  Data::operator/=(const boost::python::ob
2133    if (isProtected()) {    if (isProtected()) {
2134          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2135    }    }
2136  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2137    profData->binary++;    if (isLazy())
2138  #endif    {
2139    binaryOp(right,divides<double>());      DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV);   // for lazy / is equivalent to /=
2140    return (*this);          m_data=c->getPtr();
2141        return (*this);
2142      }
2143      else
2144      {
2145        binaryOp(tmp,divides<double>());
2146        return (*this);
2147      }
2148  }  }
2149    
2150  Data  Data
2151  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
2152  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
2153    Data left_d(left,*this);    Data left_d(left,*this);
2154    return left_d.powD(*this);    return left_d.powD(*this);
2155  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 2157  Data::rpowO(const boost::python::object&
2157  Data  Data
2158  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
2159  {  {
2160  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
2161    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
2162  }  }
2163    
2164  Data  Data
2165  Data::powD(const Data& right) const  Data::powD(const Data& right) const
2166  {  {
2167  #if defined DOPROF    if (isLazy() || right.isLazy())
2168    profData->binary++;    {
2169  #endif      DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);
2170    Data result;      return Data(c);
2171    result.copy(*this);    }
2172    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
   return result;  
2173  }  }
2174    
   
2175  //  //
2176  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
2177  Data  Data
2178  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
2179  {  {
2180    Data result;    if (left.isLazy() || right.isLazy())
2181    //    {
2182    // perform a deep copy      DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
2183    result.copy(left);      return Data(c);
2184    result+=right;    }
2185    return result;    return C_TensorBinaryOperation(left, right, plus<double>());
2186  }  }
2187    
2188  //  //
# Line 1763  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;    if (left.isLazy() || right.isLazy())
2194    //    {
2195    // perform a deep copy      DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
2196    result.copy(left);      return Data(c);
2197    result-=right;    }
2198    return result;    return C_TensorBinaryOperation(left, right, minus<double>());
2199  }  }
2200    
2201  //  //
# Line 1776  escript::operator-(const Data& left, con Line 2203  escript::operator-(const Data& left, con
2203  Data  Data
2204  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
2205  {  {
2206    Data result;    if (left.isLazy() || right.isLazy())
2207    //    {
2208    // perform a deep copy      DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
2209    result.copy(left);      return Data(c);
2210    result*=right;    }
2211    return result;    return C_TensorBinaryOperation(left, right, multiplies<double>());
2212  }  }
2213    
2214  //  //
# Line 1789  escript::operator*(const Data& left, con Line 2216  escript::operator*(const Data& left, con
2216  Data  Data
2217  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
2218  {  {
2219    Data result;    if (left.isLazy() || right.isLazy())
2220    //    {
2221    // perform a deep copy      DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);
2222    result.copy(left);      return Data(c);
2223    result/=right;    }
2224    return result;    return C_TensorBinaryOperation(left, right, divides<double>());
2225  }  }
2226    
2227  //  //
# Line 1802  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    //    if (left.isLazy())
2233    // Convert to DataArray format if possible    {
2234    DataArray temp(right);      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);
2235    Data result;      return Data(c);
2236    //    }
2237    // perform a deep copy    return left+Data(right,left.getFunctionSpace(),false);
   result.copy(left);  
   result+=right;  
   return result;  
2238  }  }
2239    
2240  //  //
# Line 1818  escript::operator+(const Data& left, con Line 2242  escript::operator+(const Data& left, con
2242  Data  Data
2243  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
2244  {  {
2245    //    if (left.isLazy())
2246    // Convert to DataArray format if possible    {
2247    DataArray temp(right);      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);
2248    Data result;      return Data(c);
2249    //    }
2250    // perform a deep copy    return left-Data(right,left.getFunctionSpace(),false);
   result.copy(left);  
   result-=right;  
   return result;  
2251  }  }
2252    
2253  //  //
# Line 1834  escript::operator-(const Data& left, con Line 2255  escript::operator-(const Data& left, con
2255  Data  Data
2256  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
2257  {  {
2258    //    if (left.isLazy())
2259    // Convert to DataArray format if possible    {
2260    DataArray temp(right);      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);
2261    Data result;      return Data(c);
2262    //    }
2263    // perform a deep copy    return left*Data(right,left.getFunctionSpace(),false);
   result.copy(left);  
   result*=right;  
   return result;  
2264  }  }
2265    
2266  //  //
# Line 1850  escript::operator*(const Data& left, con Line 2268  escript::operator*(const Data& left, con
2268  Data  Data
2269  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
2270  {  {
2271    //    if (left.isLazy())
2272    // Convert to DataArray format if possible    {
2273    DataArray temp(right);      DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);
2274    Data result;      return Data(c);
2275    //    }
2276    // perform a deep copy    return left/Data(right,left.getFunctionSpace(),false);
   result.copy(left);  
   result/=right;  
   return result;  
2277  }  }
2278    
2279  //  //
# Line 1866  escript::operator/(const Data& left, con Line 2281  escript::operator/(const Data& left, con
2281  Data  Data
2282  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
2283  {  {
2284    //    if (right.isLazy())
2285    // Construct the result using the given value and the other parameters    {
2286    // from right      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);
2287    Data result(left,right);      return Data(c);
2288    result+=right;    }
2289    return result;    return Data(left,right.getFunctionSpace(),false)+right;
2290  }  }
2291    
2292  //  //
# Line 1879  escript::operator+(const boost::python:: Line 2294  escript::operator+(const boost::python::
2294  Data  Data
2295  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
2296  {  {
2297    //    if (right.isLazy())
2298    // Construct the result using the given value and the other parameters    {
2299    // from right      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);
2300    Data result(left,right);      return Data(c);
2301    result-=right;    }
2302    return result;    return Data(left,right.getFunctionSpace(),false)-right;
2303  }  }
2304    
2305  //  //
# Line 1892  escript::operator-(const boost::python:: Line 2307  escript::operator-(const boost::python::
2307  Data  Data
2308  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
2309  {  {
2310    //    if (right.isLazy())
2311    // Construct the result using the given value and the other parameters    {
2312    // from right      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);
2313    Data result(left,right);      return Data(c);
2314    result*=right;    }
2315    return result;    return Data(left,right.getFunctionSpace(),false)*right;
2316  }  }
2317    
2318  //  //
# Line 1905  escript::operator*(const boost::python:: Line 2320  escript::operator*(const boost::python::
2320  Data  Data
2321  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
2322  {  {
2323    //    if (right.isLazy())
2324    // Construct the result using the given value and the other parameters    {
2325    // from right      DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);
2326    Data result(left,right);      return Data(c);
2327    result/=right;    }
2328    return result;    return Data(left,right.getFunctionSpace(),false)/right;
2329  }  }
2330    
 //  
 //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;  
 //}  
2331    
2332  /* TODO */  /* TODO */
2333  /* global reduction */  /* global reduction */
2334  Data  Data
2335  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
2336  {  {
   const DataArrayView& view=getPointDataView();  
2337    
2338    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2339    
2340    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=getDataPointRank()) {
2341      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2342    }    }
2343    
# Line 1977  Data::getItem(const boost::python::objec Line 2347  Data::getItem(const boost::python::objec
2347  /* TODO */  /* TODO */
2348  /* global reduction */  /* global reduction */
2349  Data  Data
2350  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataTypes::RegionType& region) const
2351  {  {
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
2352    return Data(*this,region);    return Data(*this,region);
2353  }  }
2354    
# Line 1995  Data::setItemO(const boost::python::obje Line 2362  Data::setItemO(const boost::python::obje
2362    setItemD(key,tempData);    setItemD(key,tempData);
2363  }  }
2364    
 /* TODO */  
 /* global reduction */  
2365  void  void
2366  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2367                 const Data& value)                 const Data& value)
2368  {  {
2369    const DataArrayView& view=getPointDataView();  //  const DataArrayView& view=getPointDataView();
2370    
2371    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2372    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=getDataPointRank()) {
2373      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
2374    }    }
2375    if (getFunctionSpace()!=value.getFunctionSpace()) {    if (getFunctionSpace()!=value.getFunctionSpace()) {
# Line 2014  Data::setItemD(const boost::python::obje Line 2379  Data::setItemD(const boost::python::obje
2379    }    }
2380  }  }
2381    
 /* TODO */  
 /* global reduction */  
2382  void  void
2383  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2384                 const DataArrayView::RegionType& region)                 const DataTypes::RegionType& region)
2385  {  {
2386    if (isProtected()) {    if (isProtected()) {
2387          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2388    }    }
2389  #if defined DOPROF    FORCERESOLVE;
2390    profData->slicing++;  /*  if (isLazy())
2391  #endif    {
2392        throw DataException("Error - setSlice not permitted on lazy data.");
2393      }*/
2394    Data tempValue(value);    Data tempValue(value);
2395    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2396    typeMatchRight(tempValue);    typeMatchRight(tempValue);
2397    m_data->setSlice(tempValue.m_data.get(),region);    getReady()->setSlice(tempValue.m_data.get(),region);
2398  }  }
2399    
2400  void  void
2401  Data::typeMatchLeft(Data& right) const  Data::typeMatchLeft(Data& right) const
2402  {  {
2403      if (right.isLazy() && !isLazy())
2404      {
2405        right.resolve();
2406      }
2407    if (isExpanded()){    if (isExpanded()){
2408      right.expand();      right.expand();
2409    } else if (isTagged()) {    } else if (isTagged()) {
# Line 2047  Data::typeMatchLeft(Data& right) const Line 2416  Data::typeMatchLeft(Data& right) const
2416  void  void
2417  Data::typeMatchRight(const Data& right)  Data::typeMatchRight(const Data& right)
2418  {  {
2419      if (isLazy() && !right.isLazy())
2420      {
2421        resolve();
2422      }
2423    if (isTagged()) {    if (isTagged()) {
2424      if (right.isExpanded()) {      if (right.isExpanded()) {
2425        expand();        expand();
# Line 2060  Data::typeMatchRight(const Data& right) Line 2433  Data::typeMatchRight(const Data& right)
2433    }    }
2434  }  }
2435    
2436  /* TODO */  void
2437  /* global reduction */  Data::setTaggedValueByName(std::string name,
2438                               const boost::python::object& value)
2439    {
2440         if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2441        FORCERESOLVE;
2442            int tagKey=getFunctionSpace().getDomain()->getTag(name);
2443            setTaggedValue(tagKey,value);
2444         }
2445    }
2446  void  void
2447  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2448                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2071  Data::setTaggedValue(int tagKey, Line 2452  Data::setTaggedValue(int tagKey,
2452    }    }
2453    //    //
2454    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2455    tag();    FORCERESOLVE;
2456      if (isConstant()) tag();
2457      numeric::array asNumArray(value);
2458    
2459    if (!isTagged()) {    // extract the shape of the numarray
2460      throw DataException("Error - DataTagged conversion failed!!");    DataTypes::ShapeType tempShape;
2461      for (int i=0; i < asNumArray.getrank(); i++) {
2462        tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2463    }    }
2464    
2465    //    DataVector temp_data2;
2466    // Construct DataArray from boost::python::object input value    temp_data2.copyFromNumArray(asNumArray);
   DataArray valueDataArray(value);  
2467    
2468    //    m_data->setTaggedValue(tagKey,tempShape, temp_data2);
   // Call DataAbstract::setTaggedValue  
   m_data->setTaggedValue(tagKey,valueDataArray.getView());  
2469  }  }
2470    
2471  /* TODO */  
 /* global reduction */  
2472  void  void
2473  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2474                              const DataArrayView& value)                  const DataTypes::ShapeType& pointshape,
2475                                const DataTypes::ValueType& value,
2476                    int dataOffset)
2477  {  {
2478    if (isProtected()) {    if (isProtected()) {
2479          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
2480    }    }
2481    //    //
2482    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2483    tag();    FORCERESOLVE;
2484      if (isConstant()) tag();
   if (!isTagged()) {  
     throw DataException("Error - DataTagged conversion failed!!");  
   }  
                                                                                                                 
2485    //    //
2486    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2487    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2488  }  }
2489    
 /* TODO */  
 /* global reduction */  
2490  int  int
2491  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2492  {  {
2493    return m_data->getTagNumber(dpno);    if (isEmpty())
2494  }    {
2495        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.");  
2496    }    }
2497    //    return getFunctionSpace().getTagFromDataPointNo(dpno);
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
2498  }  }
2499    
 /* 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);  
2500    
2501    //  ostream& escript::operator<<(ostream& o, const Data& data)
2502    // Load values from valueDataArray into return numarray  {
2503      o << data.toString();
2504    // extract the shape of the numarray    return o;
2505    int rank = value.getrank();  }
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
2506    
2507    if (rank==0) {  Data
2508        boost::python::numeric::array temp_numArray(valueView());  escript::C_GeneralTensorProduct(Data& arg_0,
2509        value = temp_numArray;                       Data& arg_1,
2510    }                       int axis_offset,
2511    if (rank==1) {                       int transpose)
2512      for (int i=0; i < shape[0]; i++) {  {
2513        value[i] = valueView(i);    // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2514      // SM is the product of the last axis_offset entries in arg_0.getShape().
2515    
2516      // deal with any lazy data
2517      if (arg_0.isLazy()) {arg_0.resolve();}
2518      if (arg_1.isLazy()) {arg_1.resolve();}
2519    
2520      // Interpolate if necessary and find an appropriate function space
2521      Data arg_0_Z, arg_1_Z;
2522      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2523        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2524          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2525          arg_1_Z = Data(arg_1);
2526        }
2527        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2528          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2529          arg_0_Z =Data(arg_0);
2530      }      }
2531    }      else {
2532    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);  
       }  
2533      }      }
2534      } else {
2535          arg_0_Z = Data(arg_0);
2536          arg_1_Z = Data(arg_1);
2537    }    }
2538    if (rank==3) {    // Get rank and shape of inputs
2539      for (int i=0; i < shape[0]; i++) {    int rank0 = arg_0_Z.getDataPointRank();
2540        for (int j=0; j < shape[1]; j++) {    int rank1 = arg_1_Z.getDataPointRank();
2541          for (int k=0; k < shape[2]; k++) {    const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2542            value[i][j][k] = valueView(i,j,k);    const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2543          }  
2544        }    // Prepare for the loops of the product and verify compatibility of shapes
2545      int start0=0, start1=0;
2546      if (transpose == 0)       {}
2547      else if (transpose == 1)  { start0 = axis_offset; }
2548      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2549      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2550    
2551    
2552      // Adjust the shapes for transpose
2553      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
2554      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
2555      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2556      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2557    
2558    #if 0
2559      // For debugging: show shape after transpose
2560      char tmp[100];
2561      std::string shapeStr;
2562      shapeStr = "(";
2563      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2564      shapeStr += ")";
2565      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2566      shapeStr = "(";
2567      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2568      shapeStr += ")";
2569      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2570    #endif
2571    
2572      // Prepare for the loops of the product
2573      int SL=1, SM=1, SR=1;
2574      for (int i=0; i<rank0-axis_offset; i++)   {
2575        SL *= tmpShape0[i];
2576      }
2577      for (int i=rank0-axis_offset; i<rank0; i++)   {
2578        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2579          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2580        }
2581        SM *= tmpShape0[i];
2582      }
2583      for (int i=axis_offset; i<rank1; i++)     {
2584        SR *= tmpShape1[i];
2585      }
2586    
2587      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2588      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
2589      {         // block to limit the scope of out_index
2590         int out_index=0;
2591         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2592         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2593      }
2594    
2595      // Declare output Data object
2596      Data res;
2597    
2598      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2599        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2600        double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2601        double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2602        double *ptr_2 = &(res.getDataAtOffset(0));
2603        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2604      }
2605      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2606    
2607        // Prepare the DataConstant input
2608        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2609        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2610    
2611        // Borrow DataTagged input from Data object
2612        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2613        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2614    
2615        // Prepare a DataTagged output 2
2616        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2617        res.tag();
2618        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2619        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2620    
2621        // Prepare offset into DataConstant
2622        int offset_0 = tmp_0->getPointOffset(0,0);
2623        double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2624        // Get the views
2625    //     DataArrayView view_1 = tmp_1->getDefaultValue();
2626    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2627    //     // Get the pointers to the actual data
2628    //     double *ptr_1 = &((view_1.getData())[0]);
2629    //     double *ptr_2 = &((view_2.getData())[0]);
2630    
2631        double *ptr_1 = &(tmp_1->getDefaultValue(0));
2632        double *ptr_2 = &(tmp_2->getDefaultValue(0));
2633    
2634    
2635        // Compute an MVP for the default
2636        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2637        // Compute an MVP for each tag
2638        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2639        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2640        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2641          tmp_2->addTag(i->first);
2642    //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2643    //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2644    //       double *ptr_1 = &view_1.getData(0);
2645    //       double *ptr_2 = &view_2.getData(0);
2646    
2647          double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2648          double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2649        
2650          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2651      }      }
2652    
2653    }    }
2654    if (rank==4) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2655      for (int i=0; i < shape[0]; i++) {  
2656        for (int j=0; j < shape[1]; j++) {      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2657          for (int k=0; k < shape[2]; k++) {      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2658            for (int l=0; l < shape[3]; l++) {      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2659              value[i][j][k][l] = valueView(i,j,k,l);      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2660            }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2661          }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2662        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2663        int sampleNo_1,dataPointNo_1;
2664        int numSamples_1 = arg_1_Z.getNumSamples();
2665        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2666        int offset_0 = tmp_0->getPointOffset(0,0);
2667        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2668        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2669          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2670            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2671            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2672            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2673            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2674            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2675            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2676        }        }
2677      }      }
2678    
2679    }    }
2680      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2681    
2682  }      // Borrow DataTagged input from Data object
2683        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2684        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2685    
2686  void      // Prepare the DataConstant input
2687  Data::archiveData(const std::string fileName)      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2688  {      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
   cout << "Archiving Data object to: " << fileName << endl;  
2689    
2690    //      // Prepare a DataTagged output 2
2691    // Determine type of this Data object      res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2692    int dataType = -1;      res.tag();
2693        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2694        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2695    
2696    if (isEmpty()) {      // Prepare offset into DataConstant
2697      dataType = 0;      int offset_1 = tmp_1->getPointOffset(0,0);
2698      cout << "\tdataType: DataEmpty" << endl;      double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2699    }      // Get the views
2700    if (isConstant()) {  //     DataArrayView view_0 = tmp_0->getDefaultValue();
2701      dataType = 1;  //     DataArrayView view_2 = tmp_2->getDefaultValue();
2702      cout << "\tdataType: DataConstant" << endl;  //     // Get the pointers to the actual data
2703    }  //     double *ptr_0 = &((view_0.getData())[0]);
2704    if (isTagged()) {  //     double *ptr_2 = &((view_2.getData())[0]);
     dataType = 2;  
     cout << "\tdataType: DataTagged" << endl;  
   }  
   if (isExpanded()) {  
     dataType = 3;  
     cout << "\tdataType: DataExpanded" << endl;  
   }  
2705    
2706    if (dataType == -1) {      double *ptr_0 = &(tmp_0->getDefaultValue(0));
2707      throw DataException("archiveData Error: undefined dataType");      double *ptr_2 = &(tmp_2->getDefaultValue(0));
   }  
2708    
2709    //      // Compute an MVP for the default
2710    // Collect data items common to all Data types      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2711    int noSamples = getNumSamples();      // Compute an MVP for each tag
2712    int noDPPSample = getNumDataPointsPerSample();      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2713    int functionSpaceType = getFunctionSpace().getTypeCode();      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2714    int dataPointRank = getDataPointRank();      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2715    int dataPointSize = getDataPointSize();  //      tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2716    int dataLength = getLength();  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2717    DataArrayView::ShapeType dataPointShape = getDataPointShape();  //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2718    vector<int> referenceNumbers(noSamples);  //       double *ptr_0 = &view_0.getData(0);
2719    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  //       double *ptr_2 = &view_2.getData(0);
2720      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);  
2721    }        tmp_2->addTag(i->first);
2722    vector<int> tagNumbers(noSamples);        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2723    if (isTagged()) {        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2724      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
       tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);  
2725      }      }
2726    
2727    }    }
2728      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2729    
2730    cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;      // Borrow DataTagged input from Data object
2731    cout << "\tfunctionSpaceType: " << functionSpaceType << endl;      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2732    cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2733    
2734    //      // Borrow DataTagged input from Data object
2735    // Flatten Shape to an array of integers suitable for writing to file      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2736    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;  
2737    
2738    //      // Prepare a DataTagged output 2
2739    // Open archive file      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2740    ofstream archiveFile;      res.tag();  // DataTagged output
2741    archiveFile.open(fileName.data(), ios::out);      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2742        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2743    
2744    if (!archiveFile.good()) {  //     // Get the views
2745      throw DataException("archiveData Error: problem opening archive file");  //     DataArrayView view_0 = tmp_0->getDefaultValue();
2746    }  //     DataArrayView view_1 = tmp_1->getDefaultValue();
2747    //     DataArrayView view_2 = tmp_2->getDefaultValue();
2748    //     // Get the pointers to the actual data
2749    //     double *ptr_0 = &((view_0.getData())[0]);
2750    //     double *ptr_1 = &((view_1.getData())[0]);
2751    //     double *ptr_2 = &((view_2.getData())[0]);
2752    
2753    //      double *ptr_0 = &(tmp_0->getDefaultValue(0));
2754    // Write common data items to archive file      double *ptr_1 = &(tmp_1->getDefaultValue(0));
2755    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));      double *ptr_2 = &(tmp_2->getDefaultValue(0));
   archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));  
   archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));  
   for (int dim = 0; dim < 4; dim++) {  
     archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));  
   }  
   for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
     archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));  
   }  
   if (isTagged()) {  
     for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
       archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));  
     }  
   }  
2756    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem writing to archive file");  
   }  
2757    
2758    //      // Compute an MVP for the default
2759    // Archive underlying data values for each Data type      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2760    int noValues;      // Merge the tags
2761    switch (dataType) {      DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2762      case 0:      const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2763        // DataEmpty      const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2764        noValues = 0;      for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2765        archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));        tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2766        cout << "\tnoValues: " << noValues << endl;      }
2767        break;      for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2768      case 1:        tmp_2->addTag(i->first);
2769        // DataConstant      }
2770        noValues = m_data->getLength();      // Compute an MVP for each tag
2771        archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));      const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2772        cout << "\tnoValues: " << noValues << endl;      for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2773        if (m_data->archiveData(archiveFile,noValues)) {  //       DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2774          throw DataException("archiveData Error: problem writing data to archive file");  //       DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2775        }  //       DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2776        break;  //       double *ptr_0 = &view_0.getData(0);
2777      case 2:  //       double *ptr_1 = &view_1.getData(0);
2778        // DataTagged  //       double *ptr_2 = &view_2.getData(0);
       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");  
       }  
       break;  
     case 3:  
       // DataExpanded  
       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");  
       }  
       break;  
   }  
2779    
2780    if (!archiveFile.good()) {        double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2781      throw DataException("archiveData Error: problem writing data to archive file");        double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2782    }        double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2783    
2784    //        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2785    // Close archive file      }
   archiveFile.close();  
2786    
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
2787    }    }
2788      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2789    
2790  }      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2791        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2792        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2793        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2794        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2795        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2796        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2797        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2798        int sampleNo_0,dataPointNo_0;
2799        int numSamples_0 = arg_0_Z.getNumSamples();
2800        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2801        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2802        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2803          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2804          double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2805          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_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            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2809            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2810            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2811          }
2812        }
2813    
 void  
 Data::extractData(const std::string fileName,  
                   const FunctionSpace& fspace)  
 {  
   //  
   // Can only extract Data to an object which is initially DataEmpty  
   if (!isEmpty()) {  
     throw DataException("extractData Error: can only extract to DataEmpty object");  
2814    }    }
2815      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2816    
2817    cout << "Extracting Data object from: " << fileName << endl;      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2818        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2819    int dataType;      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2820    int noSamples;      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2821    int noDPPSample;      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2822    int functionSpaceType;      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2823    int dataPointRank;      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2824    int dataPointSize;      int sampleNo_0,dataPointNo_0;
2825    int dataLength;      int numSamples_0 = arg_0_Z.getNumSamples();
2826    DataArrayView::ShapeType dataPointShape;      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2827    int flatShape[4];      int offset_1 = tmp_1->getPointOffset(0,0);
2828        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2829        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2830          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2831            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2832            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2833            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2834            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2835            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2836            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2837          }
2838        }
2839    
   //  
   // Open the archive file  
   ifstream archiveFile;  
   archiveFile.open(fileName.data(), ios::in);  
2840    
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem opening archive file");  
2841    }    }
2842      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2843    
2844    //      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2845    // Read common data items from archive file      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2846    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2847    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2848    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2849    archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2850    archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2851    archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2852    archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));      int sampleNo_0,dataPointNo_0;
2853    for (int dim = 0; dim < 4; dim++) {      int numSamples_0 = arg_0_Z.getNumSamples();
2854      archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2855      if (flatShape[dim]>0) {      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2856        dataPointShape.push_back(flatShape[dim]);      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2857          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2858          double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2859          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2860            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2861            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2862            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2863            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2864            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2865          }
2866      }      }
2867    
2868    }    }
2869    vector<int> referenceNumbers(noSamples);    else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2870    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  
2871      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2872    }      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2873    vector<int> tagNumbers(noSamples);      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2874    if (dataType==2) {      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2875      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2876        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2877        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2878        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2879        int sampleNo_0,dataPointNo_0;
2880        int numSamples_0 = arg_0_Z.getNumSamples();
2881        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2882        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2883        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2884          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2885            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2886            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2887            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2888            double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2889            double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2890            double *ptr_2 = &(res.getDataAtOffset(offset_2));
2891            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2892          }
2893      }      }
   }  
2894    
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem reading from archive file");  
2895    }    }
2896      else {
2897    //      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
   // Verify the values just read from the archive file  
   switch (dataType) {  
     case 0:  
       cout << "\tdataType: DataEmpty" << endl;  
       break;  
     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;  
2898    }    }
2899    
2900    cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;    return res;
2901    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;  
2902    
2903    //  DataAbstract*
2904    // Verify that supplied FunctionSpace object is compatible with this Data object.  Data::borrowData() const
2905    if ( (fspace.getTypeCode()!=functionSpaceType) ||  {
2906         (fspace.getNumSamples()!=noSamples) ||    return m_data.get();
2907         (fspace.getNumDPPSample()!=noDPPSample)  }
2908       ) {  
2909      throw DataException("extractData Error: incompatible FunctionSpace");  // Not all that happy about returning a non-const from a const
2910    }  DataAbstract_ptr
2911    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  Data::borrowDataPtr() const
2912      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {  {
2913        throw DataException("extractData Error: incompatible FunctionSpace");    return m_data;
2914      }  }
2915    }  
2916    if (dataType==2) {  // Not all that happy about returning a non-const from a const
2917      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {  DataReady_ptr
2918        if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {  Data::borrowReadyPtr() const
2919          throw DataException("extractData Error: incompatible FunctionSpace");  {
2920        }     DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
2921       EsysAssert((dr!=0), "Error - casting to DataReady.");
2922       return dr;
2923    }
2924    
2925    std::string
2926    Data::toString() const
2927    {
2928        if (!m_data->isEmpty() &&
2929        getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))
2930        {
2931        stringstream temp;
2932        temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
2933        return  temp.str();
2934      }      }
2935    }      return m_data->toString();
2936    }
2937    
   //  
   // Construct a DataVector to hold underlying data values  
   DataVector dataVec(dataLength);  
2938    
   //  
   // Load this DataVector with the appropriate values  
   int noValues;  
   archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));  
   cout << "\tnoValues: " << noValues << endl;  
   switch (dataType) {  
     case 0:  
       // DataEmpty  
       if (noValues != 0) {  
         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;  
   }  
2939    
2940    if (!archiveFile.good()) {  DataTypes::ValueType::const_reference
2941      throw DataException("extractData Error: problem reading from archive file");  Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
2942    }  {
2943        if (isLazy())
2944        {
2945        throw DataException("Programmer error - getDataAtOffset not permitted on lazy data (object is const which prevents resolving).");
2946        }
2947        return getReady()->getDataAtOffset(i);
2948    }
2949    
   //  
   // Close archive file  
   archiveFile.close();  
2950    
2951    if (!archiveFile.good()) {  DataTypes::ValueType::reference
2952      throw DataException("extractData Error: problem closing archive file");  Data::getDataAtOffset(DataTypes::ValueType::size_type i)
2953    }  {
2954    //     if (isLazy())
2955    //     {
2956    //  throw DataException("getDataAtOffset not permitted on lazy data.");
2957    //     }
2958        FORCERESOLVE;
2959        return getReady()->getDataAtOffset(i);
2960    }
2961    
2962    //  DataTypes::ValueType::const_reference
2963    // Construct an appropriate Data object  Data::getDataPoint(int sampleNo, int dataPointNo) const
2964    DataAbstract* tempData;  {
2965    switch (dataType) {    if (!isReady())
2966      case 0:    {
2967        // DataEmpty      throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data (object is const which prevents resolving).");
2968        tempData=new DataEmpty();    }
2969        break;    else
2970      case 1:    {
2971        // DataConstant      const DataReady* dr=getReady();
2972        tempData=new DataConstant(fspace,dataPointShape,dataVec);      return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
       break;  
     case 2:  
       // DataTagged  
       tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);  
       break;  
     case 3:  
       // DataExpanded  
       tempData=new DataExpanded(fspace,dataPointShape,dataVec);  
       break;  
2973    }    }
   shared_ptr<DataAbstract> temp_data(tempData);  
   m_data=temp_data;  
2974  }  }
2975    
2976  ostream& escript::operator<<(ostream& o, const Data& data)  
2977    DataTypes::ValueType::reference
2978    Data::getDataPoint(int sampleNo, int dataPointNo)
2979  {  {
2980    o << data.toString();    FORCERESOLVE;
2981    return o;    if (!isReady())
2982      {
2983        throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");
2984      }
2985      else
2986      {
2987        DataReady* dr=getReady();
2988        return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
2989      }
2990  }  }
2991    
2992    
2993  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2994    
2995  void  void
2996  Data::print()  Data::print()
2997  {  {
2998    int i,j;    int i,j;
2999      
3000    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
3001    for( i=0; i<getNumSamples(); i++ )    for( i=0; i<getNumSamples(); i++ )
3002    {    {
# Line 2557  Data::print() Line 3006  Data::print()
3006      printf( "\n" );      printf( "\n" );
3007    }    }
3008  }  }
3009    void
3010    Data::dump(const std::string fileName) const
3011    {
3012      try
3013      {
3014        if (isLazy())
3015        {
3016          Data temp(*this); // this is to get a non-const object which we can resolve
3017          temp.resolve();
3018          temp.dump(fileName);
3019        }
3020        else
3021        {
3022              return m_data->dump(fileName);
3023        }
3024      }
3025      catch (exception& e)
3026      {
3027            cout << e.what() << endl;
3028      }
3029    }
3030    
3031  int  int
3032  Data::get_MPISize() const  Data::get_MPISize() const
3033  {  {
3034      int error, size;      int size;
3035  #ifdef PASO_MPI  #ifdef PASO_MPI
3036        int error;
3037      error = MPI_Comm_size( get_MPIComm(), &size );      error = MPI_Comm_size( get_MPIComm(), &size );
3038  #else  #else
3039      size = 1;      size = 1;
# Line 2573  Data::get_MPISize() const Line 3044  Data::get_MPISize() const
3044  int  int
3045  Data::get_MPIRank() const  Data::get_MPIRank() const
3046  {  {
3047      int error, rank;      int rank;
3048  #ifdef PASO_MPI  #ifdef PASO_MPI
3049        int error;
3050      error = MPI_Comm_rank( get_MPIComm(), &rank );      error = MPI_Comm_rank( get_MPIComm(), &rank );
3051  #else  #else
3052      rank = 0;      rank = 0;
# Line 2584  Data::get_MPIRank() const Line 3056  Data::get_MPIRank() const
3056    
3057  MPI_Comm  MPI_Comm
3058  Data::get_MPIComm() const  Data::get_MPIComm() const
3059  {  {
3060  #ifdef PASO_MPI  #ifdef PASO_MPI
3061      return MPI_COMM_WORLD;      return MPI_COMM_WORLD;
3062  #else  #else
# Line 2592  Data::get_MPIComm() const Line 3064  Data::get_MPIComm() const
3064  #endif  #endif
3065  }  }
3066    
3067    

Legend:
Removed from v.790  
changed lines
  Added in v.2048

  ViewVC Help
Powered by ViewVC 1.1.26