/[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 711 by gross, Wed Apr 26 22:39:51 2006 UTC revision 1092 by gross, Fri Apr 13 03:39:49 2007 UTC
# Line 11  Line 11 
11   *                                                          *   *                                                          *
12   ************************************************************   ************************************************************
13  */  */
   
14  #include "Data.h"  #include "Data.h"
15    
16  #include "DataExpanded.h"  #include "DataExpanded.h"
# Line 29  Line 28 
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
# Line 40  using namespace boost::python; Line 38  using namespace boost::python;
38  using namespace boost;  using namespace boost;
39  using namespace escript;  using namespace escript;
40    
 #if defined DOPROF  
 //  
 // global table of profiling data for all Data objects  
 DataProf dataProfTable;  
 #endif  
   
41  Data::Data()  Data::Data()
42  {  {
43    //    //
# Line 53  Data::Data() Line 45  Data::Data()
45    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
46    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
47    m_data=temp_data;    m_data=temp_data;
48  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
49  }  }
50    
51  Data::Data(double value,  Data::Data(double value,
# Line 70  Data::Data(double value, Line 59  Data::Data(double value,
59    }    }
60    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
61    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
62  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
63  }  }
64    
65  Data::Data(double value,  Data::Data(double value,
# Line 84  Data::Data(double value, Line 70  Data::Data(double value,
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
72    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
73  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
74  }  }
75    
76  Data::Data(const Data& inData)  Data::Data(const Data& inData)
77  {  {
78    m_data=inData.m_data;    m_data=inData.m_data;
79  #if defined DOPROF    m_protected=inData.isProtected();
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
80  }  }
81    
82  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 107  Data::Data(const Data& inData, Line 87  Data::Data(const Data& inData,
87    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
88    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
89    m_data=temp_data;    m_data=temp_data;
90  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
91  }  }
92    
93  Data::Data(const Data& inData,  Data::Data(const Data& inData,
94             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
95  {  {
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
96    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
97      m_data=inData.m_data;      m_data=inData.m_data;
98    } else {    } else {
     #if defined DOPROF  
     profData->interpolate++;  
     #endif  
99      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
100      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
101      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 139  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109      }      }
110      m_data=tmp.m_data;      m_data=tmp.m_data;
111    }    }
112      m_protected=false;
113  }  }
114    
115  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 150  Data::Data(const DataTagged::TagListType Line 121  Data::Data(const DataTagged::TagListType
121    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
122    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
123    m_data=temp_data;    m_data=temp_data;
124      m_protected=false;
125    if (expanded) {    if (expanded) {
126      expand();      expand();
127    }    }
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
128  }  }
129    
130  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 164  Data::Data(const numeric::array& value, Line 132  Data::Data(const numeric::array& value,
132             bool expanded)             bool expanded)
133  {  {
134    initialise(value,what,expanded);    initialise(value,what,expanded);
135  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
136  }  }
137    
138  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 175  Data::Data(const DataArrayView& value, Line 140  Data::Data(const DataArrayView& value,
140             bool expanded)             bool expanded)
141  {  {
142    initialise(value,what,expanded);    initialise(value,what,expanded);
143  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
144  }  }
145    
146  Data::Data(const object& value,  Data::Data(const object& value,
# Line 187  Data::Data(const object& value, Line 149  Data::Data(const object& value,
149  {  {
150    numeric::array asNumArray(value);    numeric::array asNumArray(value);
151    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
152  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
153  }  }
154    
155  Data::Data(const object& value,  Data::Data(const object& value,
# Line 211  Data::Data(const object& value, Line 170  Data::Data(const object& value,
170      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
171      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
172    }    }
173  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
174  }  }
175    
176  Data::~Data()  Data::~Data()
# Line 356  Data::isConstant() const Line 312  Data::isConstant() const
312  }  }
313    
314  void  void
315    Data::setProtection()
316    {
317       m_protected=true;
318    }
319    
320    bool
321    Data::isProtected() const
322    {
323       return m_protected;
324    }
325    
326    
327    
328    void
329  Data::expand()  Data::expand()
330  {  {
331    if (isConstant()) {    if (isConstant()) {
# Line 397  Data::tag() Line 367  Data::tag()
367    }    }
368  }  }
369    
370  void  Data
371  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
372  {  {
373    m_data->reshapeDataPoint(shape);    return escript::unaryOp(*this,bind1st(divides<double>(),1.));
374  }  }
375    
376  Data  Data
377  Data::wherePositive() const  Data::wherePositive() const
378  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
379    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
380  }  }
381    
382  Data  Data
383  Data::whereNegative() const  Data::whereNegative() const
384  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
385    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
386  }  }
387    
388  Data  Data
389  Data::whereNonNegative() const  Data::whereNonNegative() const
390  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
391    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
392  }  }
393    
394  Data  Data
395  Data::whereNonPositive() const  Data::whereNonPositive() const
396  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
397    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
398  }  }
399    
400  Data  Data
401  Data::whereZero(double tol) const  Data::whereZero(double tol) const
402  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
403    Data dataAbs=abs();    Data dataAbs=abs();
404    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
405  }  }
# Line 452  Data::whereZero(double tol) const Line 407  Data::whereZero(double tol) const
407  Data  Data
408  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
409  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
410    Data dataAbs=abs();    Data dataAbs=abs();
411    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
412  }  }
# Line 462  Data::whereNonZero(double tol) const Line 414  Data::whereNonZero(double tol) const
414  Data  Data
415  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
416  {  {
 #if defined DOPROF  
   profData->interpolate++;  
 #endif  
417    return Data(*this,functionspace);    return Data(*this,functionspace);
418  }  }
419    
# Line 486  Data::probeInterpolation(const FunctionS Line 435  Data::probeInterpolation(const FunctionS
435  Data  Data
436  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
437  {  {
 #if defined DOPROF  
   profData->grad++;  
 #endif  
438    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
439      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
440    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 522  Data::getDataPointShape() const Line 468  Data::getDataPointShape() const
468    return getPointDataView().getShape();    return getPointDataView().getShape();
469  }  }
470    
 void  
 Data::fillFromNumArray(const boost::python::numeric::array num_array)  
 {  
   //  
   // check rank  
   if (num_array.getrank()<getDataPointRank())  
       throw DataException("Rank of numarray does not match Data object rank");  
   
   //  
   // check shape of num_array  
   for (int i=0; i<getDataPointRank(); i++) {  
     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])  
        throw DataException("Shape of numarray does not match Data object rank");  
   }  
   
   //  
   // make sure data is expanded:  
   if (!isExpanded()) {  
     expand();  
   }  
