/[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 1319 by matt, Thu Sep 27 00:27:51 2007 UTC
# Line 1  Line 1 
 // $Id$  
1    
2  /*  /* $Id$ */
3   ************************************************************  
4   *          Copyright 2006 by ACcESS MNRF                   *  /*******************************************************
5   *                                                          *   *
6   *              http://www.access.edu.au                    *   *           Copyright 2003-2007 by ACceSS MNRF
7   *       Primary Business: Queensland, Australia            *   *       Copyright 2007 by University of Queensland
8   *  Licensed under the Open Software License version 3.0    *   *
9   *     http://www.opensource.org/licenses/osl-3.0.php       *   *                http://esscc.uq.edu.au
10   *                                                          *   *        Primary Business: Queensland, Australia
11   ************************************************************   *  Licensed under the Open Software License version 3.0
12  */   *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  #include "Data.h"  #include "Data.h"
17    
18  #include "DataExpanded.h"  #include "DataExpanded.h"
19  #include "DataConstant.h"  #include "DataConstant.h"
20  #include "DataTagged.h"  #include "DataTagged.h"
21  #include "DataEmpty.h"  #include "DataEmpty.h"
 #include "DataArray.h"  
22  #include "DataArrayView.h"  #include "DataArrayView.h"
 #include "DataProf.h"  
23  #include "FunctionSpaceFactory.h"  #include "FunctionSpaceFactory.h"
24  #include "AbstractContinuousDomain.h"  #include "AbstractContinuousDomain.h"
25  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
26    extern "C" {
27    #include "escript/blocktimer.h"
28    }
29    
30  #include <fstream>  #include <fstream>
31  #include <algorithm>  #include <algorithm>
# Line 38  using namespace boost::python; Line 41  using namespace boost::python;
41  using namespace boost;  using namespace boost;
42  using namespace escript;  using namespace escript;
43    
 #if defined DOPROF  
 //  
 // global table of profiling data for all Data objects  
 DataProf dataProfTable;  
 #endif  
   
44  Data::Data()  Data::Data()
45  {  {
46    //    //
# Line 52  Data::Data() Line 49  Data::Data()
49    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
50    m_data=temp_data;    m_data=temp_data;
51    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
52  }  }
53    
54  Data::Data(double value,  Data::Data(double value,
# Line 67  Data::Data(double value, Line 60  Data::Data(double value,
60    for (int i = 0; i < shape.attr("__len__")(); ++i) {    for (int i = 0; i < shape.attr("__len__")(); ++i) {
61      dataPointShape.push_back(extract<const int>(shape[i]));      dataPointShape.push_back(extract<const int>(shape[i]));
62    }    }
63    DataArray temp(dataPointShape,value);  
64    initialise(temp.getView(),what,expanded);    int len = DataArrayView::noValues(dataPointShape);
65      DataVector temp_data(len,value,len);
66      DataArrayView temp_dataView(temp_data, dataPointShape);
67    
68      initialise(temp_dataView, what, expanded);
69    
70    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
71  }  }
72    
73  Data::Data(double value,  Data::Data(double value,
# Line 81  Data::Data(double value, Line 75  Data::Data(double value,
75         const FunctionSpace& what,         const FunctionSpace& what,
76             bool expanded)             bool expanded)
77  {  {
78    DataArray temp(dataPointShape,value);    int len = DataArrayView::noValues(dataPointShape);
79    pair<int,int> dataShape=what.getDataShape();  
80    initialise(temp.getView(),what,expanded);    DataVector temp_data(len,value,len);
81      DataArrayView temp_dataView(temp_data, dataPointShape);
82    
83      initialise(temp_dataView, what, expanded);
84    
85    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
86  }  }
87    
88  Data::Data(const Data& inData)  Data::Data(const Data& inData)
89  {  {
90    m_data=inData.m_data;    m_data=inData.m_data;
91    m_protected=inData.isProtected();    m_protected=inData.isProtected();
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
92  }  }
93    
94  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 110  Data::Data(const Data& inData, Line 100  Data::Data(const Data& inData,
100    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
101    m_data=temp_data;    m_data=temp_data;
102    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
103  }  }
104    
105  Data::Data(const Data& inData,  Data::Data(const Data& inData,
106             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
107  {  {
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
108    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
109      m_data=inData.m_data;      m_data=inData.m_data;
110    } else {    } else {
     #if defined DOPROF  
     profData->interpolate++;  
     #endif  
111      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
112      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
113      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 158  Data::Data(const DataTagged::TagListType Line 137  Data::Data(const DataTagged::TagListType
137    if (expanded) {    if (expanded) {
138      expand();      expand();
139    }    }
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
140  }  }
141    
142  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 170  Data::Data(const numeric::array& value, Line 145  Data::Data(const numeric::array& value,
145  {  {
146    initialise(value,what,expanded);    initialise(value,what,expanded);
147    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
148  }  }
149    
150  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 182  Data::Data(const DataArrayView& value, Line 153  Data::Data(const DataArrayView& value,
153  {  {
154    initialise(value,what,expanded);    initialise(value,what,expanded);
155    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
156  }  }
157    
158  Data::Data(const object& value,  Data::Data(const object& value,
# Line 195  Data::Data(const object& value, Line 162  Data::Data(const object& value,
162    numeric::array asNumArray(value);    numeric::array asNumArray(value);
163    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
164    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
165  }  }
166    
167    
168  Data::Data(const object& value,  Data::Data(const object& value,
169             const Data& other)             const Data& other)
170  {  {
171    
172      numeric::array asNumArray(value);
173    
174    
175      // extract the shape of the numarray
176      DataArrayView::ShapeType tempShape;
177      for (int i=0; i < asNumArray.getrank(); i++) {
178        tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
179      }
180      // get the space for the data vector
181      int len = DataArrayView::noValues(tempShape);
182      DataVector temp_data(len, 0.0, len);
183      DataArrayView temp_dataView(temp_data, tempShape);
184      temp_dataView.copy(asNumArray);
185    
186    //    //
187    // Create DataConstant using the given value and all other parameters    // Create DataConstant using the given value and all other parameters
188    // copied from other. If value is a rank 0 object this Data    // copied from other. If value is a rank 0 object this Data
189    // will assume the point data shape of other.    // will assume the point data shape of other.
190    DataArray temp(value);  
191    if (temp.getView().getRank()==0) {    if (temp_dataView.getRank()==0) {
192      //      int len = DataArrayView::noValues(other.getPointDataView().getShape());
193      // Create a DataArray with the scalar value for all elements  
194      DataArray temp2(other.getPointDataView().getShape(),temp.getView()());      DataVector temp2_data(len, temp_dataView(), len);
195      initialise(temp2.getView(),other.getFunctionSpace(),false);      DataArrayView temp2_dataView(temp_data, other.getPointDataView().getShape());
196        initialise(temp2_dataView, other.getFunctionSpace(), false);
197    
198    } else {    } else {
199      //      //
200      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
201      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp_dataView, other.getFunctionSpace(), false);
202    }    }
203    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
204  }  }
205    
206  Data::~Data()  Data::~Data()
# Line 266  Data::getShapeTuple() const Line 243  Data::getShapeTuple() const
243          throw DataException("Error - illegal Data rank.");          throw DataException("Error - illegal Data rank.");
244    }    }
245  }  }
   
246  void  void
247  Data::copy(const Data& other)  Data::copy(const Data& other)
248  {  {
# Line 319  Data::copy(const Data& other) Line 295  Data::copy(const Data& other)
295    throw DataException("Error - Copy not implemented for this Data type.");    throw DataException("Error - Copy not implemented for this Data type.");
296  }  }
297    
298    
299    void
300    Data::setToZero()
301    {
302      {
303        DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
304        if (temp!=0) {
305           temp->setToZero();
306           return;
307        }
308      }
309      {
310        DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
311        if (temp!=0) {
312          temp->setToZero();
313          return;
314        }
315      }
316      {
317        DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
318        if (temp!=0) {
319          temp->setToZero();
320          return;
321        }
322      }
323      throw DataException("Error - Data can not be set to zero.");
324    }
325    
326  void  void
327  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
328                     const Data& mask)                     const Data& mask)
# Line 350  Data::isTagged() const Line 354  Data::isTagged() const
354    return (temp!=0);    return (temp!=0);
355  }  }
356    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
357  bool  bool
358  Data::isEmpty() const  Data::isEmpty() const
359  {  {
# Line 367  Data::isConstant() const Line 369  Data::isConstant() const
369  }  }
370    
371  void  void
372  Data::setProtection()  Data::setProtection()
373  {  {
374     m_protected=true;     m_protected=true;
375  }  }
376    
377  bool  bool
378  Data::isProtected() const  Data::isProtected() const
379  {  {
380     return m_protected;     return m_protected;
381  }  }
382    
# Line 422  Data::tag() Line 424  Data::tag()
424    }    }
425  }  }
426    
427  void  Data
428  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
429  {  {
430    m_data->reshapeDataPoint(shape);    return escript::unaryOp(*this,bind1st(divides<double>(),1.));
431  }  }
432    
433  Data  Data
434  Data::wherePositive() const  Data::wherePositive() const
435  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
436    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
437  }  }
438    
439  Data  Data
440  Data::whereNegative() const  Data::whereNegative() const
441  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
442    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
443  }  }
444    
445  Data  Data
446  Data::whereNonNegative() const  Data::whereNonNegative() const
447  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
448    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
449  }  }
450    
451  Data  Data
452  Data::whereNonPositive() const  Data::whereNonPositive() const
453  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
454    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
455  }  }
456    
457  Data  Data
458  Data::whereZero(double tol) const  Data::whereZero(double tol) const
459  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
460    Data dataAbs=abs();    Data dataAbs=abs();
461    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
462  }  }
# Line 477  Data::whereZero(double tol) const Line 464  Data::whereZero(double tol) const
464  Data  Data
465  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
466  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
467    Data dataAbs=abs();    Data dataAbs=abs();
468    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
469  }  }
# Line 487  Data::whereNonZero(double tol) const Line 471  Data::whereNonZero(double tol) const
471  Data  Data
472  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
473  {  {
 #if defined DOPROF  
   profData->interpolate++;  
 #endif  
474    return Data(*this,functionspace);    return Data(*this,functionspace);
475  }  }
476    
# Line 511  Data::probeInterpolation(const FunctionS Line 492  Data::probeInterpolation(const FunctionS
492  Data  Data
493  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
494  {  {
495  #if defined DOPROF    double blocktimer_start = blocktimer_time();
   profData->grad++;  
 #endif  
496    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
497      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
498    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
499    grad_shape.push_back(functionspace.getDim());    grad_shape.push_back(functionspace.getDim());
500    Data out(0.0,grad_shape,functionspace,true);    Data out(0.0,grad_shape,functionspace,true);
501    getDomain().setToGradient(out,*this);    getDomain().setToGradient(out,*this);
502      blocktimer_increment("grad()", blocktimer_start);
503    return out;    return out;
504  }  }
505    
# Line 547  Data::getDataPointShape() const Line 527  Data::getDataPointShape() const
527    return getPointDataView().getShape();    return getPointDataView().getShape();
528  }  }
529    
 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");  
   }  
