/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
 // $Id$  
1    
2  /*  /* $Id$ */
3   ************************************************************  
4   *          Copyright 2006 by ACcESS MNRF                   *  /*******************************************************
5   *                                                          *   *
6   *              http://www.access.edu.au                    *   *           Copyright 2003-2007 by ACceSS MNRF
7   *       Primary Business: Queensland, Australia            *   *       Copyright 2007 by University of Queensland
8   *  Licensed under the Open Software License version 3.0    *   *
9   *     http://www.opensource.org/licenses/osl-3.0.php       *   *                http://esscc.uq.edu.au
10   *                                                          *   *        Primary Business: Queensland, Australia
11   ************************************************************   *  Licensed under the Open Software License version 3.0
12  */   *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  #include "Data.h"  #include "Data.h"
17    
18  #include "DataExpanded.h"  #include "DataExpanded.h"
# Line 19  Line 21 
21  #include "DataEmpty.h"  #include "DataEmpty.h"
22  #include "DataArray.h"  #include "DataArray.h"
23  #include "DataArrayView.h"  #include "DataArrayView.h"
 #include "DataProf.h"  
24  #include "FunctionSpaceFactory.h"  #include "FunctionSpaceFactory.h"
25  #include "AbstractContinuousDomain.h"  #include "AbstractContinuousDomain.h"
26  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
27    extern "C" {
28    #include "escript/blocktimer.h"
29    }
30    
31  #include <fstream>  #include <fstream>
32  #include <algorithm>  #include <algorithm>
# Line 38  using namespace boost::python; Line 42  using namespace boost::python;
42  using namespace boost;  using namespace boost;
43  using namespace escript;  using namespace escript;
44    
 #if defined DOPROF  
 //  
 // global table of profiling data for all Data objects  
 DataProf dataProfTable;  
 #endif  
   
45  Data::Data()  Data::Data()
46  {  {
47    //    //
# Line 52  Data::Data() Line 50  Data::Data()
50    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
51    m_data=temp_data;    m_data=temp_data;
52    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
53  }  }
54    
55  Data::Data(double value,  Data::Data(double value,
# Line 70  Data::Data(double value, Line 64  Data::Data(double value,
64    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
65    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
66    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
67  }  }
68    
69  Data::Data(double value,  Data::Data(double value,
# Line 82  Data::Data(double value, Line 72  Data::Data(double value,
72             bool expanded)             bool expanded)
73  {  {
74    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
   pair<int,int> dataShape=what.getDataShape();  
75    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
76    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
77  }  }
78    
79  Data::Data(const Data& inData)  Data::Data(const Data& inData)
80  {  {
81    m_data=inData.m_data;    m_data=inData.m_data;
82    m_protected=inData.isProtected();    m_protected=inData.isProtected();
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
83  }  }
84    
85  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 110  Data::Data(const Data& inData, Line 91  Data::Data(const Data& inData,
91    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
92    m_data=temp_data;    m_data=temp_data;
93    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
94  }  }
95    
96  Data::Data(const Data& inData,  Data::Data(const Data& inData,
97             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
98  {  {
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
99    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
100      m_data=inData.m_data;      m_data=inData.m_data;
101    } else {    } else {
     #if defined DOPROF  
     profData->interpolate++;  
     #endif  
102      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
103      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
104      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 158  Data::Data(const DataTagged::TagListType Line 128  Data::Data(const DataTagged::TagListType
128    if (expanded) {    if (expanded) {
129      expand();      expand();
130    }    }
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
131  }  }
132    
133  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 170  Data::Data(const numeric::array& value, Line 136  Data::Data(const numeric::array& value,
136  {  {
137    initialise(value,what,expanded);    initialise(value,what,expanded);
138    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
139  }  }
140    
141  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 182  Data::Data(const DataArrayView& value, Line 144  Data::Data(const DataArrayView& value,
144  {  {
145    initialise(value,what,expanded);    initialise(value,what,expanded);
146    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
147  }  }
148    
149  Data::Data(const object& value,  Data::Data(const object& value,
# Line 195  Data::Data(const object& value, Line 153  Data::Data(const object& value,
153    numeric::array asNumArray(value);    numeric::array asNumArray(value);
154    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
155    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
156  }  }
157    
158  Data::Data(const object& value,  Data::Data(const object& value,
# Line 220  Data::Data(const object& value, Line 174  Data::Data(const object& value,
174      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
175    }    }
176    m_protected=false;    m_protected=false;
 #if defined DOPROF  
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
177  }  }
178    
179  Data::~Data()  Data::~Data()
# Line 266  Data::getShapeTuple() const Line 216  Data::getShapeTuple() const
216          throw DataException("Error - illegal Data rank.");          throw DataException("Error - illegal Data rank.");
217    }    }
218  }  }
   
219  void  void
220  Data::copy(const Data& other)  Data::copy(const Data& other)
221  {  {
# Line 319  Data::copy(const Data& other) Line 268  Data::copy(const Data& other)
268    throw DataException("Error - Copy not implemented for this Data type.");    throw DataException("Error - Copy not implemented for this Data type.");
269  }  }
270    
271    
272    void
273    Data::setToZero()
274    {
275      {
276        DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
277        if (temp!=0) {
278           temp->setToZero();
279           return;
280        }
281      }
282      {
283        DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
284        if (temp!=0) {
285          temp->setToZero();
286          return;
287        }
288      }
289      {
290        DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
291        if (temp!=0) {
292          temp->setToZero();
293          return;
294        }
295      }
296      throw DataException("Error - Data can not be set to zero.");
297    }
298    
299  void  void
300  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
301                     const Data& mask)                     const Data& mask)
# Line 350  Data::isTagged() const Line 327  Data::isTagged() const
327    return (temp!=0);    return (temp!=0);
328  }  }
329    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
330  bool  bool
331  Data::isEmpty() const  Data::isEmpty() const
332  {  {
# Line 367  Data::isConstant() const Line 342  Data::isConstant() const
342  }  }
343    
344  void  void
345  Data::setProtection()  Data::setProtection()
346  {  {
347     m_protected=true;     m_protected=true;
348  }  }
349    
350  bool  bool
351  Data::isProtected() const  Data::isProtected() const
352  {  {
353     return m_protected;     return m_protected;
354  }  }
355    
# Line 422  Data::tag() Line 397  Data::tag()
397    }    }
398  }  }
399    
400  void  Data
401  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
402  {  {
403    m_data->reshapeDataPoint(shape);    return escript::unaryOp(*this,bind1st(divides<double>(),1.));
404  }  }
405    
406  Data  Data
407  Data::wherePositive() const  Data::wherePositive() const
408  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
409    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
410  }  }
411    
412  Data  Data
413  Data::whereNegative() const  Data::whereNegative() const
414  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
415    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
416  }  }
417    
418  Data  Data
419  Data::whereNonNegative() const  Data::whereNonNegative() const
420  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
421    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
422  }  }
423    
424  Data  Data
425  Data::whereNonPositive() const  Data::whereNonPositive() const
426  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
427    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
428  }  }
429    
430  Data  Data
431  Data::whereZero(double tol) const  Data::whereZero(double tol) const
432  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
433    Data dataAbs=abs();    Data dataAbs=abs();
434    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
435  }  }
# Line 477  Data::whereZero(double tol) const Line 437  Data::whereZero(double tol) const
437  Data  Data
438  Data::whereNonZero(double tol) const  Data::whereNonZero(double tol) const
439  {  {
 #if defined DOPROF  
   profData->where++;  
 #endif  
440    Data dataAbs=abs();    Data dataAbs=abs();
441    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));    return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
442  }  }
# Line 487  Data::whereNonZero(double tol) const Line 444  Data::whereNonZero(double tol) const
444  Data  Data
445  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
446  {  {
 #if defined DOPROF  
   profData->interpolate++;  
 #endif  
447    return Data(*this,functionspace);    return Data(*this,functionspace);
448  }  }
449    
# Line 511  Data::probeInterpolation(const FunctionS Line 465  Data::probeInterpolation(const FunctionS
465  Data  Data
466  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
467  {  {
468  #if defined DOPROF    double blocktimer_start = blocktimer_time();
   profData->grad++;  
 #endif  
469    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
470      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
471    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
472    grad_shape.push_back(functionspace.getDim());    grad_shape.push_back(functionspace.getDim());
473    Data out(0.0,grad_shape,functionspace,true);    Data out(0.0,grad_shape,functionspace,true);
474    getDomain().setToGradient(out,*this);    getDomain().setToGradient(out,*this);
475      blocktimer_increment("grad()", blocktimer_start);
476    return out;    return out;
477  }  }
478    
# Line 547  Data::getDataPointShape() const Line 500  Data::getDataPointShape() const
500    return getPointDataView().getShape();    return getPointDataView().getShape();
501  }  }
502    
 void  
 Data::fillFromNumArray(const boost::python::numeric::array num_array)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
   //  
   // check rank  
   if (num_array.getrank()<getDataPointRank())  
       throw DataException("Rank of numarray does not match Data object rank");  
   
   //  
   // check shape of num_array  
   for (int i=0; i<getDataPointRank(); i++) {  
     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])  
        throw DataException("Shape of numarray does not match Data object rank");  
   }  