471    
   //  
   // and copy over  
   m_data->copyAll(num_array);  
 }  
472    
473  const  const
474  boost::python::numeric::array  boost::python::numeric::array
475  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
476  {  {
477    //    size_t length=0;
478    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
479    //    //
480    // determine the rank and shape of each data point    // determine the rank and shape of each data point
481    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 568  Data::convertToNumArray() Line 486  Data::convertToNumArray()
486    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
487    
488    //    //
489    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
490    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
491    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
492    // 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]);  
   }  
493    
494    //    //
495    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
496      if (arrayRank==0) {
497        numArray.resize(1);
498      }
499    if (arrayRank==1) {    if (arrayRank==1) {
500      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
501    }    }
# Line 595  Data::convertToNumArray() Line 508  Data::convertToNumArray()
508    if (arrayRank==4) {    if (arrayRank==4) {
509      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
510    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
511    
512    //    if (getNumDataPointsPerSample()>0) {
513    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
514    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
515    int dataPoint = 0;         //
516    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
517      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
518        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
519        if (dataPointRank==0) {         }
520          numArray[dataPoint]=dataPointView();                
521        }         //
522        if (dataPointRank==1) {         // Check a valid data point number has been supplied
523          for (int i=0; i<dataPointShape[0]; i++) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
524            numArray[dataPoint][i]=dataPointView(i);             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
525          }         }
526        }         // TODO: global error handling
527        if (dataPointRank==2) {         // create a view of the data if it is stored locally
528          for (int i=0; i<dataPointShape[0]; i++) {         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
529            for (int j=0; j<dataPointShape[1]; j++) {          
530              numArray[dataPoint][i][j] = dataPointView(i,j);         switch( dataPointRank ){
531            }              case 0 :
532          }                  numArray[0] = dataPointView();
533        }                  break;
534        if (dataPointRank==3) {              case 1 :        
535          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
536            for (int j=0; j<dataPointShape[1]; j++) {                      numArray[i]=dataPointView(i);
537              for (int k=0; k<dataPointShape[2]; k++) {                  break;
538                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);              case 2 :        
539              }                  for( i=0; i<dataPointShape[0]; i++ )
540            }                      for( j=0; j<dataPointShape[1]; j++)
541          }                          numArray[make_tuple(i,j)]=dataPointView(i,j);
542        }                  break;
543        if (dataPointRank==4) {              case 3 :        
544          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
545            for (int j=0; j<dataPointShape[1]; j++) {                      for( j=0; j<dataPointShape[1]; j++ )
546              for (int k=0; k<dataPointShape[2]; k++) {                          for( k=0; k<dataPointShape[2]; k++)
547                for (int l=0; l<dataPointShape[3]; l++) {                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
548                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  break;
549                }              case 4 :
550              }                  for( i=0; i<dataPointShape[0]; i++ )
551            }                      for( j=0; j<dataPointShape[1]; j++ )
552          }                          for( k=0; k<dataPointShape[2]; k++ )
553        }                              for( l=0; l<dataPointShape[3]; l++)
554        dataPoint++;                                  numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
555      }                  break;
556        }
557    }    }
   
558    //    //
559    // return the loaded array    // return the array
560    return numArray;    return numArray;
 }  
561    
562  const  }
563  boost::python::numeric::array  void
564  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
565  {  {
566    //      // this will throw if the value cannot be represented
567    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
568    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
569    
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
570    
571    //  }
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
572    
573    void
574    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
575    {
576      if (isProtected()) {
577            throw DataException("Error - attempt to update protected Data object.");
578      }
579    //    //
580    // create the numeric array to be returned    // check rank
581    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
582          throw DataException("Rank of numarray does not match Data object rank");
583    
584    //    //
585    // the rank of the returned numeric array will be the rank of    // check shape of num_array
586    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
587    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
588    // 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]);  
589    }    }
   
590    //    //
591    // resize the numeric array to the shape just calculated    // make sure data is expanded:
592    if (arrayRank==1) {    if (!isExpanded()) {
593      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]);  
   }  
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
594    }    }
595      if (getNumDataPointsPerSample()>0) {
596    //         int sampleNo = dataPointNo/getNumDataPointsPerSample();
597    // loop through each data point in turn, loading the values for that data point         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
598    // into the numeric array.         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
599    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    } else {
600      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         m_data->copyToDataPoint(-1, 0,num_array);
     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);  
             }  
           }  
         }  
       }  
     }  
601    }    }
   
   //  
   // return the loaded array  
   return numArray;  
602  }  }
603    
604  const  void
605  boost::python::numeric::array  Data::setValueOfDataPoint(int dataPointNo, const double value)
 Data::convertToNumArrayFromDPNo(int sampleNo,  
                                 int dataPointNo)  
606  {  {
607    //    if (isProtected()) {
608    // Check a valid sample number has been supplied          throw DataException("Error - attempt to update protected Data object.");
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
609    }    }
   
610    //    //
611    // Check a valid data point number has been supplied    // make sure data is expanded:
612    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
613      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
614    }    }
615      if (getNumDataPointsPerSample()>0) {
616           int sampleNo = dataPointNo/getNumDataPointsPerSample();
617           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
618           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
619      } else {
620           m_data->copyToDataPoint(-1, 0,value);
621      }
622    }
623    
624    const
625    boost::python::numeric::array
626    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
627    {
628      size_t length=0;
629      int i, j, k, l, pos;
630    //    //
631    // determine the rank and shape of each data point    // determine the rank and shape of each data point
632    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 802  Data::convertToNumArrayFromDPNo(int samp Line 660  Data::convertToNumArrayFromDPNo(int samp
660      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
661    }    }
662    
663      // added for the MPI communication
664      length=1;
665      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
666      double *tmpData = new double[length];
667    
668    //    //
669    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
   DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
   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);  
           }  
         }  
       }  
     }  
   }  