530    
   //  
   // make sure data is expanded:  
   if (!isExpanded()) {  
     expand();  
   }  
   
   //  
   // and copy over  
   m_data->copyAll(num_array);  
 }  
531    
532  const  const
533  boost::python::numeric::array  boost::python::numeric::array
534  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
535  {  {
536    //    size_t length=0;
537    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
538    //    //
539    // determine the rank and shape of each data point    // determine the rank and shape of each data point
540    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 596  Data::convertToNumArray() Line 545  Data::convertToNumArray()
545    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
546    
547    //    //
548    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
549    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
550    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
551    // data points, whilst the first dimension will be equal to the    DataArrayView::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]);  
   }  
552    
553    //    //
554    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
555      if (arrayRank==0) {
556        numArray.resize(1);
557      }
558    if (arrayRank==1) {    if (arrayRank==1) {
559      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
560    }    }
# Line 623  Data::convertToNumArray() Line 567  Data::convertToNumArray()
567    if (arrayRank==4) {    if (arrayRank==4) {
568      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
569    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
570    
571    //    if (getNumDataPointsPerSample()>0) {
572    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
573    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
574    int dataPoint = 0;         //
575    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
576      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
577        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
578        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++;  
     }  
   }  
579    
580           //
581           // Check a valid data point number has been supplied
582           if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
583               throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
584           }
585           // TODO: global error handling
586           // create a view of the data if it is stored locally
587           DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
588    
589           switch( dataPointRank ){
590                case 0 :
591                    numArray[0] = dataPointView();
592                    break;
593                case 1 :
594                    for( i=0; i<dataPointShape[0]; i++ )
595                        numArray[i]=dataPointView(i);
596                    break;
597                case 2 :
598                    for( i=0; i<dataPointShape[0]; i++ )
599                        for( j=0; j<dataPointShape[1]; j++)
600                            numArray[make_tuple(i,j)]=dataPointView(i,j);
601                    break;
602                case 3 :
603                    for( i=0; i<dataPointShape[0]; i++ )
604                        for( j=0; j<dataPointShape[1]; j++ )
605                            for( k=0; k<dataPointShape[2]; k++)
606                                numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
607                    break;
608                case 4 :
609                    for( i=0; i<dataPointShape[0]; i++ )
610                        for( j=0; j<dataPointShape[1]; j++ )
611                            for( k=0; k<dataPointShape[2]; k++ )
612                                for( l=0; l<dataPointShape[3]; l++)
613                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
614                    break;
615        }
616      }
617    //    //
618    // return the loaded array    // return the array
619    return numArray;    return numArray;
 }  