503    
   //  
   // make sure data is expanded:  
   if (!isExpanded()) {  
     expand();  
   }  
   
   //  
   // and copy over  
   m_data->copyAll(num_array);  
 }  
504    
505  const  const
506  boost::python::numeric::array  boost::python::numeric::array
507  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
508  {  {
509    //    size_t length=0;
510    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
511    //    //
512    // determine the rank and shape of each data point    // determine the rank and shape of each data point
513    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 596  Data::convertToNumArray() Line 518  Data::convertToNumArray()
518    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
519    
520    //    //
521    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
522    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
523    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
524    // 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]);  
   }  
525    
526    //    //
527    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
528      if (arrayRank==0) {
529        numArray.resize(1);
530      }
531    if (arrayRank==1) {    if (arrayRank==1) {
532      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
533    }    }
# Line 623  Data::convertToNumArray() Line 540  Data::convertToNumArray()
540    if (arrayRank==4) {    if (arrayRank==4) {
541      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
542    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
543    
544    //    if (getNumDataPointsPerSample()>0) {
545    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
546    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
547    int dataPoint = 0;         //
548    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
549      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
550        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
551        if (dataPointRank==0) {         }
         numArray[dataPoint]=dataPointView();  
       }  
       if (dataPointRank==1) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           numArray[dataPoint][i]=dataPointView(i);  
         }  
       }  
       if (dataPointRank==2) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             numArray[dataPoint][i][j] = dataPointView(i,j);  
           }  
         }  
       }  
       if (dataPointRank==3) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             for (int k=0; k<dataPointShape[2]; k++) {  
               numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
             }  
           }  
         }  
       }  
       if (dataPointRank==4) {  
         for (int i=0; i<dataPointShape[0]; i++) {  
           for (int j=0; j<dataPointShape[1]; j++) {  
             for (int k=0; k<dataPointShape[2]; k++) {  
               for (int l=0; l<dataPointShape[3]; l++) {  
                 numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
               }  
             }  
           }  
         }  
       }  
       dataPoint++;  
     }  
   }  
552    
553           //
554           // Check a valid data point number has been supplied
555           if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
556               throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
557           }
558           // TODO: global error handling
559           // create a view of the data if it is stored locally
560           DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
561    
562           switch( dataPointRank ){
563                case 0 :
564                    numArray[0] = dataPointView();
565                    break;
566                case 1 :
567                    for( i=0; i<dataPointShape[0]; i++ )
568                        numArray[i]=dataPointView(i);
569                    break;
570                case 2 :
571                    for( i=0; i<dataPointShape[0]; i++ )
572                        for( j=0; j<dataPointShape[1]; j++)
573                            numArray[make_tuple(i,j)]=dataPointView(i,j);
574                    break;
575                case 3 :
576                    for( i=0; i<dataPointShape[0]; i++ )
577                        for( j=0; j<dataPointShape[1]; j++ )
578                            for( k=0; k<dataPointShape[2]; k++)
579                                numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
580                    break;
581                case 4 :
582                    for( i=0; i<dataPointShape[0]; i++ )
583                        for( j=0; j<dataPointShape[1]; j++ )
584                            for( k=0; k<dataPointShape[2]; k++ )
585                                for( l=0; l<dataPointShape[3]; l++)
586                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
587                    break;
588        }
589      }
590    //    //
591    // return the loaded array    // return the array
592    return numArray;    return numArray;
 }  