670    
671        // updated for the MPI case
672        if( get_MPIRank()==procNo ){
673                 if (getNumDataPointsPerSample()>0) {
674                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
675                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
676                    //
677                    // Check a valid sample number has been supplied
678                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
679                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
680                    }
681                  
682                    //
683                    // Check a valid data point number has been supplied
684                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
685                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
686                    }
687                    // TODO: global error handling
688            // create a view of the data if it is stored locally
689            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
690            
691            // pack the data from the view into tmpData for MPI communication
692            pos=0;
693            switch( dataPointRank ){
694                case 0 :
695                    tmpData[0] = dataPointView();
696                    break;
697                case 1 :        
698                    for( i=0; i<dataPointShape[0]; i++ )
699                        tmpData[i]=dataPointView(i);
700                    break;
701                case 2 :        
702                    for( i=0; i<dataPointShape[0]; i++ )
703                        for( j=0; j<dataPointShape[1]; j++, pos++ )
704                            tmpData[pos]=dataPointView(i,j);
705                    break;
706                case 3 :        
707                    for( i=0; i<dataPointShape[0]; i++ )
708                        for( j=0; j<dataPointShape[1]; j++ )
709                            for( k=0; k<dataPointShape[2]; k++, pos++ )
710                                tmpData[pos]=dataPointView(i,j,k);
711                    break;
712                case 4 :
713                    for( i=0; i<dataPointShape[0]; i++ )
714                        for( j=0; j<dataPointShape[1]; j++ )
715                            for( k=0; k<dataPointShape[2]; k++ )
716                                for( l=0; l<dataPointShape[3]; l++, pos++ )
717                                    tmpData[pos]=dataPointView(i,j,k,l);
718                    break;
719            }
720                }
721        }
722            #ifdef PASO_MPI
723            // broadcast the data to all other processes
724        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
725            #endif
726    
727        // unpack the data
728        switch( dataPointRank ){
729            case 0 :
730                numArray[0]=tmpData[0];
731                break;
732            case 1 :        
733                for( i=0; i<dataPointShape[0]; i++ )
734                    numArray[i]=tmpData[i];
735                break;
736            case 2 :        
737                for( i=0; i<dataPointShape[0]; i++ )
738                    for( j=0; j<dataPointShape[1]; j++ )
739                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
740                break;
741            case 3 :        
742                for( i=0; i<dataPointShape[0]; i++ )
743                    for( j=0; j<dataPointShape[1]; j++ )
744                        for( k=0; k<dataPointShape[2]; k++ )
745                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
746                break;
747            case 4 :
748                for( i=0; i<dataPointShape[0]; i++ )
749                    for( j=0; j<dataPointShape[1]; j++ )
750                        for( k=0; k<dataPointShape[2]; k++ )
751                            for( l=0; l<dataPointShape[3]; l++ )
752                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
753                break;
754        }
755    
756        delete [] tmpData;  
757    //    //
758    // return the loaded array    // return the loaded array
759    return numArray;    return numArray;
760  }  }
761    
762    
763    
764  boost::python::numeric::array  boost::python::numeric::array
765  Data::integrate() const  Data::integrate() const
766  {  {
# Line 853  Data::integrate() const Line 768  Data::integrate() const
768    int rank = getDataPointRank();    int rank = getDataPointRank();
769    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
770    
 #if defined DOPROF  
   profData->integrate++;  
 #endif  
771    
772    //    //
773    // calculate the integral values    // calculate the integral values
# Line 920  Data::integrate() const Line 832  Data::integrate() const
832  Data  Data
833  Data::sin() const  Data::sin() const
834  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
835    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
836  }  }
837    
838  Data  Data
839  Data::cos() const  Data::cos() const
840  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
841    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
842  }  }
843    
844  Data  Data
845  Data::tan() const  Data::tan() const
846  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
847    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
848  }  }
849    
850  Data  Data
851  Data::asin() const  Data::asin() const
852  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
853    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
854  }  }
855    
856  Data  Data
857  Data::acos() const  Data::acos() const
858  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
859    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
860  }  }
861    
862    
863  Data  Data
864  Data::atan() const  Data::atan() const
865  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
866    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
867  }  }
868    
869  Data  Data
870  Data::sinh() const  Data::sinh() const
871  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
872    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
873  }  }
874    
875  Data  Data
876  Data::cosh() const  Data::cosh() const
877  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
878    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
879  }  }
880    
881  Data  Data
882  Data::tanh() const  Data::tanh() const
883  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
884    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
885  }  }
886    
887    
888  Data  Data
889  Data::asinh() const  Data::erf() const
890  {  {
891  #if defined DOPROF  #ifdef _WIN32
892    profData->unary++;    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893    #else
894      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
895  #endif  #endif
896    }
897    
898    Data
899    Data::asinh() const
900    {
901    #ifdef _WIN32
902      return escript::unaryOp(*this,escript::asinh_substitute);
903    #else
904    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
905    #endif
906  }  }
907    
908  Data  Data
909  Data::acosh() const  Data::acosh() const
910  {  {
911  #if defined DOPROF  #ifdef _WIN32
912    profData->unary++;    return escript::unaryOp(*this,escript::acosh_substitute);
913  #endif  #else
914    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
915    #endif
916  }  }
917    
918  Data  Data
919  Data::atanh() const  Data::atanh() const
920  {  {
921  #if defined DOPROF  #ifdef _WIN32
922    profData->unary++;    return escript::unaryOp(*this,escript::atanh_substitute);
923  #endif  #else
924    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
925    #endif
926  }  }
927    
928  Data  Data
929  Data::log10() const  Data::log10() const
930  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
931    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
932  }  }
933    
934  Data  Data
935  Data::log() const  Data::log() const
936  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
937    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
938  }  }
939    
940  Data  Data
941  Data::sign() const  Data::sign() const
942  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
943    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
944  }  }
945    
946  Data  Data
947  Data::abs() const  Data::abs() const
948  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
949    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
950  }  }
951    
952  Data  Data
953  Data::neg() const  Data::neg() const
954  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
955    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
956  }  }
957    
958  Data  Data
959  Data::pos() const  Data::pos() const
960  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
961    Data result;    Data result;
962    // perform a deep copy    // perform a deep copy
963    result.copy(*this);    result.copy(*this);
# Line 1085  Data::pos() const Line 967  Data::pos() const
967  Data  Data
968  Data::exp() const  Data::exp() const
969  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
970    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
971  }  }
972    
973  Data  Data
974  Data::sqrt() const  Data::sqrt() const
975  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
976    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
977  }  }
978    
979  double  double
980  Data::Lsup() const  Data::Lsup() const
981  {  {
982  #if defined DOPROF    double localValue, globalValue;
   profData->reduction1++;  
 #endif  
983    //    //
984    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
985    
986    AbsMax abs_max_func;    AbsMax abs_max_func;
987    return algorithm(abs_max_func,0);    localValue = algorithm(abs_max_func,0);
988    #ifdef PASO_MPI
989      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
990      return globalValue;
991    #else
992      return localValue;
993    #endif
994  }  }
995    
996  double  double
997  Data::Linf() const  Data::Linf() const
998  {  {
999  #if defined DOPROF    double localValue, globalValue;
   profData->reduction1++;  
 #endif  
1000    //    //
1001    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1002    AbsMin abs_min_func;    AbsMin abs_min_func;
1003    return algorithm(abs_min_func,numeric_limits<double>::max());    localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1004    
1005    #ifdef PASO_MPI
1006      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1007      return globalValue;
1008    #else
1009      return localValue;
1010    #endif
1011  }  }
1012    
1013  double  double
1014  Data::sup() const  Data::sup() const
1015  {  {
1016  #if defined DOPROF    double localValue, globalValue;
   profData->reduction1++;  
 #endif  
1017    //    //
1018    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1019    FMax fmax_func;    FMax fmax_func;
1020    return algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1021    #ifdef PASO_MPI
1022      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1023      return globalValue;
1024    #else
1025      return localValue;
1026    #endif
1027  }  }
1028    
1029  double  double
1030  Data::inf() const  Data::inf() const
1031  {  {
1032  #if defined DOPROF    double localValue, globalValue;
   profData->reduction1++;  
 #endif  
1033    //    //
1034    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1035    FMin fmin_func;    FMin fmin_func;
1036    return algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1037    #ifdef PASO_MPI
1038      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1039      return globalValue;
1040    #else
1041      return localValue;
1042    #endif
1043  }  }
1044    
1045    /* TODO */
1046    /* global reduction */
1047  Data  Data
1048  Data::maxval() const  Data::maxval() const
1049  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1050    //    //
1051    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1052    FMax fmax_func;    FMax fmax_func;
# Line 1163  Data::maxval() const Line 1056  Data::maxval() const
1056  Data  Data
1057  Data::minval() const  Data::minval() const
1058  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1059    //    //
1060    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1061    FMin fmin_func;    FMin fmin_func;
# Line 1173  Data::minval() const Line 1063  Data::minval() const
1063  }  }
1064    
1065  Data  Data
1066  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1067  {  {
1068  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1069    profData->reduction2++;       DataArrayView::ShapeType s=getDataPointShape();
1070  #endif       DataArrayView::ShapeType ev_shape;
1071    Trace trace_func;       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1072    return dp_algorithm(trace_func,0);       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1073         int rank=getDataPointRank();
1074         if (rank<2) {
1075            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1076         }
1077         if (axis0<0 || axis0>rank-1) {
1078            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1079         }
1080         if (axis1<0 || axis1>rank-1) {
1081             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1082         }
1083         if (axis0 == axis1) {
1084             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1085         }
1086         if (axis0 > axis1) {
1087             axis0_tmp=axis1;
1088             axis1_tmp=axis0;
1089         } else {
1090             axis0_tmp=axis0;
1091             axis1_tmp=axis1;
1092         }
1093         for (int i=0; i<rank; i++) {
1094           if (i == axis0_tmp) {
1095              ev_shape.push_back(s[axis1_tmp]);
1096           } else if (i == axis1_tmp) {
1097              ev_shape.push_back(s[axis0_tmp]);
1098           } else {
1099              ev_shape.push_back(s[i]);
1100           }
1101         }
1102         Data ev(0.,ev_shape,getFunctionSpace());
1103         ev.typeMatchRight(*this);
1104         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1105         return ev;
1106    
1107  }  }
1108    
1109  Data  Data
1110  Data::transpose(int axis) const  Data::symmetric() const
1111  {  {
1112  #if defined DOPROF       // check input
1113    profData->reduction2++;       DataArrayView::ShapeType s=getDataPointShape();
1114  #endif       if (getDataPointRank()==2) {
1115            if(s[0] != s[1])
1116               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1117         }
1118         else if (getDataPointRank()==4) {
1119            if(!(s[0] == s[2] && s[1] == s[3]))
1120               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1121         }
1122         else {
1123            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1124         }
1125         Data ev(0.,getDataPointShape(),getFunctionSpace());
1126         ev.typeMatchRight(*this);
1127         m_data->symmetric(ev.m_data.get());
1128         return ev;
1129    }
1130    
1131    // not implemented  Data
1132    throw DataException("Error - Data::transpose not implemented yet.");  Data::nonsymmetric() const
1133    return Data();  {
1134         // check input
1135         DataArrayView::ShapeType s=getDataPointShape();
1136         if (getDataPointRank()==2) {
1137            if(s[0] != s[1])
1138               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1139            DataArrayView::ShapeType ev_shape;
1140            ev_shape.push_back(s[0]);
1141            ev_shape.push_back(s[1]);
1142            Data ev(0.,ev_shape,getFunctionSpace());
1143            ev.typeMatchRight(*this);
1144            m_data->nonsymmetric(ev.m_data.get());
1145            return ev;
1146         }
1147         else if (getDataPointRank()==4) {
1148            if(!(s[0] == s[2] && s[1] == s[3]))
1149               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1150            DataArrayView::ShapeType ev_shape;
1151            ev_shape.push_back(s[0]);
1152            ev_shape.push_back(s[1]);
1153            ev_shape.push_back(s[2]);
1154            ev_shape.push_back(s[3]);
1155            Data ev(0.,ev_shape,getFunctionSpace());
1156            ev.typeMatchRight(*this);
1157            m_data->nonsymmetric(ev.m_data.get());
1158            return ev;
1159         }
1160         else {
1161            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1162         }
1163    }
1164    
1165    Data
1166    Data::trace(int axis_offset) const
1167    {
1168         DataArrayView::ShapeType s=getDataPointShape();
1169         if (getDataPointRank()==2) {
1170            DataArrayView::ShapeType ev_shape;
1171            Data ev(0.,ev_shape,getFunctionSpace());
1172            ev.typeMatchRight(*this);
1173            m_data->trace(ev.m_data.get(), axis_offset);
1174            return ev;
1175         }
1176         if (getDataPointRank()==3) {
1177            DataArrayView::ShapeType ev_shape;
1178            if (axis_offset==0) {
1179              int s2=s[2];
1180              ev_shape.push_back(s2);
1181            }
1182            else if (axis_offset==1) {
1183              int s0=s[0];
1184              ev_shape.push_back(s0);
1185            }
1186            Data ev(0.,ev_shape,getFunctionSpace());
1187            ev.typeMatchRight(*this);
1188            m_data->trace(ev.m_data.get(), axis_offset);
1189            return ev;
1190         }
1191         if (getDataPointRank()==4) {
1192            DataArrayView::ShapeType ev_shape;
1193            if (axis_offset==0) {
1194              ev_shape.push_back(s[2]);
1195              ev_shape.push_back(s[3]);
1196            }
1197            else if (axis_offset==1) {
1198              ev_shape.push_back(s[0]);
1199              ev_shape.push_back(s[3]);
1200            }
1201        else if (axis_offset==2) {
1202          ev_shape.push_back(s[0]);
1203          ev_shape.push_back(s[1]);
1204        }
1205            Data ev(0.,ev_shape,getFunctionSpace());
1206            ev.typeMatchRight(*this);
1207        m_data->trace(ev.m_data.get(), axis_offset);
1208            return ev;
1209         }
1210         else {
1211            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1212         }
1213    }
1214    
1215    Data
1216    Data::transpose(int axis_offset) const
1217    {
1218         DataArrayView::ShapeType s=getDataPointShape();
1219         DataArrayView::ShapeType ev_shape;
1220         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1221         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1222         int rank=getDataPointRank();
1223         if (axis_offset<0 || axis_offset>rank) {
1224            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1225         }
1226         for (int i=0; i<rank; i++) {
1227           int index = (axis_offset+i)%rank;
1228           ev_shape.push_back(s[index]); // Append to new shape
1229         }
1230         Data ev(0.,ev_shape,getFunctionSpace());
1231         ev.typeMatchRight(*this);
1232         m_data->transpose(ev.m_data.get(), axis_offset);
1233         return ev;
1234  }  }
1235    
1236  Data  Data
1237  Data::eigenvalues() const  Data::eigenvalues() const
1238  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1239       // check input       // check input
1240       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1241       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
# Line 1217  Data::eigenvalues() const Line 1253  Data::eigenvalues() const
1253  const boost::python::tuple  const boost::python::tuple
1254  Data::eigenvalues_and_eigenvectors(const double tol) const  Data::eigenvalues_and_eigenvectors(const double tol) const
1255  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1256       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1257       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1258          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.");
# Line 1237  Data::eigenvalues_and_eigenvectors(const Line 1270  Data::eigenvalues_and_eigenvectors(const
1270  }  }
1271    
1272  const boost::python::tuple  const boost::python::tuple
1273  Data::mindp() const  Data::minGlobalDataPoint() const
1274  {  {
1275    // 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
1276    // 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
1277    // surrounding function    // surrounding function
1278    
   int SampleNo;  
1279    int DataPointNo;    int DataPointNo;
1280      int ProcNo;
1281    calc_mindp(SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1282      return make_tuple(ProcNo,DataPointNo);
   return make_tuple(SampleNo,DataPointNo);  
1283  }  }
1284    
1285  void  void
1286  Data::calc_mindp(int& SampleNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1287                   int& DataPointNo) const                          int& DataPointNo) const
1288  {  {
1289    int i,j;    int i,j;
1290    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1289  Data::calc_mindp(int& SampleNo, Line 1320  Data::calc_mindp(int& SampleNo,
1320      }      }
1321    }    }
1322    
1323    SampleNo = lowi;  #ifdef PASO_MPI
1324    DataPointNo = lowj;      // determine the processor on which the minimum occurs
1325        next = temp.getDataPoint(lowi,lowj)();
1326        int lowProc = 0;
1327        double *globalMins = new double[get_MPISize()+1];
1328        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1329        
1330        if( get_MPIRank()==0 ){
1331            next = globalMins[lowProc];
1332            for( i=1; i<get_MPISize(); i++ )
1333                if( next>globalMins[i] ){
1334                    lowProc = i;
1335                    next = globalMins[i];
1336                }
1337        }
1338        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1339    
1340        delete [] globalMins;
1341        ProcNo = lowProc;
1342    #else
1343        ProcNo = 0;
1344    #endif
1345      DataPointNo = lowj + lowi * numDPPSample;
1346  }  }
1347    
1348  void  void
# Line 1314  Data::saveVTK(std::string fileName) cons Line 1366  Data::saveVTK(std::string fileName) cons
1366  Data&  Data&
1367  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1368  {  {
1369  #if defined DOPROF    if (isProtected()) {
1370    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1371  #endif    }
1372    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1373    return (*this);    return (*this);
1374  }  }
# Line 1324  Data::operator+=(const Data& right) Line 1376  Data::operator+=(const Data& right)
1376  Data&  Data&
1377  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1378  {  {
1379  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1380    profData->binary++;    binaryOp(tmp,plus<double>());
 #endif  
   binaryOp(right,plus<double>());  
1381    return (*this);    return (*this);
1382  }  }
1383    
1384  Data&  Data&
1385  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1386  {  {
1387  #if defined DOPROF    if (isProtected()) {
1388    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1389  #endif    }
1390    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1391    return (*this);    return (*this);
1392  }  }
# Line 1344  Data::operator-=(const Data& right) Line 1394  Data::operator-=(const Data& right)
1394  Data&  Data&
1395  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1396  {  {
1397  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1398    profData->binary++;    binaryOp(tmp,minus<double>());
 #endif  
   binaryOp(right,minus<double>());  
1399    return (*this);    return (*this);
1400  }  }
1401    
1402  Data&  Data&
1403  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1404  {  {
1405  #if defined DOPROF    if (isProtected()) {
1406    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1407  #endif    }
1408    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1409    return (*this);    return (*this);
1410  }  }
# Line 1364  Data::operator*=(const Data& right) Line 1412  Data::operator*=(const Data& right)
1412  Data&  Data&
1413  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1414  {  {
1415  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1416    profData->binary++;    binaryOp(tmp,multiplies<double>());
 #endif  
   binaryOp(right,multiplies<double>());  
1417    return (*this);    return (*this);
1418  }  }
1419    
1420  Data&  Data&
1421  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1422  {  {
1423  #if defined DOPROF    if (isProtected()) {
1424    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1425  #endif    }
1426    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1427    return (*this);    return (*this);
1428  }  }
# Line 1384  Data::operator/=(const Data& right) Line 1430  Data::operator/=(const Data& right)
1430  Data&  Data&
1431  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1432  {  {
1433  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1434    profData->binary++;    binaryOp(tmp,divides<double>());
 #endif  
   binaryOp(right,divides<double>());  
1435    return (*this);    return (*this);
1436  }  }
1437    
1438  Data  Data
1439  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1440  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1441    Data left_d(left,*this);    Data left_d(left,*this);
1442    return left_d.powD(*this);    return left_d.powD(*this);
1443  }  }
# Line 1404  Data::rpowO(const boost::python::object& Line 1445  Data::rpowO(const boost::python::object&
1445  Data  Data
1446  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1447  {  {
1448  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1449    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1450  }  }
1451    
1452  Data  Data
1453  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1454  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1455    Data result;    Data result;
1456    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1457    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1458         result.binaryOp(*this,escript::rpow);
1459      } else {
1460         result.copy(*this);
1461         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1462      }
1463    return result;    return result;
1464  }  }
1465    
1466    
1467  //  //
1468  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1469  Data  Data
# Line 1433  escript::operator+(const Data& left, con Line 1472  escript::operator+(const Data& left, con
1472    Data result;    Data result;
1473    //    //
1474    // perform a deep copy    // perform a deep copy
1475    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1476    result+=right;       result.copy(right);
1477         result+=left;
1478      } else {
1479         result.copy(left);
1480         result+=right;
1481      }
1482    return result;    return result;
1483  }  }
1484    
# Line 1446  escript::operator-(const Data& left, con Line 1490  escript::operator-(const Data& left, con
1490    Data result;    Data result;
1491    //    //
1492    // perform a deep copy    // perform a deep copy
1493    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1494    result-=right;       result=right.neg();
1495         result+=left;
1496      } else {
1497         result.copy(left);
1498         result-=right;
1499      }
1500    return result;    return result;
1501  }  }
1502    
# Line 1459  escript::operator*(const Data& left, con Line 1508  escript::operator*(const Data& left, con
1508    Data result;    Data result;
1509    //    //
1510    // perform a deep copy    // perform a deep copy
1511    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1512    result*=right;       result.copy(right);
1513         result*=left;
1514      } else {
1515         result.copy(left);
1516         result*=right;
1517      }
1518    return result;    return result;
1519  }  }
1520    
# Line 1472  escript::operator/(const Data& left, con Line 1526  escript::operator/(const Data& left, con
1526    Data result;    Data result;
1527    //    //
1528    // perform a deep copy    // perform a deep copy
1529    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1530    result/=right;       result=right.oneOver();
1531         result*=left;
1532      } else {
1533         result.copy(left);
1534         result/=right;
1535      }
1536    return result;    return result;
1537  }  }
1538    
# Line 1482  escript::operator/(const Data& left, con Line 1541  escript::operator/(const Data& left, con
1541  Data  Data
1542  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1543  {  {
1544    //    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;  
1545  }  }
1546    
1547  //  //
# Line 1498  escript::operator+(const Data& left, con Line 1549  escript::operator+(const Data& left, con
1549  Data  Data
1550  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1551  {  {
1552    //    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;  
1553  }  }
1554    
1555  //  //
# Line 1514  escript::operator-(const Data& left, con Line 1557  escript::operator-(const Data& left, con
1557  Data  Data
1558  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1559  {  {
1560    //    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;  
1561  }  }
1562    
1563  //  //
# Line 1530  escript::operator*(const Data& left, con Line 1565  escript::operator*(const Data& left, con
1565  Data  Data
1566  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1567  {  {
1568    //    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;  
1569  }  }
1570    
1571  //  //
# Line 1546  escript::operator/(const Data& left, con Line 1573  escript::operator/(const Data& left, con
1573  Data  Data
1574  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1575  {  {
1576    //    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;  
1577  }  }
1578    
1579  //  //
# Line 1559  escript::operator+(const boost::python:: Line 1581  escript::operator+(const boost::python::
1581  Data  Data
1582  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1583  {  {
1584    //    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;  
1585  }  }
1586    
1587  //  //
# Line 1572  escript::operator-(const boost::python:: Line 1589  escript::operator-(const boost::python::
1589  Data  Data
1590  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1591  {  {
1592    //    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;  
1593  }  }
1594    
1595  //  //
# Line 1585  escript::operator*(const boost::python:: Line 1597  escript::operator*(const boost::python::
1597  Data  Data
1598  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1599  {  {
1600    //    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;  
1601  }  }
1602    
1603  //  //
# Line 1638  escript::operator/(const boost::python:: Line 1645  escript::operator/(const boost::python::
1645  //  return ret;  //  return ret;
1646  //}  //}
1647    
1648    /* TODO */
1649    /* global reduction */
1650  Data  Data
1651  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1652  {  {
# Line 1652  Data::getItem(const boost::python::objec Line 1661  Data::getItem(const boost::python::objec
1661    return getSlice(slice_region);    return getSlice(slice_region);
1662  }  }
1663    
1664    /* TODO */
1665    /* global reduction */
1666  Data  Data
1667  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1668  {  {
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
1669    return Data(*this,region);    return Data(*this,region);
1670  }  }
1671    
1672    /* TODO */
1673    /* global reduction */
1674  void  void
1675  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1676                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1669  Data::setItemO(const boost::python::obje Line 1679  Data::setItemO(const boost::python::obje
1679    setItemD(key,tempData);    setItemD(key,tempData);
1680  }  }
1681    
1682    /* TODO */
1683    /* global reduction */
1684  void  void
1685  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1686                 const Data& value)                 const Data& value)
# Line 1686  Data::setItemD(const boost::python::obje Line 1698  Data::setItemD(const boost::python::obje
1698    }    }
1699  }  }
1700    
1701    /* TODO */
1702    /* global reduction */
1703  void  void
1704  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1705                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1706  {  {
1707  #if defined DOPROF    if (isProtected()) {
1708    profData->slicing++;          throw DataException("Error - attempt to update protected Data object.");
1709  #endif    }
1710    Data tempValue(value);    Data tempValue(value);
1711    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1712    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1728  Data::typeMatchRight(const Data& right) Line 1742  Data::typeMatchRight(const Data& right)
1742  }  }
1743    
1744  void  void
1745    Data::setTaggedValueByName(std::string name,
1746                               const boost::python::object& value)
1747    {
1748         if (getFunctionSpace().getDomain().isValidTagName(name)) {
1749            int tagKey=getFunctionSpace().getDomain().getTag(name);
1750            setTaggedValue(tagKey,value);
1751         }
1752    }
1753    void
1754  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1755                       const boost::python::object& value)                       const boost::python::object& value)
1756  {  {
1757      if (isProtected()) {
1758            throw DataException("Error - attempt to update protected Data object.");
1759      }
1760    //    //
1761    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1762    tag();    tag();
# Line 1752  void Line 1778  void
1778  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1779                              const DataArrayView& value)                              const DataArrayView& value)
1780  {  {
1781      if (isProtected()) {
1782            throw DataException("Error - attempt to update protected Data object.");
1783      }
1784    //    //
1785    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1786    tag();    tag();
# Line 1772  Data::getTagNumber(int dpno) Line 1801  Data::getTagNumber(int dpno)
1801  }  }
1802    
1803  void  void
 Data::setRefValue(int ref,  
                   const boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 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);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
 void  