620    
621  const  }
622  boost::python::numeric::array  void
623  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
624  {  {
625    //      // this will throw if the value cannot be represented
626    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
627    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
628    
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
629    
630    //  }
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
631    
632    void
633    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
634    {
635      if (isProtected()) {
636            throw DataException("Error - attempt to update protected Data object.");
637      }
638    //    //
639    // create the numeric array to be returned    // check rank
640    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
641          throw DataException("Rank of numarray does not match Data object rank");
642    
643    //    //
644    // the rank of the returned numeric array will be the rank of    // check shape of num_array
645    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
646    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
647    // 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]);  
648    }    }
   
649    //    //
650    // resize the numeric array to the shape just calculated    // make sure data is expanded:
651    if (arrayRank==1) {    if (!isExpanded()) {
652      numArray.resize(arrayShape[0]);      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]);  
653    }    }
654    if (arrayRank==5) {    if (getNumDataPointsPerSample()>0) {
655      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);         int sampleNo = dataPointNo/getNumDataPointsPerSample();
656           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
657           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
658      } else {
659           m_data->copyToDataPoint(-1, 0,num_array);
660    }    }
661    }
662    
663    //  void
664    // loop through each data point in turn, loading the values for that data point  Data::setValueOfDataPoint(int dataPointNo, const double value)
665    // into the numeric array.  {
666    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    if (isProtected()) {
667      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);  
             }  
           }  
         }  
       }  
     }  
668    }    }
   
669    //    //
670    // return the loaded array    // make sure data is expanded:
671    return numArray;    if (!isExpanded()) {
672        expand();
673      }
674      if (getNumDataPointsPerSample()>0) {
675           int sampleNo = dataPointNo/getNumDataPointsPerSample();
676           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
677           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
678      } else {
679           m_data->copyToDataPoint(-1, 0,value);
680      }
681  }  }
682    
683  const  const
684  boost::python::numeric::array  boost::python::numeric::array
685  Data::convertToNumArrayFromDPNo(int procNo,  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
                                 int sampleNo,  
                                                                 int dataPointNo)  
   
686  {  {
687      size_t length=0;    size_t length=0;
688      int i, j, k, l, pos;    int i, j, k, l, pos;
   
   //  
   // 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.");  
   }  
   
689    //    //
690    // determine the rank and shape of each data point    // determine the rank and shape of each data point
691    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 719  Data::convertToNumArrayFromDPNo(int proc
719      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
720    }    }
721    
722      // added for the MPI communication    // added for the MPI communication
723      length=1;    length=1;
724      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
725          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
726    
727    //    //
728    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
729    
730      // updated for the MPI case      // updated for the MPI case
731      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
732                 if (getNumDataPointsPerSample()>0) {
733                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
734                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
735                    //
736                    // Check a valid sample number has been supplied
737                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
738                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
739                    }
740    
741                    //
742                    // Check a valid data point number has been supplied
743                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
744                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
745                    }
746                    // TODO: global error handling
747          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
748          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
749            
750          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
751          pos=0;          pos=0;
752          switch( dataPointRank ){          switch( dataPointRank ){
753              case 0 :              case 0 :
754                  tmpData[0] = dataPointView();                  tmpData[0] = dataPointView();
755                  break;                  break;
756              case 1 :                      case 1 :
757                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
758                      tmpData[i]=dataPointView(i);                      tmpData[i]=dataPointView(i);
759                  break;                  break;
760              case 2 :                      case 2 :
761                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
762                      for( j=0; j<dataPointShape[1]; j++, pos++ )                      for( j=0; j<dataPointShape[1]; j++, pos++ )
763                          tmpData[pos]=dataPointView(i,j);                          tmpData[pos]=dataPointView(i,j);
764                  break;                  break;
765              case 3 :                      case 3 :
766                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
767                      for( j=0; j<dataPointShape[1]; j++ )                      for( j=0; j<dataPointShape[1]; j++ )
768                          for( k=0; k<dataPointShape[2]; k++, pos++ )                          for( k=0; k<dataPointShape[2]; k++, pos++ )
# Line 878  Data::convertToNumArrayFromDPNo(int proc Line 776  Data::convertToNumArrayFromDPNo(int proc
776                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=dataPointView(i,j,k,l);
777                  break;                  break;
778          }          }
779                }
780      }      }
781  #ifdef PASO_MPI          #ifdef PASO_MPI
782          // broadcast the data to all other processes          // broadcast the data to all other processes
783          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
784  #endif          #endif
785    
786      // unpack the data      // unpack the data
787      switch( dataPointRank ){      switch( dataPointRank ){
788          case 0 :          case 0 :
789              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
790              break;              break;
791          case 1 :                  case 1 :
792              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
793                  numArray[i]=tmpData[i];                  numArray[i]=tmpData[i];
794              break;              break;
795          case 2 :                  case 2 :
796              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
797                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
798                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
799              break;              break;
800          case 3 :                  case 3 :
801              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
802                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
803                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
804                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
805              break;              break;
806          case 4 :          case 4 :
807              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
808                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
809                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
810                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
811                              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]))];
812              break;              break;
813      }      }
814    
815      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);  
           }  
         }  
       }  
     }  
   }  
 */  
   