593    
594  const  }
595  boost::python::numeric::array  void
596  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
597  {  {
598    //      // this will throw if the value cannot be represented
599    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
600    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
601    
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
602    
603    //  }
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
604    
605    void
606    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
607    {
608      if (isProtected()) {
609            throw DataException("Error - attempt to update protected Data object.");
610      }
611    //    //
612    // create the numeric array to be returned    // check rank
613    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
614          throw DataException("Rank of numarray does not match Data object rank");
615    
616    //    //
617    // the rank of the returned numeric array will be the rank of    // check shape of num_array
618    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
619    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
620    // 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]);  
621    }    }
   
622    //    //
623    // resize the numeric array to the shape just calculated    // make sure data is expanded:
624    if (arrayRank==1) {    if (!isExpanded()) {
625      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]);  
626    }    }
627    if (arrayRank==5) {    if (getNumDataPointsPerSample()>0) {
628      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);         int sampleNo = dataPointNo/getNumDataPointsPerSample();
629           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
630           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
631      } else {
632           m_data->copyToDataPoint(-1, 0,num_array);
633    }    }
634    }
635    
636    //  void
637    // loop through each data point in turn, loading the values for that data point  Data::setValueOfDataPoint(int dataPointNo, const double value)
638    // into the numeric array.  {
639    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    if (isProtected()) {
640      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);          throw DataException("Error - attempt to update protected Data object.");
     if (dataPointRank==0) {  
       numArray[dataPoint]=dataPointView();  
     }  
     if (dataPointRank==1) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         numArray[dataPoint][i]=dataPointView(i);  
       }  
     }  
     if (dataPointRank==2) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           numArray[dataPoint][i][j] = dataPointView(i,j);  
         }  
       }  
     }  
     if (dataPointRank==3) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
           }  
         }  
       }  
     }  
     if (dataPointRank==4) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             for (int l=0; l<dataPointShape[3]; l++) {  
               numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
             }  
           }  
         }  
       }  
     }  
641    }    }
   
642    //    //
643    // return the loaded array    // make sure data is expanded:
644    return numArray;    if (!isExpanded()) {
645        expand();
646      }
647      if (getNumDataPointsPerSample()>0) {
648           int sampleNo = dataPointNo/getNumDataPointsPerSample();
649           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
650           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
651      } else {
652           m_data->copyToDataPoint(-1, 0,value);
653      }
654  }  }
655    
656  const  const
657  boost::python::numeric::array  boost::python::numeric::array
658  Data::convertToNumArrayFromDPNo(int procNo,  Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
                                 int sampleNo,  
                                                                 int dataPointNo)  
   
659  {  {
660      size_t length=0;    size_t length=0;
661      int i, j, k, l, pos;    int i, j, k, l, pos;
   
   //  
   // Check a valid sample number has been supplied  
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // Check a valid data point number has been supplied  
   if (dataPointNo >= getNumDataPointsPerSample()) {  
     throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");  
   }  
   
662    //    //
663    // determine the rank and shape of each data point    // determine the rank and shape of each data point
664    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 692  Data::convertToNumArrayFromDPNo(int proc
692      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
693    }    }
694    
695      // added for the MPI communication    // added for the MPI communication
696      length=1;    length=1;
697      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
698          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
699    
700    //    //
701    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
702    
703      // updated for the MPI case      // updated for the MPI case
704      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
705                 if (getNumDataPointsPerSample()>0) {
706                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
707                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
708                    //
709                    // Check a valid sample number has been supplied
710                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
711                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
712                    }
713    
714                    //
715                    // Check a valid data point number has been supplied
716                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
717                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
718                    }
719                    // TODO: global error handling
720          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
721          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
722            
723          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
724          pos=0;          pos=0;
725          switch( dataPointRank ){          switch( dataPointRank ){
726              case 0 :              case 0 :
727                  tmpData[0] = dataPointView();                  tmpData[0] = dataPointView();
728                  break;                  break;
729              case 1 :                      case 1 :
730                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
731                      tmpData[i]=dataPointView(i);                      tmpData[i]=dataPointView(i);
732                  break;                  break;
733              case 2 :                      case 2 :
734                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
735                      for( j=0; j<dataPointShape[1]; j++, pos++ )                      for( j=0; j<dataPointShape[1]; j++, pos++ )
736                          tmpData[pos]=dataPointView(i,j);                          tmpData[pos]=dataPointView(i,j);
737                  break;                  break;
738              case 3 :                      case 3 :
739                  for( i=0; i<dataPointShape[0]; i++ )                  for( i=0; i<dataPointShape[0]; i++ )
740                      for( j=0; j<dataPointShape[1]; j++ )                      for( j=0; j<dataPointShape[1]; j++ )
741                          for( k=0; k<dataPointShape[2]; k++, pos++ )                          for( k=0; k<dataPointShape[2]; k++, pos++ )
# Line 878  Data::convertToNumArrayFromDPNo(int proc Line 749  Data::convertToNumArrayFromDPNo(int proc
749                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=dataPointView(i,j,k,l);
750                  break;                  break;
751          }          }
752                }
753      }      }
754  #ifdef PASO_MPI          #ifdef PASO_MPI
755          // broadcast the data to all other processes          // broadcast the data to all other processes
756          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
757  #endif          #endif
758    
759      // unpack the data      // unpack the data
760      switch( dataPointRank ){      switch( dataPointRank ){
761          case 0 :          case 0 :
762              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
763              break;              break;
764          case 1 :                  case 1 :
765              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
766                  numArray[i]=tmpData[i];                  numArray[i]=tmpData[i];
767              break;              break;
768          case 2 :                  case 2 :
769              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
770                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
771                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
772              break;              break;
773          case 3 :                  case 3 :
774              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
775                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
776                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
777                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
778              break;              break;
779          case 4 :          case 4 :
780              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
781                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
782                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
783                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
784                              tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];                                  numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
785              break;              break;
786      }      }
787    
788      delete [] tmpData;        delete [] tmpData;
 /*  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
 */  
   