1804  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
1805  {  {
1806    cout << "Archiving Data object to: " << fileName << endl;    cout << "Archiving Data object to: " << fileName << endl;
# Line 1887  Data::archiveData(const std::string file Line 1839  Data::archiveData(const std::string file
1839    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1840    int dataLength = getLength();    int dataLength = getLength();
1841    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
1842    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
1843    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1844      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1845    }    }
1846    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
1847    if (isTagged()) {    if (isTagged()) {
1848      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1849        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 2045  Data::extractData(const std::string file Line 1997  Data::extractData(const std::string file
1997        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
1998      }      }
1999    }    }
2000    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2001    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2002      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2003    }    }
2004    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2005    if (dataType==2) {    if (dataType==2) {
2006      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2007        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 2098  Data::extractData(const std::string file Line 2050  Data::extractData(const std::string file
2050      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2051    }    }
2052    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2053      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2054        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2055      }      }
2056    }    }
# Line 2188  ostream& escript::operator<<(ostream& o, Line 2140  ostream& escript::operator<<(ostream& o,
2140    o << data.toString();    o << data.toString();
2141    return o;    return o;
2142  }  }
2143    
2144    Data
2145    escript::C_GeneralTensorProduct(Data& arg_0,
2146                         Data& arg_1,
2147                         int axis_offset,
2148                         int transpose)
2149    {
2150      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2151      // SM is the product of the last axis_offset entries in arg_0.getShape().
2152    
2153    
2154      // Interpolate if necessary and find an appropriate function space
2155      Data arg_0_Z, arg_1_Z;
2156      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2157        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2158          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2159          arg_1_Z = Data(arg_1);
2160        }
2161        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2162          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2163          arg_0_Z =Data(arg_0);
2164        }
2165        else {
2166          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2167        }
2168      } else {
2169          arg_0_Z = Data(arg_0);
2170          arg_1_Z = Data(arg_1);
2171      }
2172      // Get rank and shape of inputs
2173      int rank0 = arg_0_Z.getDataPointRank();
2174      int rank1 = arg_1_Z.getDataPointRank();
2175      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2176      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2177    
2178      // Prepare for the loops of the product and verify compatibility of shapes
2179      int start0=0, start1=0;
2180      if (transpose == 0)       {}
2181      else if (transpose == 1)  { start0 = axis_offset; }
2182      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2183      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2184    
2185      // Adjust the shapes for transpose
2186      DataArrayView::ShapeType tmpShape0;
2187      DataArrayView::ShapeType tmpShape1;
2188      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2189      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2190    
2191    #if 0
2192      // For debugging: show shape after transpose
2193      char tmp[100];
2194      std::string shapeStr;
2195      shapeStr = "(";
2196      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2197      shapeStr += ")";
2198      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2199      shapeStr = "(";
2200      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2201      shapeStr += ")";
2202      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2203    #endif
2204    
2205      // Prepare for the loops of the product
2206      int SL=1, SM=1, SR=1;
2207      for (int i=0; i<rank0-axis_offset; i++)   {
2208        SL *= tmpShape0[i];
2209      }
2210      for (int i=rank0-axis_offset; i<rank0; i++)   {
2211        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2212          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2213        }
2214        SM *= tmpShape0[i];
2215      }
2216      for (int i=axis_offset; i<rank1; i++)     {
2217        SR *= tmpShape1[i];
2218      }
2219    
2220      // Define the shape of the output
2221      DataArrayView::ShapeType shape2;
2222      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2223      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2224    
2225      // Declare output Data object
2226      Data res;
2227    
2228      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2229        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2230        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2231        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2232        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2233        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2234      }
2235      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2236    
2237        // Prepare the DataConstant input
2238        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2239        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2240    
2241        // Borrow DataTagged input from Data object
2242        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2243        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2244    
2245        // Prepare a DataTagged output 2
2246        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2247        res.tag();
2248        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2249        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2250    
2251        // Prepare offset into DataConstant
2252        int offset_0 = tmp_0->getPointOffset(0,0);
2253        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2254        // Get the views
2255        DataArrayView view_1 = tmp_1->getDefaultValue();
2256        DataArrayView view_2 = tmp_2->getDefaultValue();
2257        // Get the pointers to the actual data
2258        double *ptr_1 = &((view_1.getData())[0]);
2259        double *ptr_2 = &((view_2.getData())[0]);
2260        // Compute an MVP for the default
2261        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2262        // Compute an MVP for each tag
2263        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2264        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2265        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2266          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2267          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2268          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2269          double *ptr_1 = &view_1.getData(0);
2270          double *ptr_2 = &view_2.getData(0);
2271          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2272        }
2273    
2274      }
2275      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2276    
2277        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2278        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2279        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2280        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2281        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2282        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2283        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2284        int sampleNo_1,dataPointNo_1;
2285        int numSamples_1 = arg_1_Z.getNumSamples();
2286        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2287        int offset_0 = tmp_0->getPointOffset(0,0);
2288        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2289        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2290          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2291            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2292            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2293            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2294            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2295            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2296            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2297          }
2298        }
2299    
2300      }
2301      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2302    
2303        // Borrow DataTagged input from Data object
2304        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2305        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2306    
2307        // Prepare the DataConstant input
2308        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2309        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2310    
2311        // Prepare a DataTagged output 2
2312        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2313        res.tag();
2314        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2315        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2316    
2317        // Prepare offset into DataConstant
2318        int offset_1 = tmp_1->getPointOffset(0,0);
2319        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2320        // Get the views
2321        DataArrayView view_0 = tmp_0->getDefaultValue();
2322        DataArrayView view_2 = tmp_2->getDefaultValue();
2323        // Get the pointers to the actual data
2324        double *ptr_0 = &((view_0.getData())[0]);
2325        double *ptr_2 = &((view_2.getData())[0]);
2326        // Compute an MVP for the default
2327        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2328        // Compute an MVP for each tag
2329        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2330        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2331        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2332          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2333          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2334          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2335          double *ptr_0 = &view_0.getData(0);
2336          double *ptr_2 = &view_2.getData(0);
2337          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2338        }
2339    
2340      }
2341      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2342    
2343        // Borrow DataTagged input from Data object
2344        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2345        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2346    
2347        // Borrow DataTagged input from Data object
2348        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2349        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2350    
2351        // Prepare a DataTagged output 2
2352        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2353        res.tag();  // DataTagged output
2354        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2355        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2356    
2357        // Get the views
2358        DataArrayView view_0 = tmp_0->getDefaultValue();
2359        DataArrayView view_1 = tmp_1->getDefaultValue();
2360        DataArrayView view_2 = tmp_2->getDefaultValue();
2361        // Get the pointers to the actual data
2362        double *ptr_0 = &((view_0.getData())[0]);
2363        double *ptr_1 = &((view_1.getData())[0]);
2364        double *ptr_2 = &((view_2.getData())[0]);
2365        // Compute an MVP for the default
2366        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2367        // Merge the tags
2368        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2369        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2370        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2371        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2372          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2373        }
2374        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2375          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2376        }
2377        // Compute an MVP for each tag
2378        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2379        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2380          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2381          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2382          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2383          double *ptr_0 = &view_0.getData(0);
2384          double *ptr_1 = &view_1.getData(0);
2385          double *ptr_2 = &view_2.getData(0);
2386          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2387        }
2388    
2389      }
2390      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2391    
2392        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2393        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2394        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2395        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2396        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2397        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2398        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2399        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2400        int sampleNo_0,dataPointNo_0;
2401        int numSamples_0 = arg_0_Z.getNumSamples();
2402        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2403        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2404        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2405          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2406          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2407          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2408            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2409            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2410            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2411            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2412            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2413          }
2414        }
2415    
2416      }
2417      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2418    
2419        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2420        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2421        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2422        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2423        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2424        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2425        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2426        int sampleNo_0,dataPointNo_0;
2427        int numSamples_0 = arg_0_Z.getNumSamples();
2428        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2429        int offset_1 = tmp_1->getPointOffset(0,0);
2430        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2431        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2432          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2433            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2434            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2435            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2436            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2437            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2438            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2439          }
2440        }
2441    
2442    
2443      }
2444      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2445    
2446        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2447        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2448        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2449        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2450        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2451        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2452        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2453        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2454        int sampleNo_0,dataPointNo_0;
2455        int numSamples_0 = arg_0_Z.getNumSamples();
2456        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2457        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2458        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2459          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2460          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2461          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2462            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2463            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2464            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2465            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2466            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2467          }
2468        }
2469    
2470      }
2471      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2472    
2473        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2474        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2475        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2476        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2477        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2478        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2479        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2480        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2481        int sampleNo_0,dataPointNo_0;
2482        int numSamples_0 = arg_0_Z.getNumSamples();
2483        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2484        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2485        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2486          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2487            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2488            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2489            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2490            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2491            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2492            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2493            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2494          }
2495        }
2496    
2497      }
2498      else {
2499        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2500      }
2501    
2502      return res;
2503    }
2504    
2505    DataAbstract*
2506    Data::borrowData() const
2507    {
2508      return m_data.get();
2509    }
2510    
2511    /* Member functions specific to the MPI implementation */
2512    
2513    void
2514    Data::print()
2515    {
2516      int i,j;
2517      
2518      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2519      for( i=0; i<getNumSamples(); i++ )
2520      {
2521        printf( "[%6d]", i );
2522        for( j=0; j<getNumDataPointsPerSample(); j++ )
2523          printf( "\t%10.7g", (getSampleData(i))[j] );
2524        printf( "\n" );
2525      }
2526    }
2527    
2528    int
2529    Data::get_MPISize() const
2530    {
2531        int error, size;
2532    #ifdef PASO_MPI
2533        error = MPI_Comm_size( get_MPIComm(), &size );
2534    #else
2535        size = 1;
2536    #endif
2537        return size;
2538    }
2539    
2540    int
2541    Data::get_MPIRank() const
2542    {
2543        int error, rank;
2544    #ifdef PASO_MPI
2545        error = MPI_Comm_rank( get_MPIComm(), &rank );
2546    #else
2547        rank = 0;
2548    #endif
2549        return rank;
2550    }
2551    
2552    MPI_Comm
2553    Data::get_MPIComm() const
2554    {
2555    #ifdef PASO_MPI
2556        return MPI_COMM_WORLD;
2557    #else
2558        return -1;
2559    #endif
2560    }
2561    

Legend:
Removed from v.711  
changed lines
  Added in v.1092

  ViewVC Help
Powered by ViewVC 1.1.26