816    //    //
817    // return the loaded array    // return the loaded array
818    return numArray;    return numArray;
819  }  }
820    
821    
822    
823  boost::python::numeric::array  boost::python::numeric::array
824  Data::integrate() const  Data::integrate() const
825  {  {
826    int index;    int index;
827    int rank = getDataPointRank();    int rank = getDataPointRank();
828    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
829      int dataPointSize = getDataPointSize();
 #if defined DOPROF  
   profData->integrate++;  
 #endif  
830    
831    //    //
832    // calculate the integral values    // calculate the integral values
833    vector<double> integrals(getDataPointSize());    vector<double> integrals(dataPointSize);
834      vector<double> integrals_local(dataPointSize);
835    #ifdef PASO_MPI
836      AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
837      // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
838      double *tmp = new double[dataPointSize];
839      double *tmp_local = new double[dataPointSize];
840      for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
841      MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
842      for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
843      delete[] tmp;
844      delete[] tmp_local;
845    #else
846    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
847    #endif
848    
849    //    //
850    // create the numeric array to be returned    // create the numeric array to be returned
# Line 1031  Data::integrate() const Line 904  Data::integrate() const
904  Data  Data
905  Data::sin() const  Data::sin() const
906  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
907    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
908  }  }
909    
910  Data  Data
911  Data::cos() const  Data::cos() const
912  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
913    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
914  }  }
915    
916  Data  Data
917  Data::tan() const  Data::tan() const
918  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
919    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
920  }  }
921    
922  Data  Data
923  Data::asin() const  Data::asin() const
924  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
925    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
926  }  }
927    
928  Data  Data
929  Data::acos() const  Data::acos() const
930  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
931    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
932  }  }
933    
934    
935  Data  Data
936  Data::atan() const  Data::atan() const
937  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
938    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
939  }  }
940    
941  Data  Data
942  Data::sinh() const  Data::sinh() const
943  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
944    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
945  }  }
946    
947  Data  Data
948  Data::cosh() const  Data::cosh() const
949  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
950    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
951  }  }
952    
953  Data  Data
954  Data::tanh() const  Data::tanh() const
955  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
956    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
957  }  }
958    
959    
960  Data  Data
961  Data::asinh() const  Data::erf() const
962  {  {
963  #if defined DOPROF  #ifdef _WIN32
964    profData->unary++;    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
965    #else
966      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
967  #endif  #endif
968    }
969    
970    Data
971    Data::asinh() const
972    {
973    #ifdef _WIN32
974      return escript::unaryOp(*this,escript::asinh_substitute);
975    #else
976    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
977    #endif
978  }  }
979    
980  Data  Data
981  Data::acosh() const  Data::acosh() const
982  {  {
983  #if defined DOPROF  #ifdef _WIN32
984    profData->unary++;    return escript::unaryOp(*this,escript::acosh_substitute);
985  #endif  #else
986    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
987    #endif
988  }  }
989    
990  Data  Data
991  Data::atanh() const  Data::atanh() const
992  {  {
993  #if defined DOPROF  #ifdef _WIN32
994    profData->unary++;    return escript::unaryOp(*this,escript::atanh_substitute);
995  #endif  #else
996    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
997    #endif
998  }  }
999    
1000  Data  Data
1001  Data::log10() const  Data::log10() const
1002  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1003    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1004  }  }
1005    
1006  Data  Data
1007  Data::log() const  Data::log() const
1008  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1009    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1010  }  }
1011    
1012  Data  Data
1013  Data::sign() const  Data::sign() const
1014  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1015    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1016  }  }
1017    
1018  Data  Data
1019  Data::abs() const  Data::abs() const
1020  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1021    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1022  }  }
1023    
1024  Data  Data
1025  Data::neg() const  Data::neg() const
1026  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1027    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1028  }  }
1029    
1030  Data  Data
1031  Data::pos() const  Data::pos() const
1032  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1033    Data result;    Data result;
1034    // perform a deep copy    // perform a deep copy
1035    result.copy(*this);    result.copy(*this);
# Line 1196  Data::pos() const Line 1039  Data::pos() const
1039  Data  Data
1040  Data::exp() const  Data::exp() const
1041  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1042    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1043  }  }
1044    
1045  Data  Data
1046  Data::sqrt() const  Data::sqrt() const
1047  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1048    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1049  }  }
1050    
# Line 1215  double Line 1052  double
1052  Data::Lsup() const  Data::Lsup() const
1053  {  {
1054    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1055    //    //
1056    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1057    
# Line 1232  Data::Lsup() const Line 1066  Data::Lsup() const
1066  }  }
1067    
1068  double  double
 Data::Linf() const  
 {  
   double localValue, globalValue;  
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
   //  
   // set the initial absolute minimum value to max double  
   AbsMin abs_min_func;  
   localValue = algorithm(abs_min_func,numeric_limits<double>::max());  
   
 #ifdef PASO_MPI  
   MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );  
   return globalValue;  
 #else  
   return localValue;  
 #endif  
 }  
   
 double  