789    //    //
790    // return the loaded array    // return the loaded array
791    return numArray;    return numArray;
792  }  }
793    
794    
795    
796  boost::python::numeric::array  boost::python::numeric::array
797  Data::integrate() const  Data::integrate() const
798  {  {
799    int index;    int index;
800    int rank = getDataPointRank();    int rank = getDataPointRank();
801    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
802      int dataPointSize = getDataPointSize();
 #if defined DOPROF  
   profData->integrate++;  
 #endif  
803    
804    //    //
805    // calculate the integral values    // calculate the integral values
806    vector<double> integrals(getDataPointSize());    vector<double> integrals(dataPointSize);
807      vector<double> integrals_local(dataPointSize);
808    #ifdef PASO_MPI
809      AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
810      // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
811      double *tmp = new double[dataPointSize];
812      double *tmp_local = new double[dataPointSize];
813      for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
814      MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
815      for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
816      delete[] tmp;
817      delete[] tmp_local;
818    #else
819    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);    AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
820    #endif
821    
822    //    //
823    // create the numeric array to be returned    // create the numeric array to be returned
# Line 1031  Data::integrate() const Line 877  Data::integrate() const
877  Data  Data
878  Data::sin() const  Data::sin() const
879  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
880    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
881  }  }
882    
883  Data  Data
884  Data::cos() const  Data::cos() const
885  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
886    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
887  }  }
888    
889  Data  Data
890  Data::tan() const  Data::tan() const
891  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
892    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
893  }  }
894    
895  Data  Data
896  Data::asin() const  Data::asin() const
897  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
898    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
899  }  }
900    
901  Data  Data
902  Data::acos() const  Data::acos() const
903  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
904    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
905  }  }
906    
907    
908  Data  Data
909  Data::atan() const  Data::atan() const
910  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
911    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
912  }  }
913    
914  Data  Data
915  Data::sinh() const  Data::sinh() const
916  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
917    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
918  }  }
919    
920  Data  Data
921  Data::cosh() const  Data::cosh() const
922  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
923    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
924  }  }
925    
926  Data  Data
927  Data::tanh() const  Data::tanh() const
928  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
929    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
930  }  }
931    
932    
933  Data  Data
934  Data::asinh() const  Data::erf() const
935  {  {
936  #if defined DOPROF  #ifdef _WIN32
937    profData->unary++;    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
938    #else
939      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
940  #endif  #endif
941    }
942    
943    Data
944    Data::asinh() const
945    {
946    #ifdef _WIN32
947      return escript::unaryOp(*this,escript::asinh_substitute);
948    #else
949    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
950    #endif
951  }  }
952    
953  Data  Data
954  Data::acosh() const  Data::acosh() const
955  {  {
956  #if defined DOPROF  #ifdef _WIN32
957    profData->unary++;    return escript::unaryOp(*this,escript::acosh_substitute);
958  #endif  #else
959    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
960    #endif
961  }  }
962    
963  Data  Data
964  Data::atanh() const  Data::atanh() const
965  {  {
966  #if defined DOPROF  #ifdef _WIN32
967    profData->unary++;    return escript::unaryOp(*this,escript::atanh_substitute);
968  #endif  #else
969    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
970    #endif
971  }  }
972    
973  Data  Data
974  Data::log10() const  Data::log10() const
975  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
976    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
977  }  }
978    
979  Data  Data
980  Data::log() const  Data::log() const
981  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
982    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
983  }  }
984    
985  Data  Data
986  Data::sign() const  Data::sign() const
987  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
988    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
989  }  }
990    
991  Data  Data
992  Data::abs() const  Data::abs() const
993  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
994    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
995  }  }
996    
997  Data  Data
998  Data::neg() const  Data::neg() const
999  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1000    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1001  }  }
1002    
1003  Data  Data
1004  Data::pos() const  Data::pos() const
1005  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1006    Data result;    Data result;
1007    // perform a deep copy    // perform a deep copy
1008    result.copy(*this);    result.copy(*this);
# Line 1196  Data::pos() const Line 1012  Data::pos() const
1012  Data  Data
1013  Data::exp() const  Data::exp() const
1014  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1015    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1016  }  }
1017    
1018  Data  Data
1019  Data::sqrt() const  Data::sqrt() const
1020  {  {
 #if defined DOPROF  
   profData->unary++;  
 #endif  
1021    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1022  }  }
1023    
# Line 1215  double Line 1025  double
1025  Data::Lsup() const  Data::Lsup() const
1026  {  {
1027    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1028    //    //
1029    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1030    
# Line 1232  Data::Lsup() const Line 1039  Data::Lsup() const
1039  }  }
1040    
1041  double  double
 Data::Linf() const  
 {  
   double localValue, globalValue;  
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
   //  
   // set the initial absolute minimum value to max double  
   AbsMin abs_min_func;  
   localValue = algorithm(abs_min_func,numeric_limits<double>::max());  
   
 #ifdef PASO_MPI  
   MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );  
   return globalValue;  
 #else  
   return localValue;  
 #endif  
 }  
   
 double  
1042  Data::sup() const  Data::sup() const
1043  {  {
1044    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1045    //    //
1046    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1047    FMax fmax_func;    FMax fmax_func;
# Line 1274  double Line 1058  double
1058  Data::inf() const  Data::inf() const
1059  {  {
1060    double localValue, globalValue;    double localValue, globalValue;
 #if defined DOPROF  
   profData->reduction1++;  
 #endif  
1061    //    //
1062    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1063    FMin fmin_func;    FMin fmin_func;
# Line 1294  Data::inf() const Line 1075  Data::inf() const
1075  Data  Data
1076  Data::maxval() const  Data::maxval() const
1077  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1078    //    //
1079    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1080    FMax fmax_func;    FMax fmax_func;
# Line 1306  Data::maxval() const Line 1084  Data::maxval() const
1084  Data  Data
1085  Data::minval() const  Data::minval() const
1086  {  {
 #if defined DOPROF  
   profData->reduction2++;  
 #endif  
1087    //    //
1088    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1089    FMin fmin_func;    FMin fmin_func;
# Line 1316  Data::minval() const Line 1091  Data::minval() const
1091  }  }
1092    
1093  Data  Data
1094  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1095  {  {
1096  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1097    profData->reduction2++;       DataArrayView::ShapeType s=getDataPointShape();
1098  #endif       DataArrayView::ShapeType ev_shape;
1099    Trace trace_func;       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1100    return dp_algorithm(trace_func,0);       // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1101         int rank=getDataPointRank();
1102         if (rank<2) {
1103            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1104         }
1105         if (axis0<0 || axis0>rank-1) {
1106            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1107         }
1108         if (axis1<0 || axis1>rank-1) {
1109             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1110         }
1111         if (axis0 == axis1) {
1112             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1113         }
1114         if (axis0 > axis1) {
1115             axis0_tmp=axis1;
1116             axis1_tmp=axis0;
1117         } else {
1118             axis0_tmp=axis0;
1119             axis1_tmp=axis1;
1120         }
1121         for (int i=0; i<rank; i++) {
1122           if (i == axis0_tmp) {
1123              ev_shape.push_back(s[axis1_tmp]);
1124           } else if (i == axis1_tmp) {
1125              ev_shape.push_back(s[axis0_tmp]);
1126           } else {
1127              ev_shape.push_back(s[i]);
1128           }
1129         }
1130         Data ev(0.,ev_shape,getFunctionSpace());
1131         ev.typeMatchRight(*this);
1132         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1133         return ev;
1134    
1135  }  }
1136    
1137  Data  Data
1138  Data::symmetric() const  Data::symmetric() const
1139  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1140       // check input       // check input
1141       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1142       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1143          if(s[0] != s[1])          if(s[0] != s[1])
1144             throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");             throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1145       }       }
1146       else if (getDataPointRank()==4) {       else if (getDataPointRank()==4) {
# Line 1353  Data::symmetric() const Line 1159  Data::symmetric() const
1159  Data  Data
1160  Data::nonsymmetric() const  Data::nonsymmetric() const
1161  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1162       // check input       // check input
1163       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1164       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1165          if(s[0] != s[1])          if(s[0] != s[1])
1166             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");             throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1167          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1168          ev_shape.push_back(s[0]);          ev_shape.push_back(s[0]);
# Line 1388  Data::nonsymmetric() const Line 1191  Data::nonsymmetric() const
1191  }  }
1192    
1193  Data  Data
1194  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1195  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1196       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1197       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
1198          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1199          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1200          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1201          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1202          return ev;          return ev;
1203       }       }
1204       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1213  Data::matrixtrace(int axis_offset) const
1213          }          }
1214          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1215          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1216          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1217          return ev;          return ev;
1218       }       }
1219       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1232  Data::matrixtrace(int axis_offset) const
1232      }      }
1233          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1234          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1235      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1236          return ev;          return ev;
1237       }       }
1238       else {       else {
1239          throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");          throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1240       }       }
1241  }  }
1242    
1243  Data  Data
1244  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1245  {  {
 #if defined DOPROF  
      profData->reduction2++;  
 #endif  
1246       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1247       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1248       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
# Line 1467  Data::transpose(int axis_offset) const Line 1264  Data::transpose(int axis_offset) const
1264  Data  Data
1265  Data::eigenvalues() const  Data::eigenvalues() const
1266  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1267       // check input       // check input
1268       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1269       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1270          throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");          throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1271       if(s[0] != s[1])       if(s[0] != s[1])
1272          throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");          throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1273       // create return       // create return
1274       DataArrayView::ShapeType ev_shape(1,s[0]);       DataArrayView::ShapeType ev_shape(1,s[0]);
# Line 1487  Data::eigenvalues() const Line 1281  Data::eigenvalues() const
1281  const boost::python::tuple  const boost::python::tuple
1282  Data::eigenvalues_and_eigenvectors(const double tol) const  Data::eigenvalues_and_eigenvectors(const double tol) const
1283  {  {
      #if defined DOPROF  
         profData->unary++;  
      #endif  
1284       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1285       if (getDataPointRank()!=2)       if (getDataPointRank()!=2)
1286          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.");
1287       if(s[0] != s[1])       if(s[0] != s[1])
1288          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1289       // create return       // create return
1290       DataArrayView::ShapeType ev_shape(1,s[0]);       DataArrayView::ShapeType ev_shape(1,s[0]);
# Line 1507  Data::eigenvalues_and_eigenvectors(const Line 1298  Data::eigenvalues_and_eigenvectors(const
1298  }  }
1299    
1300  const boost::python::tuple  const boost::python::tuple
1301  Data::mindp() const  Data::minGlobalDataPoint() const
1302  {  {
1303    // 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
1304    // 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
1305    // surrounding function    // surrounding function
1306    
   int SampleNo;  
1307    int DataPointNo;    int DataPointNo;
1308      int ProcNo;    int ProcNo;
1309    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1310    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1311  }  }
1312    
1313  void  void
1314  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1315                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1316  {  {
1317    int i,j;    int i,j;
1318    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1565  Data::calc_mindp(  int& ProcNo, Line 1354  Data::calc_mindp(  int& ProcNo,
1354      int lowProc = 0;      int lowProc = 0;
1355      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1356      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1357        
1358      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1359          next = globalMins[lowProc];          next = globalMins[lowProc];
1360          for( i=1; i<get_MPISize(); i++ )          for( i=1; i<get_MPISize(); i++ )
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1370  Data::calc_mindp(  int& ProcNo,
1370  #else  #else
1371      ProcNo = 0;      ProcNo = 0;
1372  #endif  #endif
1373    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1374  }  }
1375    
1376  void  void
# Line 1609  Data::operator+=(const Data& right) Line 1397  Data::operator+=(const Data& right)
1397    if (isProtected()) {    if (isProtected()) {
1398          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1399    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1400    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1401    return (*this);    return (*this);
1402  }  }
# Line 1619  Data::operator+=(const Data& right) Line 1404  Data::operator+=(const Data& right)
1404  Data&  Data&
1405  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1406  {  {
1407    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1408          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,plus<double>());
1409    }    return (*this);
1410  #if defined DOPROF  }
1411    profData->binary++;  Data&
1412  #endif  Data::operator=(const Data& other)
1413    binaryOp(right,plus<double>());  {
1414      copy(other);
1415    return (*this);    return (*this);
1416  }  }
1417    
# Line 1635  Data::operator-=(const Data& right) Line 1421  Data::operator-=(const Data& right)
1421    if (isProtected()) {    if (isProtected()) {
1422          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1423    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1424    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1425    return (*this);    return (*this);
1426  }  }
# Line 1645  Data::operator-=(const Data& right) Line 1428  Data::operator-=(const Data& right)
1428  Data&  Data&
1429  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1430  {  {
1431    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1432          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1433    return (*this);    return (*this);
1434  }  }
1435    
# Line 1661  Data::operator*=(const Data& right) Line 1439  Data::operator*=(const Data& right)
1439    if (isProtected()) {    if (isProtected()) {
1440          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1441    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1442    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1443    return (*this);    return (*this);
1444  }  }
# Line 1671  Data::operator*=(const Data& right) Line 1446  Data::operator*=(const Data& right)
1446  Data&  Data&
1447  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1448  {  {
1449    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1450          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1451    return (*this);    return (*this);
1452  }  }
1453    
# Line 1687  Data::operator/=(const Data& right) Line 1457  Data::operator/=(const Data& right)
1457    if (isProtected()) {    if (isProtected()) {
1458          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1459    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1460    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1461    return (*this);    return (*this);
1462  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1464  Data::operator/=(const Data& right)
1464  Data&  Data&
1465  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1466  {  {
1467    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1468          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1469    return (*this);    return (*this);
1470  }  }
1471    
1472  Data  Data
1473  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1474  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1475    Data left_d(left,*this);    Data left_d(left,*this);
1476    return left_d.powD(*this);    return left_d.powD(*this);
1477  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1479  Data::rpowO(const boost::python::object&
1479  Data  Data
1480  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1481  {  {
1482  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1483    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1484  }  }
1485    
1486  Data  Data
1487  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1488  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1489    Data result;    Data result;
1490    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1491    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1492         result.binaryOp(*this,escript::rpow);
1493      } else {
1494         result.copy(*this);
1495         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1496      }
1497    return result;    return result;
1498  }  }
1499    
# Line 1753  escript::operator+(const Data& left, con Line 1506  escript::operator+(const Data& left, con
1506    Data result;    Data result;
1507    //    //
1508    // perform a deep copy    // perform a deep copy
1509    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1510    result+=right;       result.copy(right);
1511         result+=left;
1512      } else {
1513         result.copy(left);
1514         result+=right;
1515      }
1516    return result;    return result;
1517  }  }
1518    
# Line 1766  escript::operator-(const Data& left, con Line 1524  escript::operator-(const Data& left, con
1524    Data result;    Data result;
1525    //    //
1526    // perform a deep copy    // perform a deep copy
1527    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1528    result-=right;       result=right.neg();
1529         result+=left;
1530      } else {
1531         result.copy(left);
1532         result-=right;
1533      }
1534    return result;    return result;
1535  }  }
1536    
# Line 1779  escript::operator*(const Data& left, con Line 1542  escript::operator*(const Data& left, con
1542    Data result;    Data result;
1543    //    //
1544    // perform a deep copy    // perform a deep copy
1545    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1546    result*=right;       result.copy(right);
1547         result*=left;
1548      } else {
1549         result.copy(left);
1550         result*=right;
1551      }
1552    return result;    return result;
1553  }  }
1554    
# Line 1792  escript::operator/(const Data& left, con Line 1560  escript::operator/(const Data& left, con
1560    Data result;    Data result;
1561    //    //
1562    // perform a deep copy    // perform a deep copy
1563    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1564    result/=right;       result=right.oneOver();
1565         result*=left;
1566      } else {
1567         result.copy(left);
1568         result/=right;
1569      }
1570    return result;    return result;
1571  }  }
1572    
# Line 1802  escript::operator/(const Data& left, con Line 1575  escript::operator/(const Data& left, con
1575  Data  Data
1576  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1577  {  {
1578    //    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;  
1579  }  }
1580    
1581  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1583  escript::operator+(const Data& left, con
1583  Data  Data
1584  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1585  {  {
1586    //    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;  
1587  }  }
1588    
1589  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1591  escript::operator-(const Data& left, con
1591  Data  Data
1592  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1593  {  {
1594    //    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;  
1595  }  }
1596    
1597  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1599  escript::operator*(const Data& left, con
1599  Data  Data
1600  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1601  {  {
1602    //    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;  
1603  }  }
1604    
1605  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1607  escript::operator/(const Data& left, con
1607  Data  Data
1608  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1609  {  {
1610    //    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;  
1611  }  }
1612    
1613  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1615  escript::operator+(const boost::python::
1615  Data  Data
1616  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1617  {  {
1618    //    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;  
1619  }  }
1620    
1621  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1623  escript::operator-(const boost::python::
1623  Data  Data
1624  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1625  {  {
1626    //    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;  
1627  }  }
1628    
1629  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1631  escript::operator*(const boost::python::
1631  Data  Data
1632  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1633  {  {
1634    //    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;  
1635  }  }
1636    
1637  //  //
# Line 1961  escript::operator/(const boost::python:: Line 1682  escript::operator/(const boost::python::
1682  /* TODO */  /* TODO */
1683  /* global reduction */  /* global reduction */
1684  Data  Data
1685  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1686  {  {
1687    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1688    
# Line 1979  Data::getItem(const boost::python::objec Line 1700  Data::getItem(const boost::python::objec
1700  Data  Data
1701  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1702  {  {
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
1703    return Data(*this,region);    return Data(*this,region);
1704  }  }
1705    
# Line 1995  Data::setItemO(const boost::python::obje Line 1713  Data::setItemO(const boost::python::obje
1713    setItemD(key,tempData);    setItemD(key,tempData);
1714  }  }
1715    
 /* TODO */  
 /* global reduction */  
1716  void  void
1717  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1718                 const Data& value)                 const Data& value)
# Line 2014  Data::setItemD(const boost::python::obje Line 1730  Data::setItemD(const boost::python::obje
1730    }    }
1731  }  }
1732    
 /* TODO */  
 /* global reduction */  
1733  void  void
1734  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1735                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
# Line 2023  Data::setSlice(const Data& value, Line 1737  Data::setSlice(const Data& value,
1737    if (isProtected()) {    if (isProtected()) {
1738          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1739    }    }
 #if defined DOPROF  
   profData->slicing++;  
 #endif  
1740    Data tempValue(value);    Data tempValue(value);
1741    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1742    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 2060  Data::typeMatchRight(const Data& right) Line 1771  Data::typeMatchRight(const Data& right)
1771    }    }
1772  }  }
1773    
1774  /* TODO */  void
1775  /* global reduction */  Data::setTaggedValueByName(std::string name,
1776                               const boost::python::object& value)
1777    {
1778         if (getFunctionSpace().getDomain().isValidTagName(name)) {
1779            int tagKey=getFunctionSpace().getDomain().getTag(name);
1780            setTaggedValue(tagKey,value);
1781         }
1782    }
1783  void  void
1784  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1785                       const boost::python::object& value)                       const boost::python::object& value)
# Line 2086  Data::setTaggedValue(int tagKey, Line 1804  Data::setTaggedValue(int tagKey,
1804    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1805  }  }
1806    
 /* TODO */  
 /* global reduction */  