1069  Data::sup() const  Data::sup() const
1070  {  {
1071    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1072    //    //
1073    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1074    FMax fmax_func;    FMax fmax_func;
# Line 1274  double Line 1085  double
1085  Data::inf() const  Data::inf() const
1086  {  {
1087    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1088    //    //
1089    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1090    FMin fmin_func;    FMin fmin_func;
# Line 1294  Data::inf() const Line 1102  Data::inf() const
1102  Data  Data
1103  Data::maxval() const  Data::maxval() const
1104  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1105    //    //
1106    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1107    FMax fmax_func;    FMax fmax_func;
# Line 1306  Data::maxval() const Line 1111  Data::maxval() const
1111  Data  Data
1112  Data::minval() const  Data::minval() const
1113  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1114    //    //
1115    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1116    FMin fmin_func;    FMin fmin_func;
# Line 1316  Data::minval() const Line 1118  Data::minval() const
1118  }  }
1119    
1120  Data  Data
1121  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1122  {  {
1123  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1124    profData->reduction2++;       DataArrayView::ShapeType s=getDataPointShape();
1125  #endif       DataArrayView::ShapeType ev_shape;
1126    Trace trace_func;       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1127    return dp_algorithm(trace_func,0);       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1128         int rank=getDataPointRank();
1129         if (rank<2) {
1130            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1131         }
1132         if (axis0<0 || axis0>rank-1) {
1133            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1134         }
1135         if (axis1<0 || axis1>rank-1) {
1136             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1137         }
1138         if (axis0 == axis1) {
1139             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1140         }
1141         if (axis0 > axis1) {
1142             axis0_tmp=axis1;
1143             axis1_tmp=axis0;
1144         } else {
1145             axis0_tmp=axis0;
1146             axis1_tmp=axis1;
1147         }
1148         for (int i=0; i<rank; i++) {
1149           if (i == axis0_tmp) {
1150              ev_shape.push_back(s[axis1_tmp]);
1151           } else if (i == axis1_tmp) {
1152              ev_shape.push_back(s[axis0_tmp]);
1153           } else {
1154              ev_shape.push_back(s[i]);
1155           }
1156         }
1157         Data ev(0.,ev_shape,getFunctionSpace());
1158         ev.typeMatchRight(*this);
1159         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1160         return ev;
1161    
1162  }  }
1163    
1164  Data  Data
1165  Data::symmetric() const  Data::symmetric() const
1166  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1167       // check input       // check input
1168       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1169       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1170          if(s[0] != s[1])          if(s[0] != s[1])
1171             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.");
1172       }       }
1173       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
# Line 1353  Data::symmetric() const Line 1186  Data::symmetric() const
1186  Data  Data
1187  Data::nonsymmetric() const  Data::nonsymmetric() const
1188  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1189       // check input       // check input
1190       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1191       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1192          if(s[0] != s[1])          if(s[0] != s[1])
1193             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.");
1194          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1195          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
# Line 1388  Data::nonsymmetric() const Line 1218  Data::nonsymmetric() const
1218  }  }
1219    
1220  Data  Data
1221  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1222  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1223       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1224       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1225          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1226          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1227          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1228          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1229          return ev;          return ev;
1230       }       }
1231       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1240  Data::matrixtrace(int axis_offset) const
1240          }          }
1241          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1242          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1243          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1244          return ev;          return ev;
1245       }       }
1246       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1259  Data::matrixtrace(int axis_offset) const
1259      }      }
1260          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1261          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1262      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1263          return ev;          return ev;
1264       }       }
1265       else {       else {
1266          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.");
1267       }       }
1268  }  }
1269    
1270  Data  Data
1271  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1272  {  {
 #if defined DOPROF  
      profData->reduction2++;  
 #endif  
1273       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1274       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1275       // 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]
# Line 1467  Data::transpose(int axis_offset) const Line 1291  Data::transpose(int axis_offset) const
1291  Data  Data
1292  Data::eigenvalues() const  Data::eigenvalues() const
1293  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1294       // check input       // check input
1295       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1296       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1297          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.");
1298       if(s[0] != s[1])       if(s[0] != s[1])
1299          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.");
1300       // create return       // create return
1301       DataArrayView::ShapeType ev_shape(1,s[0]);       DataArrayView::ShapeType ev_shape(1,s[0]);
# Line 1487  Data::eigenvalues() const Line 1308  Data::eigenvalues() const
1308  const boost::python::tuple  const boost::python::tuple
1309  Data::eigenvalues_and_eigenvectors(const double tol) const  Data::eigenvalues_and_eigenvectors(const double tol) const
1310  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1311       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1312       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1313          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.");
1314       if(s[0] != s[1])       if(s[0] != s[1])
1315          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.");
1316       // create return       // create return
1317       DataArrayView::ShapeType ev_shape(1,s[0]);       DataArrayView::ShapeType ev_shape(1,s[0]);
# Line 1507  Data::eigenvalues_and_eigenvectors(const Line 1325  Data::eigenvalues_and_eigenvectors(const
1325  }  }
1326    
1327  const boost::python::tuple  const boost::python::tuple
1328  Data::mindp() const  Data::minGlobalDataPoint() const
1329  {  {
1330    // 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
1331    // 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
1332    // surrounding function    // surrounding function
1333    
   int SampleNo;  
1334    int DataPointNo;    int DataPointNo;
1335      int ProcNo;    int ProcNo;
1336    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1337    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1338  }  }
1339    
1340  void  void
1341  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1342                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1343  {  {
1344    int i,j;    int i,j;
1345    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1565  Data::calc_mindp(  int& ProcNo, Line 1381  Data::calc_mindp(  int& ProcNo,
1381      int lowProc = 0;      int lowProc = 0;
1382      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1383      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() );
1384        
1385      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1386          next = globalMins[lowProc];          next = globalMins[lowProc];
1387          for( i=1; i<get_MPISize(); i++ )          for( i=1; i<get_MPISize(); i++ )
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1397  Data::calc_mindp(  int& ProcNo,
1397  #else  #else
1398      ProcNo = 0;      ProcNo = 0;
1399  #endif  #endif
1400    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1401  }  }
1402    
1403  void  void
# Line 1609  Data::operator+=(const Data& right) Line 1424  Data::operator+=(const Data& right)
1424    if (isProtected()) {    if (isProtected()) {
1425          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1426    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1427    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1428    return (*this);    return (*this);
1429  }  }
# Line 1619  Data::operator+=(const Data& right) Line 1431  Data::operator+=(const Data& right)
1431  Data&  Data&
1432  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1433  {  {
1434    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1435          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,plus<double>());
1436    }    return (*this);
1437  #if defined DOPROF  }
1438    profData->binary++;  Data&
1439  #endif  Data::operator=(const Data& other)
1440    binaryOp(right,plus<double>());  {
1441      copy(other);
1442    return (*this);    return (*this);
1443  }  }
1444    
# Line 1635  Data::operator-=(const Data& right) Line 1448  Data::operator-=(const Data& right)
1448    if (isProtected()) {    if (isProtected()) {
1449          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1450    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1451    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1452    return (*this);    return (*this);
1453  }  }
# Line 1645  Data::operator-=(const Data& right) Line 1455  Data::operator-=(const Data& right)
1455  Data&  Data&
1456  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1457  {  {
1458    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1459          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1460    return (*this);    return (*this);
1461  }  }
1462    
# Line 1661  Data::operator*=(const Data& right) Line 1466  Data::operator*=(const Data& right)
1466    if (isProtected()) {    if (isProtected()) {
1467          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1468    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1469    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1470    return (*this);    return (*this);
1471  }  }
# Line 1671  Data::operator*=(const Data& right) Line 1473  Data::operator*=(const Data& right)
1473  Data&  Data&
1474  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1475  {  {
1476    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1477          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1478    return (*this);    return (*this);
1479  }  }
1480    
# Line 1687  Data::operator/=(const Data& right) Line 1484  Data::operator/=(const Data& right)
1484    if (isProtected()) {    if (isProtected()) {
1485          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1486    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1487    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1488    return (*this);    return (*this);
1489  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1491  Data::operator/=(const Data& right)
1491  Data&  Data&
1492  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1493  {  {
1494    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1495          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1496    return (*this);    return (*this);
1497  }  }
1498    
1499  Data  Data
1500  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1501  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1502    Data left_d(left,*this);    Data left_d(left,*this);
1503    return left_d.powD(*this);    return left_d.powD(*this);
1504  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1506  Data::rpowO(const boost::python::object&
1506  Data  Data
1507  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1508  {  {
1509  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1510    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1511  }  }
1512    
1513  Data  Data
1514  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1515  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1516    Data result;    Data result;
1517    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1518    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1519         result.binaryOp(*this,escript::rpow);
1520      } else {
1521         result.copy(*this);
1522         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1523      }
1524    return result;    return result;
1525  }  }
1526    
# Line 1753  escript::operator+(const Data& left, con Line 1533  escript::operator+(const Data& left, con
1533    Data result;    Data result;
1534    //    //
1535    // perform a deep copy    // perform a deep copy
1536    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1537    result+=right;       result.copy(right);
1538         result+=left;
1539      } else {
1540         result.copy(left);
1541         result+=right;
1542      }
1543    return result;    return result;
1544  }  }
1545    
# Line 1766  escript::operator-(const Data& left, con Line 1551  escript::operator-(const Data& left, con
1551    Data result;    Data result;
1552    //    //
1553    // perform a deep copy    // perform a deep copy
1554    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1555    result-=right;       result=right.neg();
1556         result+=left;
1557      } else {
1558         result.copy(left);
1559         result-=right;
1560      }
1561    return result;    return result;
1562  }  }
1563    
# Line 1779  escript::operator*(const Data& left, con Line 1569  escript::operator*(const Data& left, con
1569    Data result;    Data result;
1570    //    //
1571    // perform a deep copy    // perform a deep copy
1572    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1573    result*=right;       result.copy(right);
1574         result*=left;
1575      } else {
1576         result.copy(left);
1577         result*=right;
1578      }
1579    return result;    return result;
1580  }  }
1581    
# Line 1792  escript::operator/(const Data& left, con Line 1587  escript::operator/(const Data& left, con
1587    Data result;    Data result;
1588    //    //
1589    // perform a deep copy    // perform a deep copy
1590    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1591    result/=right;       result=right.oneOver();
1592         result*=left;
1593      } else {
1594         result.copy(left);
1595         result/=right;
1596      }
1597    return result;    return result;
1598  }  }
1599    
# Line 1802  escript::operator/(const Data& left, con Line 1602  escript::operator/(const Data& left, con
1602  Data  Data
1603  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1604  {  {
1605    //    return left+Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
1606  }  }
1607    
1608  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1610  escript::operator+(const Data& left, con
1610  Data  Data
1611  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1612  {  {
1613    //    return left-Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
1614  }  }
1615    
1616  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1618  escript::operator-(const Data& left, con
1618  Data  Data
1619  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1620  {  {
1621    //    return left*Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
1622  }  }
1623    
1624  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1626  escript::operator*(const Data& left, con
1626  Data  Data
1627  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1628  {  {
1629    //    return left/Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
1630  }  }
1631    
1632  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1634  escript::operator/(const Data& left, con
1634  Data  Data
1635  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1636  {  {
1637    //    return Data(left,right.getFunctionSpace(),false)+right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result+=right;  
   return result;  
1638  }  }
1639    
1640  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1642  escript::operator+(const boost::python::
1642  Data  Data
1643  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1644  {  {
1645    //    return Data(left,right.getFunctionSpace(),false)-right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result-=right;  
   return result;  
1646  }  }
1647    
1648  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1650  escript::operator-(const boost::python::
1650  Data  Data
1651  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1652  {  {
1653    //    return Data(left,right.getFunctionSpace(),false)*right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result*=right;  
   return result;  
1654  }  }
1655    
1656  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1658  escript::operator*(const boost::python::
1658  Data  Data
1659  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1660  {  {
1661    //    return Data(left,right.getFunctionSpace(),false)/right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result/=right;  
   return result;  
1662  }  }
1663    
1664  //  //
# Line 1961  escript::operator/(const boost::python:: Line 1709  escript::operator/(const boost::python::
1709  /* TODO */  /* TODO */
1710  /* global reduction */  /* global reduction */
1711  Data  Data
1712  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1713  {  {
1714    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1715    
# Line 1979  Data::getItem(const boost::python::objec Line 1727  Data::getItem(const boost::python::objec
1727  Data  Data
1728  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1729  {  {
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
1730    return Data(*this,region);    return Data(*this,region);
1731  }  }
1732    
# Line 1995  Data::setItemO(const boost::python::obje Line 1740  Data::setItemO(const boost::python::obje
1740    setItemD(key,tempData);    setItemD(key,tempData);
1741  }  }
1742    
 /* TODO */  
 /* global reduction */  
1743  void  void
1744  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1745                 const Data& value)                 const Data& value)
# Line 2014  Data::setItemD(const boost::python::obje Line 1757  Data::setItemD(const boost::python::obje
1757    }    }
1758  }  }
1759    
 /* TODO */  
 /* global reduction */  
1760  void  void
1761  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1762                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
# Line 2023  Data::setSlice(const Data& value, Line 1764  Data::setSlice(const Data& value,
1764    if (isProtected()) {    if (isProtected()) {
1765          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1766    }    }
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
1767    Data tempValue(value);    Data tempValue(value);
1768    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1769    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2060  Data::typeMatchRight(const Data& right) Line 1798  Data::typeMatchRight(const Data& right)
1798    }    }
1799  }  }
1800    
1801  /* TODO */  void
1802  /* global reduction */  Data::setTaggedValueByName(std::string name,
1803                               const boost::python::object& value)
1804    {
1805         if (getFunctionSpace().getDomain().isValidTagName(name)) {
1806            int tagKey=getFunctionSpace().getDomain().getTag(name);
1807            setTaggedValue(tagKey,value);
1808         }
1809    }
1810  void  void
1811  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1812                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2077  Data::setTaggedValue(int tagKey, Line 1822  Data::setTaggedValue(int tagKey,
1822      throw DataException("Error - DataTagged conversion failed!!");      throw DataException("Error - DataTagged conversion failed!!");
1823    }    }
1824    
1825    //    numeric::array asNumArray(value);
1826    // Construct DataArray from boost::python::object input value  
1827    DataArray valueDataArray(value);  
1828      // extract the shape of the numarray
1829      DataArrayView::ShapeType tempShape;
1830      for (int i=0; i < asNumArray.getrank(); i++) {
1831        tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
1832      }
1833    
1834      // get the space for the data vector
1835      int len = DataArrayView::noValues(tempShape);
1836      DataVector temp_data(len, 0.0, len);
1837      DataArrayView temp_dataView(temp_data, tempShape);
1838      temp_dataView.copy(asNumArray);
1839    
1840    //    //
1841    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1842    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,temp_dataView);
1843  }  }
1844    
 /* TODO */  
 /* global reduction */  