1807  void  void
1808  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1809                              const DataArrayView& value)                              const DataArrayView& value)
# Line 2102  Data::setTaggedValueFromCPP(int tagKey, Line 1818  Data::setTaggedValueFromCPP(int tagKey,
1818    if (!isTagged()) {    if (!isTagged()) {
1819      throw DataException("Error - DataTagged conversion failed!!");      throw DataException("Error - DataTagged conversion failed!!");
1820    }    }
1821                                                                                                                  
1822    //    //
1823    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1824    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1825  }  }
1826    
 /* TODO */  
 /* global reduction */  
1827  int  int
1828  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
1829  {  {
1830    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
1831  }  }
1832    
 /* TODO */  
 /* global reduction */  
 void  
 Data::setRefValue(int ref,  
                   const boost::python::numeric::array& value)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
   //  
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 /* TODO */  
 /* global reduction */  
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
   
   //  
   // Load DataArray with values from data-points specified by ref  
   m_data->getRefValue(ref,valueDataArray);  
   
   //  
   // Load values from valueDataArray into return numarray  
   
   // extract the shape of the numarray  
   int rank = value.getrank();  
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
   
   if (rank==0) {  
       boost::python::numeric::array temp_numArray(valueView());  
       value = temp_numArray;  
   }  
   if (rank==1) {  
     for (int i=0; i < shape[0]; i++) {  
       value[i] = valueView(i);  
     }  
   }  
   if (rank==2) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         value[i][j] = valueView(i,j);  
       }  
     }  
   }  
   if (rank==3) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           value[i][j][k] = valueView(i,j,k);  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           for (int l=0; l < shape[3]; l++) {  
             value[i][j][k][l] = valueView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
1833  void  void
1834  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
1835  {  {
# Line 2241  Data::archiveData(const std::string file Line 1871  Data::archiveData(const std::string file
1871    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
1872    vector<int> referenceNumbers(noSamples);    vector<int> referenceNumbers(noSamples);
1873    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1874      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1875    }    }
1876    vector<int> tagNumbers(noSamples);    vector<int> tagNumbers(noSamples);
1877    if (isTagged()) {    if (isTagged()) {
# Line 2450  Data::extractData(const std::string file Line 2080  Data::extractData(const std::string file
2080      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2081    }    }
2082    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2083      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2084        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2085      }      }
2086    }    }
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2171  ostream& escript::operator<<(ostream& o,
2171    return o;    return o;
2172  }  }
2173    
2174    Data
2175    escript::C_GeneralTensorProduct(Data& arg_0,
2176                         Data& arg_1,
2177                         int axis_offset,
2178                         int transpose)
2179    {
2180      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2181      // SM is the product of the last axis_offset entries in arg_0.getShape().
2182    
2183      // Interpolate if necessary and find an appropriate function space
2184      Data arg_0_Z, arg_1_Z;
2185      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2186        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2187          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2188          arg_1_Z = Data(arg_1);
2189        }
2190        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2191          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2192          arg_0_Z =Data(arg_0);
2193        }
2194        else {
2195          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2196        }
2197      } else {
2198          arg_0_Z = Data(arg_0);
2199          arg_1_Z = Data(arg_1);
2200      }
2201      // Get rank and shape of inputs
2202      int rank0 = arg_0_Z.getDataPointRank();
2203      int rank1 = arg_1_Z.getDataPointRank();
2204      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2205      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2206    
2207      // Prepare for the loops of the product and verify compatibility of shapes
2208      int start0=0, start1=0;
2209      if (transpose == 0)       {}
2210      else if (transpose == 1)  { start0 = axis_offset; }
2211      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2212      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2213    
2214      // Adjust the shapes for transpose
2215      DataArrayView::ShapeType tmpShape0;
2216      DataArrayView::ShapeType tmpShape1;
2217      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2218      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2219    
2220    #if 0
2221      // For debugging: show shape after transpose
2222      char tmp[100];
2223      std::string shapeStr;
2224      shapeStr = "(";
2225      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2226      shapeStr += ")";
2227      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2228      shapeStr = "(";
2229      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2230      shapeStr += ")";
2231      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2232    #endif
2233    
2234      // Prepare for the loops of the product
2235      int SL=1, SM=1, SR=1;
2236      for (int i=0; i<rank0-axis_offset; i++)   {
2237        SL *= tmpShape0[i];
2238      }
2239      for (int i=rank0-axis_offset; i<rank0; i++)   {
2240        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2241          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2242        }
2243        SM *= tmpShape0[i];
2244      }
2245      for (int i=axis_offset; i<rank1; i++)     {
2246        SR *= tmpShape1[i];
2247      }
2248    
2249      // Define the shape of the output
2250      DataArrayView::ShapeType shape2;
2251      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2252      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2253    
2254      // Declare output Data object
2255      Data res;
2256    
2257      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2258        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2259        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2260        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2261        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2262        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2263      }
2264      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2265    
2266        // Prepare the DataConstant input
2267        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2268        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2269    
2270        // Borrow DataTagged input from Data object
2271        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2272        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2273    
2274        // Prepare a DataTagged output 2
2275        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2276        res.tag();
2277        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2278        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2279    
2280        // Prepare offset into DataConstant
2281        int offset_0 = tmp_0->getPointOffset(0,0);
2282        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2283        // Get the views
2284        DataArrayView view_1 = tmp_1->getDefaultValue();
2285        DataArrayView view_2 = tmp_2->getDefaultValue();
2286        // Get the pointers to the actual data
2287        double *ptr_1 = &((view_1.getData())[0]);
2288        double *ptr_2 = &((view_2.getData())[0]);
2289        // Compute an MVP for the default
2290        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2291        // Compute an MVP for each tag
2292        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2293        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2294        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2295          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2296          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2297          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2298          double *ptr_1 = &view_1.getData(0);
2299          double *ptr_2 = &view_2.getData(0);
2300          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2301        }
2302    
2303      }
2304      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2305    
2306        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2307        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2308        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2309        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2310        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2311        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2312        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2313        int sampleNo_1,dataPointNo_1;
2314        int numSamples_1 = arg_1_Z.getNumSamples();
2315        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2316        int offset_0 = tmp_0->getPointOffset(0,0);
2317        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2318        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2319          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2320            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2321            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2322            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2323            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2324            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2325            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2326          }
2327        }
2328    
2329      }
2330      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2331    
2332        // Borrow DataTagged input from Data object
2333        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2334        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2335    
2336        // Prepare the DataConstant input
2337        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2338        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2339    
2340        // Prepare a DataTagged output 2
2341        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2342        res.tag();
2343        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2344        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2345    
2346        // Prepare offset into DataConstant
2347        int offset_1 = tmp_1->getPointOffset(0,0);
2348        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2349        // Get the views
2350        DataArrayView view_0 = tmp_0->getDefaultValue();
2351        DataArrayView view_2 = tmp_2->getDefaultValue();
2352        // Get the pointers to the actual data
2353        double *ptr_0 = &((view_0.getData())[0]);
2354        double *ptr_2 = &((view_2.getData())[0]);
2355        // Compute an MVP for the default
2356        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2357        // Compute an MVP for each tag
2358        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2359        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2360        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2361          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2362          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2363          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2364          double *ptr_0 = &view_0.getData(0);
2365          double *ptr_2 = &view_2.getData(0);
2366          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2367        }
2368    
2369      }
2370      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2371    
2372        // Borrow DataTagged input from Data object
2373        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2374        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2375    
2376        // Borrow DataTagged input from Data object
2377        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2378        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2379    
2380        // Prepare a DataTagged output 2
2381        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2382        res.tag();  // DataTagged output
2383        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2384        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2385    
2386        // Get the views
2387        DataArrayView view_0 = tmp_0->getDefaultValue();
2388        DataArrayView view_1 = tmp_1->getDefaultValue();
2389        DataArrayView view_2 = tmp_2->getDefaultValue();
2390        // Get the pointers to the actual data
2391        double *ptr_0 = &((view_0.getData())[0]);
2392        double *ptr_1 = &((view_1.getData())[0]);
2393        double *ptr_2 = &((view_2.getData())[0]);
2394        // Compute an MVP for the default
2395        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2396        // Merge the tags
2397        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2398        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2399        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2400        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2401          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2402        }
2403        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2404          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2405        }
2406        // Compute an MVP for each tag
2407        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2408        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2409          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2410          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2411          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2412          double *ptr_0 = &view_0.getData(0);
2413          double *ptr_1 = &view_1.getData(0);
2414          double *ptr_2 = &view_2.getData(0);
2415          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2416        }
2417    
2418      }
2419      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2420    
2421        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2422        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2423        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2424        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2425        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2426        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2427        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2428        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2429        int sampleNo_0,dataPointNo_0;
2430        int numSamples_0 = arg_0_Z.getNumSamples();
2431        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2432        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2433        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2434          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2435          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2436          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2437            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2438            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2439            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2440            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2441            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2442          }
2443        }
2444    
2445      }
2446      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2447    
2448        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2449        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2450        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2451        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2452        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2453        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2454        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2455        int sampleNo_0,dataPointNo_0;
2456        int numSamples_0 = arg_0_Z.getNumSamples();
2457        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2458        int offset_1 = tmp_1->getPointOffset(0,0);
2459        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2460        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
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_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2466            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2467            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2468          }
2469        }
2470    
2471    
2472      }
2473      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2474    
2475        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2476        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2477        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2478        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2479        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2480        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2481        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2482        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2483        int sampleNo_0,dataPointNo_0;
2484        int numSamples_0 = arg_0_Z.getNumSamples();
2485        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2486        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2487        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2488          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2489          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2490          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2491            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2492            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2493            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2494            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2495            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2496          }
2497        }
2498    
2499      }
2500      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2501    
2502        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2503        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2504        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2505        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2506        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2507        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2508        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2509        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2510        int sampleNo_0,dataPointNo_0;
2511        int numSamples_0 = arg_0_Z.getNumSamples();
2512        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2513        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2514        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2515          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2516            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2517            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2518            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2519            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2520            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2521            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2522            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2523          }
2524        }
2525    
2526      }
2527      else {
2528        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2529      }
2530    
2531      return res;
2532    }
2533    
2534    DataAbstract*
2535    Data::borrowData() const
2536    {
2537      return m_data.get();
2538    }
2539    
2540  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2541    
2542  void  void
2543  Data::print()  Data::print()
2544  {  {
2545    int i,j;    int i,j;
2546      
2547    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );    printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2548    for( i=0; i<getNumSamples(); i++ )    for( i=0; i<getNumSamples(); i++ )
2549    {    {
# Line 2557  Data::print() Line 2553  Data::print()
2553      printf( "\n" );      printf( "\n" );
2554    }    }
2555  }  }
2556    void
2557    Data::dump(const std::string fileName) const
2558    {
2559      try
2560         {
2561            return m_data->dump(fileName);
2562         }
2563         catch (exception& e)
2564         {
2565            cout << e.what() << endl;
2566         }
2567    }
2568    
2569  int  int
2570  Data::get_MPISize() const  Data::get_MPISize() const
# Line 2584  Data::get_MPIRank() const Line 2592  Data::get_MPIRank() const
2592    
2593  MPI_Comm  MPI_Comm
2594  Data::get_MPIComm() const  Data::get_MPIComm() const
2595  {  {
2596  #ifdef PASO_MPI  #ifdef PASO_MPI
2597      return MPI_COMM_WORLD;      return MPI_COMM_WORLD;
2598  #else  #else

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

  ViewVC Help
Powered by ViewVC 1.1.26