1845  void  void
1846  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1847                              const DataArrayView& value)                              const DataArrayView& value)
# Line 2102  Data::setTaggedValueFromCPP(int tagKey, Line 1856  Data::setTaggedValueFromCPP(int tagKey,
1856    if (!isTagged()) {    if (!isTagged()) {
1857      throw DataException("Error - DataTagged conversion failed!!");      throw DataException("Error - DataTagged conversion failed!!");
1858    }    }
1859                                                                                                                  
1860    //    //
1861    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1862    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1863  }  }
1864    
 /* TODO */  
 /* global reduction */  
1865  int  int
1866  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
1867  {  {
1868    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
1869  }  }
1870    
 /* 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.");  
   }  
   //  
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 /* 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);  
   
   //  
   // Load values from valueDataArray into return numarray  
   
   // extract the shape of the numarray  
   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();  
   
   if (rank==0) {  
       boost::python::numeric::array temp_numArray(valueView());  
       value = temp_numArray;  
   }  
   if (rank==1) {  
     for (int i=0; i < shape[0]; i++) {  
       value[i] = valueView(i);  
     }  
   }  
   if (rank==2) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         value[i][j] = valueView(i,j);  
       }  
     }  
   }  
   if (rank==3) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           value[i][j][k] = valueView(i,j,k);  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           for (int l=0; l < shape[3]; l++) {  
             value[i][j][k][l] = valueView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
1871  void  void
1872  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
1873  {  {
# Line 2241  Data::archiveData(const std::string file Line 1909  Data::archiveData(const std::string file
1909    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
1910    vector<int> referenceNumbers(noSamples);    vector<int> referenceNumbers(noSamples);
1911    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1912      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1913    }    }
1914    vector<int> tagNumbers(noSamples);    vector<int> tagNumbers(noSamples);
1915    if (isTagged()) {    if (isTagged()) {
# Line 2450  Data::extractData(const std::string file Line 2118  Data::extractData(const std::string file
2118      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2119    }    }
2120    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2121      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2122        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2123      }      }
2124    }    }
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2209  ostream& escript::operator<<(ostream& o,
2209    return o;    return o;
2210  }  }
2211    
2212    Data
2213    escript::C_GeneralTensorProduct(Data& arg_0,
2214                         Data& arg_1,
2215                         int axis_offset,
2216                         int transpose)
2217    {
2218      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2219      // SM is the product of the last axis_offset entries in arg_0.getShape().
2220    
2221      // Interpolate if necessary and find an appropriate function space
2222      Data arg_0_Z, arg_1_Z;
2223      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2224        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2225          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2226          arg_1_Z = Data(arg_1);
2227        }
2228        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2229          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2230          arg_0_Z =Data(arg_0);
2231        }
2232        else {
2233          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2234        }
2235      } else {
2236          arg_0_Z = Data(arg_0);
2237          arg_1_Z = Data(arg_1);
2238      }
2239      // Get rank and shape of inputs
2240      int rank0 = arg_0_Z.getDataPointRank();
2241      int rank1 = arg_1_Z.getDataPointRank();
2242      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2243      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2244    
2245      // Prepare for the loops of the product and verify compatibility of shapes
2246      int start0=0, start1=0;
2247      if (transpose == 0)       {}
2248      else if (transpose == 1)  { start0 = axis_offset; }
2249      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2250      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2251    
2252      // Adjust the shapes for transpose
2253      DataArrayView::ShapeType tmpShape0;
2254      DataArrayView::ShapeType tmpShape1;
2255      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2256      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2257    
2258    #if 0
2259      // For debugging: show shape after transpose
2260      char tmp[100];
2261      std::string shapeStr;
2262      shapeStr = "(";
2263      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2264      shapeStr += ")";
2265      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2266      shapeStr = "(";
2267      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2268      shapeStr += ")";
2269      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2270    #endif
2271    
2272      // Prepare for the loops of the product
2273      int SL=1, SM=1, SR=1;
2274      for (int i=0; i<rank0-axis_offset; i++)   {
2275        SL *= tmpShape0[i];
2276      }
2277      for (int i=rank0-axis_offset; i<rank0; i++)   {
2278        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2279          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2280        }
2281        SM *= tmpShape0[i];
2282      }
2283      for (int i=axis_offset; i<rank1; i++)     {
2284        SR *= tmpShape1[i];
2285      }
2286    
2287      // Define the shape of the output
2288      DataArrayView::ShapeType shape2;
2289      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2290      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2291    
2292      // Declare output Data object
2293      Data res;
2294    
2295      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2296        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2297        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2298        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2299        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2300        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2301      }
2302      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2303    
2304        // Prepare the DataConstant input
2305        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2306        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2307    
2308        // Borrow DataTagged input from Data object
2309        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2310        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2311    
2312        // Prepare a DataTagged output 2
2313        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2314        res.tag();
2315        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2316        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2317    
2318        // Prepare offset into DataConstant
2319        int offset_0 = tmp_0->getPointOffset(0,0);
2320        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2321        // Get the views
2322        DataArrayView view_1 = tmp_1->getDefaultValue();
2323        DataArrayView view_2 = tmp_2->getDefaultValue();
2324        // Get the pointers to the actual data
2325        double *ptr_1 = &((view_1.getData())[0]);
2326        double *ptr_2 = &((view_2.getData())[0]);
2327        // Compute an MVP for the default
2328        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2329        // Compute an MVP for each tag
2330        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2331        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2332        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2333          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2334          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2335          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2336          double *ptr_1 = &view_1.getData(0);
2337          double *ptr_2 = &view_2.getData(0);
2338          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2339        }
2340    
2341      }
2342      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2343    
2344        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2345        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2346        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2347        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2348        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2349        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2350        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2351        int sampleNo_1,dataPointNo_1;
2352        int numSamples_1 = arg_1_Z.getNumSamples();
2353        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2354        int offset_0 = tmp_0->getPointOffset(0,0);
2355        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2356        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2357          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2358            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2359            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2360            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2361            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2362            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2363            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2364          }
2365        }
2366    
2367      }
2368      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2369    
2370        // Borrow DataTagged input from Data object
2371        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2372        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2373    
2374        // Prepare the DataConstant input
2375        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2376        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2377    
2378        // Prepare a DataTagged output 2
2379        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2380        res.tag();
2381        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2382        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2383    
2384        // Prepare offset into DataConstant
2385        int offset_1 = tmp_1->getPointOffset(0,0);
2386        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2387        // Get the views
2388        DataArrayView view_0 = tmp_0->getDefaultValue();
2389        DataArrayView view_2 = tmp_2->getDefaultValue();
2390        // Get the pointers to the actual data
2391        double *ptr_0 = &((view_0.getData())[0]);
2392        double *ptr_2 = &((view_2.getData())[0]);
2393        // Compute an MVP for the default
2394        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2395        // Compute an MVP for each tag
2396        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2397        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2398        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2399          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2400          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2401          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2402          double *ptr_0 = &view_0.getData(0);
2403          double *ptr_2 = &view_2.getData(0);
2404          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2405        }
2406    
2407      }
2408      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2409    
2410        // Borrow DataTagged input from Data object
2411        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2412        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2413    
2414        // Borrow DataTagged input from Data object
2415        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2416        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2417    
2418        // Prepare a DataTagged output 2
2419        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2420        res.tag();  // DataTagged output
2421        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2422        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2423    
2424        // Get the views
2425        DataArrayView view_0 = tmp_0->getDefaultValue();
2426        DataArrayView view_1 = tmp_1->getDefaultValue();
2427        DataArrayView view_2 = tmp_2->getDefaultValue();
2428        // Get the pointers to the actual data
2429        double *ptr_0 = &((view_0.getData())[0]);
2430        double *ptr_1 = &((view_1.getData())[0]);
2431        double *ptr_2 = &((view_2.getData())[0]);
2432        // Compute an MVP for the default
2433        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2434        // Merge the tags
2435        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2436        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2437        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2438        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2439          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2440        }
2441        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2442          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2443        }
2444        // Compute an MVP for each tag
2445        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2446        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2447          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2448          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2449          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2450          double *ptr_0 = &view_0.getData(0);
2451          double *ptr_1 = &view_1.getData(0);
2452          double *ptr_2 = &view_2.getData(0);
2453          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2454        }
2455    
2456      }
2457      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2458    
2459        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2460        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2461        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2462        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2463        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2464        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2465        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2466        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2467        int sampleNo_0,dataPointNo_0;
2468        int numSamples_0 = arg_0_Z.getNumSamples();
2469        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2470        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2471        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2472          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2473          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2474          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2475            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2476            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2477            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2478            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2479            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2480          }
2481        }
2482    
2483      }
2484      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2485    
2486        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2487        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2488        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2489        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2490        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2491        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2492        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2493        int sampleNo_0,dataPointNo_0;
2494        int numSamples_0 = arg_0_Z.getNumSamples();
2495        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2496        int offset_1 = tmp_1->getPointOffset(0,0);
2497        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2498        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2499          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2500            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2501            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2502            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2503            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2504            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2505            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2506          }
2507        }
2508    
2509    
2510      }
2511      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2512    
2513        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2514        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2515        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2516        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2517        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2518        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2519        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2520        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2521        int sampleNo_0,dataPointNo_0;
2522        int numSamples_0 = arg_0_Z.getNumSamples();
2523        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2524        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2525        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2526          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2527          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2528          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2529            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2530            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2531            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2532            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2533            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2534          }
2535        }
2536    
2537      }
2538      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2539    
2540        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2541        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2542        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2543        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2544        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2545        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2546        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2547        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2548        int sampleNo_0,dataPointNo_0;
2549        int numSamples_0 = arg_0_Z.getNumSamples();
2550        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2551        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2552        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2553          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2554            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2555            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2556            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2557            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2558            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2559            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2560            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2561          }
2562        }
2563    
2564      }
2565      else {
2566        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2567      }
2568    
2569      return res;
2570    }
2571    
2572    DataAbstract*
2573    Data::borrowData() const
2574    {
2575      return m_data.get();
2576    }
2577    
2578  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2579    
2580  void  void
2581  Data::print()  Data::print()
2582  {  {
2583    int i,j;    int i,j;
2584      
2585    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2586    for( i=0; i<getNumSamples(); i++ )    for( i=0; i<getNumSamples(); i++ )
2587    {    {
# Line 2557  Data::print() Line 2591  Data::print()
2591      printf( "\n" );      printf( "\n" );
2592    }    }
2593  }  }
2594    void
2595    Data::dump(const std::string fileName) const
2596    {
2597      try
2598         {
2599            return m_data->dump(fileName);
2600         }
2601         catch (exception& e)
2602         {
2603            cout << e.what() << endl;
2604         }
2605    }
2606    
2607  int  int
2608  Data::get_MPISize() const  Data::get_MPISize() const
# Line 2584  Data::get_MPIRank() const Line 2630  Data::get_MPIRank() const
2630    
2631  MPI_Comm  MPI_Comm
2632  Data::get_MPIComm() const  Data::get_MPIComm() const
2633  {  {
2634  #ifdef PASO_MPI  #ifdef PASO_MPI
2635      return MPI_COMM_WORLD;      return MPI_COMM_WORLD;
2636  #else  #else

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

  ViewVC Help
Powered by ViewVC 1.1.26