/[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 474 by jgs, Mon Jan 30 04:23:44 2006 UTC revision 1031 by phornby, Wed Mar 14 06:03:21 2007 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
   
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  
  *                                                                            *  
  * This software is the property of ACcESS.  No part of this code             *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that                          *  
  * person has a software license agreement with ACcESS.                       *  
  *                                                                            *  
  ******************************************************************************  
   
 ******************************************************************************/  
2    
3    /*
4     ************************************************************
5     *          Copyright 2006 by ACcESS MNRF                   *
6     *                                                          *
7     *              http://www.access.edu.au                    *
8     *       Primary Business: Queensland, Australia            *
9     *  Licensed under the Open Software License version 3.0    *
10     *     http://www.opensource.org/licenses/osl-3.0.php       *
11     *                                                          *
12     ************************************************************
13    */
14  #include "Data.h"  #include "Data.h"
15    
 #include <iostream>  
 #include <fstream>  
 #include <algorithm>  
 #include <vector>  
 #include <exception>  
 #include <functional>  
 #include <math.h>  
   
 #include <boost/python/dict.hpp>  
 #include <boost/python/str.hpp>  
 #include <boost/python/extract.hpp>  
 #include <boost/python/long.hpp>  
 #include <boost/python/tuple.hpp>  
   
 #include "DataException.h"  
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataConstant.h"  #include "DataConstant.h"
18  #include "DataTagged.h"  #include "DataTagged.h"
19  #include "DataEmpty.h"  #include "DataEmpty.h"
20  #include "DataArray.h"  #include "DataArray.h"
21    #include "DataArrayView.h"
22  #include "DataProf.h"  #include "DataProf.h"
23  #include "FunctionSpaceFactory.h"  #include "FunctionSpaceFactory.h"
24  #include "AbstractContinuousDomain.h"  #include "AbstractContinuousDomain.h"
25  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
26    
27    #include <fstream>
28    #include <algorithm>
29    #include <vector>
30    #include <functional>
31    
32    #include <boost/python/dict.hpp>
33    #include <boost/python/extract.hpp>
34    #include <boost/python/long.hpp>
35    
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
# Line 60  Data::Data() Line 51  Data::Data()
51    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
52    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
53    m_data=temp_data;    m_data=temp_data;
54      m_protected=false;
55  #if defined DOPROF  #if defined DOPROF
56    // create entry in global profiling table for this object    // create entry in global profiling table for this object
57    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 77  Data::Data(double value, Line 69  Data::Data(double value,
69    }    }
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
72      m_protected=false;
73  #if defined DOPROF  #if defined DOPROF
74    // create entry in global profiling table for this object    // create entry in global profiling table for this object
75    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 91  Data::Data(double value, Line 84  Data::Data(double value,
84    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
85    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
86    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
87      m_protected=false;
88  #if defined DOPROF  #if defined DOPROF
89    // create entry in global profiling table for this object    // create entry in global profiling table for this object
90    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 100  Data::Data(double value, Line 94  Data::Data(double value,
94  Data::Data(const Data& inData)  Data::Data(const Data& inData)
95  {  {
96    m_data=inData.m_data;    m_data=inData.m_data;
97      m_protected=inData.isProtected();
98  #if defined DOPROF  #if defined DOPROF
99    // create entry in global profiling table for this object    // create entry in global profiling table for this object
100    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 114  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
110    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
111    m_data=temp_data;    m_data=temp_data;
112      m_protected=false;
113  #if defined DOPROF  #if defined DOPROF
114    // create entry in global profiling table for this object    // create entry in global profiling table for this object
115    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 123  Data::Data(const Data& inData, Line 119  Data::Data(const Data& inData,
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122    #if defined DOPROF
123      // create entry in global profiling table for this object
124      profData = dataProfTable.newData();
125    #endif
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    } else {
129        #if defined DOPROF
130        profData->interpolate++;
131        #endif
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
133      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
134      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 139  Data::Data(const Data& inData, Line 142  Data::Data(const Data& inData,
142      }      }
143      m_data=tmp.m_data;      m_data=tmp.m_data;
144    }    }
145  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 154  Data::Data(const DataTagged::TagListType Line 154  Data::Data(const DataTagged::TagListType
154    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
155    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
156    m_data=temp_data;    m_data=temp_data;
157      m_protected=false;
158    if (expanded) {    if (expanded) {
159      expand();      expand();
160    }    }
# Line 168  Data::Data(const numeric::array& value, Line 169  Data::Data(const numeric::array& value,
169             bool expanded)             bool expanded)
170  {  {
171    initialise(value,what,expanded);    initialise(value,what,expanded);
172      m_protected=false;
173  #if defined DOPROF  #if defined DOPROF
174    // create entry in global profiling table for this object    // create entry in global profiling table for this object
175    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 179  Data::Data(const DataArrayView& value, Line 181  Data::Data(const DataArrayView& value,
181             bool expanded)             bool expanded)
182  {  {
183    initialise(value,what,expanded);    initialise(value,what,expanded);
184      m_protected=false;
185  #if defined DOPROF  #if defined DOPROF
186    // create entry in global profiling table for this object    // create entry in global profiling table for this object
187    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 191  Data::Data(const object& value, Line 194  Data::Data(const object& value,
194  {  {
195    numeric::array asNumArray(value);    numeric::array asNumArray(value);
196    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
197      m_protected=false;
198  #if defined DOPROF  #if defined DOPROF
199    // create entry in global profiling table for this object    // create entry in global profiling table for this object
200    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 215  Data::Data(const object& value, Line 219  Data::Data(const object& value,
219      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
220      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
221    }    }
222      m_protected=false;
223  #if defined DOPROF  #if defined DOPROF
224    // create entry in global profiling table for this object    // create entry in global profiling table for this object
225    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 360  Data::isConstant() const Line 365  Data::isConstant() const
365  }  }
366    
367  void  void
368    Data::setProtection()
369    {
370       m_protected=true;
371    }
372    
373    bool
374    Data::isProtected() const
375    {
376       return m_protected;
377    }
378    
379    
380    
381    void
382  Data::expand()  Data::expand()
383  {  {
384    if (isConstant()) {    if (isConstant()) {
# Line 401  Data::tag() Line 420  Data::tag()
420    }    }
421  }  }
422    
423  void  Data
424  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
425  {  {
426    m_data->reshapeDataPoint(shape);  #if defined DOPROF
427      profData->where++;
428    #endif
429      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
430  }  }
431    
432  Data  Data
# Line 444  Data::whereNonPositive() const Line 466  Data::whereNonPositive() const
466  }  }
467    
468  Data  Data
469  Data::whereZero() const  Data::whereZero(double tol) const
470  {  {
471  #if defined DOPROF  #if defined DOPROF
472    profData->where++;    profData->where++;
473  #endif  #endif
474    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
475      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
476  }  }
477    
478  Data  Data
479  Data::whereNonZero() const  Data::whereNonZero(double tol) const
480  {  {
481  #if defined DOPROF  #if defined DOPROF
482    profData->where++;    profData->where++;
483  #endif  #endif
484    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
485      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
486  }  }
487    
488  Data  Data
# Line 524  Data::getDataPointShape() const Line 548  Data::getDataPointShape() const
548    return getPointDataView().getShape();    return getPointDataView().getShape();
549  }  }
550    
551    
552    
553  void  void
554  Data::fillFromNumArray(const boost::python::numeric::array num_array)  Data::fillFromNumArray(const boost::python::numeric::array num_array)
555  {  {
556      if (isProtected()) {
557            throw DataException("Error - attempt to update protected Data object.");
558      }
559    //    //
560    // check rank    // check rank
561    if (num_array.getrank()<getDataPointRank())    if (num_array.getrank()<getDataPointRank())
# Line 613  Data::convertToNumArray() Line 642  Data::convertToNumArray()
642        }        }
643        if (dataPointRank==1) {        if (dataPointRank==1) {
644          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
645            numArray[dataPoint][i]=dataPointView(i);            numArray[make_tuple(dataPoint,i)]=dataPointView(i);
646          }          }
647        }        }
648        if (dataPointRank==2) {        if (dataPointRank==2) {
649          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
650            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
651              numArray[dataPoint][i][j] = dataPointView(i,j);              numArray[make_tuple(dataPoint,i,j)] = dataPointView(i,j);
652            }            }
653          }          }
654        }        }
# Line 627  Data::convertToNumArray() Line 656  Data::convertToNumArray()
656          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
657            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
658              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
659                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                numArray[make_tuple(dataPoint,i,j,k)]=dataPointView(i,j,k);
660              }              }
661            }            }
662          }          }
# Line 637  Data::convertToNumArray() Line 666  Data::convertToNumArray()
666            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
667              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
668                for (int l=0; l<dataPointShape[3]; l++) {                for (int l=0; l<dataPointShape[3]; l++) {
669                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  numArray[make_tuple(dataPoint,i,j,k,l)]=dataPointView(i,j,k,l);
670                }                }
671              }              }
672            }            }
# Line 652  Data::convertToNumArray() Line 681  Data::convertToNumArray()
681    return numArray;    return numArray;
682  }  }
683    
684  const  
685    const
686  boost::python::numeric::array  boost::python::numeric::array
687  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data:: getValueOfDataPoint(int dataPointNo)
688  {  {
689    //    size_t length=0;
690    // Check a valid sample number has been supplied    int i, j, k, l;
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   
691    //    //
692    // determine the rank and shape of each data point    // determine the rank and shape of each data point
693    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 676  Data::convertToNumArrayFromSampleNo(int Line 698  Data::convertToNumArrayFromSampleNo(int
698    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
699    
700    //    //
701    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
702    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
703    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
704    // 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(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
705    
706    //    //
707    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
708      if (arrayRank==0) {
709        numArray.resize(1);
710      }
711    if (arrayRank==1) {    if (arrayRank==1) {
712      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
713    }    }
# Line 703  Data::convertToNumArrayFromSampleNo(int Line 720  Data::convertToNumArrayFromSampleNo(int
720    if (arrayRank==4) {    if (arrayRank==4) {
721      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
722    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
723    
724    //    if (getNumDataPointsPerSample()>0) {
725    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
726    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
727    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {         //
728      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         // Check a valid sample number has been supplied
729      if (dataPointRank==0) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
730        numArray[dataPoint]=dataPointView();             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
731      }         }
732      if (dataPointRank==1) {                
733        for (int i=0; i<dataPointShape[0]; i++) {         //
734          numArray[dataPoint][i]=dataPointView(i);         // Check a valid data point number has been supplied
735        }         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
736      }             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
737      if (dataPointRank==2) {         }
738        for (int i=0; i<dataPointShape[0]; i++) {         // TODO: global error handling
739          for (int j=0; j<dataPointShape[1]; j++) {         // create a view of the data if it is stored locally
740            numArray[dataPoint][i][j] = dataPointView(i,j);         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
741          }          
742        }         switch( dataPointRank ){
743      }              case 0 :
744      if (dataPointRank==3) {                  numArray[0] = dataPointView();
745        for (int i=0; i<dataPointShape[0]; i++) {                  break;
746          for (int j=0; j<dataPointShape[1]; j++) {              case 1 :        
747            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
748              numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                      numArray[i]=dataPointView(i);
749            }                  break;
750          }              case 2 :        
751        }                  for( i=0; i<dataPointShape[0]; i++ )
752      }                      for( j=0; j<dataPointShape[1]; j++)
753      if (dataPointRank==4) {                          numArray[make_tuple(i,j)]=dataPointView(i,j);
754        for (int i=0; i<dataPointShape[0]; i++) {                  break;
755          for (int j=0; j<dataPointShape[1]; j++) {              case 3 :        
756            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
757              for (int l=0; l<dataPointShape[3]; l++) {                      for( j=0; j<dataPointShape[1]; j++ )
758                numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                          for( k=0; k<dataPointShape[2]; k++)
759              }                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
760            }                  break;
761          }              case 4 :
762        }                  for( i=0; i<dataPointShape[0]; i++ )
763      }                      for( j=0; j<dataPointShape[1]; j++ )
764                            for( k=0; k<dataPointShape[2]; k++ )
765                                for( l=0; l<dataPointShape[3]; l++)
766                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
767                    break;
768        }
769    }    }
   
770    //    //
771    // return the loaded array    // return the array
772    return numArray;    return numArray;
 }  
773    
774  const  }
775  boost::python::numeric::array  void
776  Data::convertToNumArrayFromDPNo(int sampleNo,  Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array num_array)
                                 int dataPointNo)  
777  {  {
778      if (isProtected()) {
779            throw DataException("Error - attempt to update protected Data object.");
780      }
781      //
782      // check rank
783      if (num_array.getrank()<getDataPointRank())
784          throw DataException("Rank of numarray does not match Data object rank");
785    
786    //    //
787    // Check a valid sample number has been supplied    // check shape of num_array
788    if (sampleNo >= getNumSamples()) {    for (int i=0; i<getDataPointRank(); i++) {
789      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
790           throw DataException("Shape of numarray does not match Data object rank");
791    }    }
792      //
793      // make sure data is expanded:
794      if (!isExpanded()) {
795        expand();
796      }
797      if (getNumDataPointsPerSample()>0) {
798           int sampleNo = dataPointNo/getNumDataPointsPerSample();
799           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
800           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
801      } else {
802           m_data->copyToDataPoint(-1, 0,num_array);
803      }
804    }
805    
806    void
807    Data::setValueOfDataPoint(int dataPointNo, const double value)
808    {
809      if (isProtected()) {
810            throw DataException("Error - attempt to update protected Data object.");
811      }
812    //    //
813    // Check a valid data point number has been supplied    // make sure data is expanded:
814    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
815      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
816      }
817      if (getNumDataPointsPerSample()>0) {
818           int sampleNo = dataPointNo/getNumDataPointsPerSample();
819           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
820           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
821      } else {
822           m_data->copyToDataPoint(-1, 0,value);
823    }    }
824    }
825    
826    const
827    boost::python::numeric::array
828    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
829    {
830      size_t length=0;
831      int i, j, k, l, pos;
832    //    //
833    // determine the rank and shape of each data point    // determine the rank and shape of each data point
834    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 804  Data::convertToNumArrayFromDPNo(int samp Line 862  Data::convertToNumArrayFromDPNo(int samp
862      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
863    }    }
864    
865      // added for the MPI communication
866      length=1;
867      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
868      double *tmpData = new double[length];
869    
870    //    //
871    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
   DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
872    
873        // updated for the MPI case
874        if( get_MPIRank()==procNo ){
875                 if (getNumDataPointsPerSample()>0) {
876                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
877                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
878                    //
879                    // Check a valid sample number has been supplied
880                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
881                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
882                    }
883                  
884                    //
885                    // Check a valid data point number has been supplied
886                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
887                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
888                    }
889                    // TODO: global error handling
890            // create a view of the data if it is stored locally
891            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
892            
893            // pack the data from the view into tmpData for MPI communication
894            pos=0;
895            switch( dataPointRank ){
896                case 0 :
897                    tmpData[0] = dataPointView();
898                    break;
899                case 1 :        
900                    for( i=0; i<dataPointShape[0]; i++ )
901                        tmpData[i]=dataPointView(i);
902                    break;
903                case 2 :        
904                    for( i=0; i<dataPointShape[0]; i++ )
905                        for( j=0; j<dataPointShape[1]; j++, pos++ )
906                            tmpData[pos]=dataPointView(i,j);
907                    break;
908                case 3 :        
909                    for( i=0; i<dataPointShape[0]; i++ )
910                        for( j=0; j<dataPointShape[1]; j++ )
911                            for( k=0; k<dataPointShape[2]; k++, pos++ )
912                                tmpData[pos]=dataPointView(i,j,k);
913                    break;
914                case 4 :
915                    for( i=0; i<dataPointShape[0]; i++ )
916                        for( j=0; j<dataPointShape[1]; j++ )
917                            for( k=0; k<dataPointShape[2]; k++ )
918                                for( l=0; l<dataPointShape[3]; l++, pos++ )
919                                    tmpData[pos]=dataPointView(i,j,k,l);
920                    break;
921            }
922                }
923        }
924            #ifdef PASO_MPI
925            // broadcast the data to all other processes
926        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
927            #endif
928    
929        // unpack the data
930        switch( dataPointRank ){
931            case 0 :
932                numArray[0]=tmpData[0];
933                break;
934            case 1 :        
935                for( i=0; i<dataPointShape[0]; i++ )
936                    numArray[i]=tmpData[i];
937                break;
938            case 2 :        
939                for( i=0; i<dataPointShape[0]; i++ )
940                    for( j=0; j<dataPointShape[1]; j++ )
941                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
942                break;
943            case 3 :        
944                for( i=0; i<dataPointShape[0]; i++ )
945                    for( j=0; j<dataPointShape[1]; j++ )
946                        for( k=0; k<dataPointShape[2]; k++ )
947                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
948                break;
949            case 4 :
950                for( i=0; i<dataPointShape[0]; i++ )
951                    for( j=0; j<dataPointShape[1]; j++ )
952                        for( k=0; k<dataPointShape[2]; k++ )
953                            for( l=0; l<dataPointShape[3]; l++ )
954                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
955                break;
956        }
957    
958        delete [] tmpData;  
959    //    //
960    // return the loaded array    // return the loaded array
961    return numArray;    return numArray;
962  }  }
963    
964    
965    
966  boost::python::numeric::array  boost::python::numeric::array
967  Data::integrate() const  Data::integrate() const
968  {  {
# Line 964  Data::acos() const Line 1079  Data::acos() const
1079    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1080  }  }
1081    
1082    
1083  Data  Data
1084  Data::atan() const  Data::atan() const
1085  {  {
# Line 1000  Data::tanh() const Line 1116  Data::tanh() const
1116    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1117  }  }
1118    
1119    
1120    Data
1121    Data::erf() const
1122    {
1123    #if defined DOPROF
1124      profData->unary++;
1125    #endif
1126    #ifdef _WIN32
1127      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1128    #else
1129      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1130    #endif
1131    }
1132    
1133  Data  Data
1134  Data::asinh() const  Data::asinh() const
1135  {  {
# Line 1105  Data::sqrt() const Line 1235  Data::sqrt() const
1235  double  double
1236  Data::Lsup() const  Data::Lsup() const
1237  {  {
1238      double localValue, globalValue;
1239  #if defined DOPROF  #if defined DOPROF
1240    profData->reduction1++;    profData->reduction1++;
1241  #endif  #endif
1242    //    //
1243    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1244    
1245    AbsMax abs_max_func;    AbsMax abs_max_func;
1246    return algorithm(abs_max_func,0);    localValue = algorithm(abs_max_func,0);
1247    #ifdef PASO_MPI
1248      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1249      return globalValue;
1250    #else
1251      return localValue;
1252    #endif
1253  }  }
1254    
1255  double  double
1256  Data::Linf() const  Data::Linf() const
1257  {  {
1258      double localValue, globalValue;
1259  #if defined DOPROF  #if defined DOPROF
1260    profData->reduction1++;    profData->reduction1++;
1261  #endif  #endif
1262    //    //
1263    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1264    AbsMin abs_min_func;    AbsMin abs_min_func;
1265    return algorithm(abs_min_func,numeric_limits<double>::max());    localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1266    
1267    #ifdef PASO_MPI
1268      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1269      return globalValue;
1270    #else
1271      return localValue;
1272    #endif
1273  }  }
1274    
1275  double  double
1276  Data::sup() const  Data::sup() const
1277  {  {
1278      double localValue, globalValue;
1279  #if defined DOPROF  #if defined DOPROF
1280    profData->reduction1++;    profData->reduction1++;
1281  #endif  #endif
1282    //    //
1283    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1284    FMax fmax_func;    FMax fmax_func;
1285    return algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1286    #ifdef PASO_MPI
1287      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1288      return globalValue;
1289    #else
1290      return localValue;
1291    #endif
1292  }  }
1293    
1294  double  double
1295  Data::inf() const  Data::inf() const
1296  {  {
1297      double localValue, globalValue;
1298  #if defined DOPROF  #if defined DOPROF
1299    profData->reduction1++;    profData->reduction1++;
1300  #endif  #endif
1301    //    //
1302    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1303    FMin fmin_func;    FMin fmin_func;
1304    return algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1305    #ifdef PASO_MPI
1306      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1307      return globalValue;
1308    #else
1309      return localValue;
1310    #endif
1311  }  }
1312    
1313    /* TODO */
1314    /* global reduction */
1315  Data  Data
1316  Data::maxval() const  Data::maxval() const
1317  {  {
# Line 1175  Data::minval() const Line 1337  Data::minval() const
1337  }  }
1338    
1339  Data  Data
1340  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1341  {  {
1342  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1343    profData->reduction2++;       #if defined DOPROF
1344  #endif       profData->unary++;
1345    Trace trace_func;       #endif
1346    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1347         DataArrayView::ShapeType ev_shape;
1348         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1349         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1350         int rank=getDataPointRank();
1351         if (rank<2) {
1352            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1353         }
1354         if (axis0<0 || axis0>rank-1) {
1355            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1356         }
1357         if (axis1<0 || axis1>rank-1) {
1358             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1359         }
1360         if (axis0 == axis1) {
1361             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1362         }
1363         if (axis0 > axis1) {
1364             axis0_tmp=axis1;
1365             axis1_tmp=axis0;
1366         } else {
1367             axis0_tmp=axis0;
1368             axis1_tmp=axis1;
1369         }
1370         for (int i=0; i<rank; i++) {
1371           if (i == axis0_tmp) {
1372              ev_shape.push_back(s[axis1_tmp]);
1373           } else if (i == axis1_tmp) {
1374              ev_shape.push_back(s[axis0_tmp]);
1375           } else {
1376              ev_shape.push_back(s[i]);
1377           }
1378         }
1379         Data ev(0.,ev_shape,getFunctionSpace());
1380         ev.typeMatchRight(*this);
1381         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1382         return ev;
1383    
1384    }
1385    
1386    Data
1387    Data::symmetric() const
1388    {
1389         #if defined DOPROF
1390            profData->unary++;
1391         #endif
1392         // check input
1393         DataArrayView::ShapeType s=getDataPointShape();
1394         if (getDataPointRank()==2) {
1395            if(s[0] != s[1])
1396               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1397         }
1398         else if (getDataPointRank()==4) {
1399            if(!(s[0] == s[2] && s[1] == s[3]))
1400               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1401         }
1402         else {
1403            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1404         }
1405         Data ev(0.,getDataPointShape(),getFunctionSpace());
1406         ev.typeMatchRight(*this);
1407         m_data->symmetric(ev.m_data.get());
1408         return ev;
1409    }
1410    
1411    Data
1412    Data::nonsymmetric() const
1413    {
1414         #if defined DOPROF
1415            profData->unary++;
1416         #endif
1417         // check input
1418         DataArrayView::ShapeType s=getDataPointShape();
1419         if (getDataPointRank()==2) {
1420            if(s[0] != s[1])
1421               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1422            DataArrayView::ShapeType ev_shape;
1423            ev_shape.push_back(s[0]);
1424            ev_shape.push_back(s[1]);
1425            Data ev(0.,ev_shape,getFunctionSpace());
1426            ev.typeMatchRight(*this);
1427            m_data->nonsymmetric(ev.m_data.get());
1428            return ev;
1429         }
1430         else if (getDataPointRank()==4) {
1431            if(!(s[0] == s[2] && s[1] == s[3]))
1432               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1433            DataArrayView::ShapeType ev_shape;
1434            ev_shape.push_back(s[0]);
1435            ev_shape.push_back(s[1]);
1436            ev_shape.push_back(s[2]);
1437            ev_shape.push_back(s[3]);
1438            Data ev(0.,ev_shape,getFunctionSpace());
1439            ev.typeMatchRight(*this);
1440            m_data->nonsymmetric(ev.m_data.get());
1441            return ev;
1442         }
1443         else {
1444            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1445         }
1446    }
1447    
1448    Data
1449    Data::trace(int axis_offset) const
1450    {
1451         #if defined DOPROF
1452            profData->unary++;
1453         #endif
1454         DataArrayView::ShapeType s=getDataPointShape();
1455         if (getDataPointRank()==2) {
1456            DataArrayView::ShapeType ev_shape;
1457            Data ev(0.,ev_shape,getFunctionSpace());
1458            ev.typeMatchRight(*this);
1459            m_data->trace(ev.m_data.get(), axis_offset);
1460            return ev;
1461         }
1462         if (getDataPointRank()==3) {
1463            DataArrayView::ShapeType ev_shape;
1464            if (axis_offset==0) {
1465              int s2=s[2];
1466              ev_shape.push_back(s2);
1467            }
1468            else if (axis_offset==1) {
1469              int s0=s[0];
1470              ev_shape.push_back(s0);
1471            }
1472            Data ev(0.,ev_shape,getFunctionSpace());
1473            ev.typeMatchRight(*this);
1474            m_data->trace(ev.m_data.get(), axis_offset);
1475            return ev;
1476         }
1477         if (getDataPointRank()==4) {
1478            DataArrayView::ShapeType ev_shape;
1479            if (axis_offset==0) {
1480              ev_shape.push_back(s[2]);
1481              ev_shape.push_back(s[3]);
1482            }
1483            else if (axis_offset==1) {
1484              ev_shape.push_back(s[0]);
1485              ev_shape.push_back(s[3]);
1486            }
1487        else if (axis_offset==2) {
1488          ev_shape.push_back(s[0]);
1489          ev_shape.push_back(s[1]);
1490        }
1491            Data ev(0.,ev_shape,getFunctionSpace());
1492            ev.typeMatchRight(*this);
1493        m_data->trace(ev.m_data.get(), axis_offset);
1494            return ev;
1495         }
1496         else {
1497            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1498         }
1499  }  }
1500    
1501  Data  Data
1502  Data::transpose(int axis) const  Data::transpose(int axis_offset) const
1503    {
1504         #if defined DOPROF
1505         profData->unary++;
1506         #endif
1507         DataArrayView::ShapeType s=getDataPointShape();
1508         DataArrayView::ShapeType ev_shape;
1509         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1510         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1511         int rank=getDataPointRank();
1512         if (axis_offset<0 || axis_offset>rank) {
1513            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1514         }
1515         for (int i=0; i<rank; i++) {
1516           int index = (axis_offset+i)%rank;
1517           ev_shape.push_back(s[index]); // Append to new shape
1518         }
1519         Data ev(0.,ev_shape,getFunctionSpace());
1520         ev.typeMatchRight(*this);
1521         m_data->transpose(ev.m_data.get(), axis_offset);
1522         return ev;
1523    }
1524    
1525    Data
1526    Data::eigenvalues() const
1527    {
1528         #if defined DOPROF
1529            profData->unary++;
1530         #endif
1531         // check input
1532         DataArrayView::ShapeType s=getDataPointShape();
1533         if (getDataPointRank()!=2)
1534            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1535         if(s[0] != s[1])
1536            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1537         // create return
1538         DataArrayView::ShapeType ev_shape(1,s[0]);
1539         Data ev(0.,ev_shape,getFunctionSpace());
1540         ev.typeMatchRight(*this);
1541         m_data->eigenvalues(ev.m_data.get());
1542         return ev;
1543    }
1544    
1545    const boost::python::tuple
1546    Data::eigenvalues_and_eigenvectors(const double tol) const
1547  {  {
1548  #if defined DOPROF       #if defined DOPROF
1549    profData->reduction2++;          profData->unary++;
1550  #endif       #endif
1551    // not implemented       DataArrayView::ShapeType s=getDataPointShape();
1552    throw DataException("Error - Data::transpose not implemented yet.");       if (getDataPointRank()!=2)
1553    return Data();          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1554         if(s[0] != s[1])
1555            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1556         // create return
1557         DataArrayView::ShapeType ev_shape(1,s[0]);
1558         Data ev(0.,ev_shape,getFunctionSpace());
1559         ev.typeMatchRight(*this);
1560         DataArrayView::ShapeType V_shape(2,s[0]);
1561         Data V(0.,V_shape,getFunctionSpace());
1562         V.typeMatchRight(*this);
1563         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1564         return make_tuple(boost::python::object(ev),boost::python::object(V));
1565  }  }
1566    
1567  const boost::python::tuple  const boost::python::tuple
1568  Data::mindp() const  Data::minGlobalDataPoint() const
1569  {  {
1570    // 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
1571    // 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
1572    // surrounding function    // surrounding function
1573    
   int SampleNo;  
1574    int DataPointNo;    int DataPointNo;
1575      int ProcNo;
1576    calc_mindp(SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1577      return make_tuple(ProcNo,DataPointNo);
   return make_tuple(SampleNo,DataPointNo);  
1578  }  }
1579    
1580  void  void
1581  Data::calc_mindp(int& SampleNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1582                   int& DataPointNo) const                          int& DataPointNo) const
1583  {  {
1584    int i,j;    int i,j;
1585    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1248  Data::calc_mindp(int& SampleNo, Line 1615  Data::calc_mindp(int& SampleNo,
1615      }      }
1616    }    }
1617    
1618    SampleNo = lowi;  #ifdef PASO_MPI
1619    DataPointNo = lowj;      // determine the processor on which the minimum occurs
1620        next = temp.getDataPoint(lowi,lowj)();
1621        int lowProc = 0;
1622        double *globalMins = new double[get_MPISize()+1];
1623        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1624        
1625        if( get_MPIRank()==0 ){
1626            next = globalMins[lowProc];
1627            for( i=1; i<get_MPISize(); i++ )
1628                if( next>globalMins[i] ){
1629                    lowProc = i;
1630                    next = globalMins[i];
1631                }
1632        }
1633        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1634    
1635        delete [] globalMins;
1636        ProcNo = lowProc;
1637    #else
1638        ProcNo = 0;
1639    #endif
1640      DataPointNo = lowj + lowi * numDPPSample;
1641  }  }
1642    
1643  void  void
# Line 1273  Data::saveVTK(std::string fileName) cons Line 1661  Data::saveVTK(std::string fileName) cons
1661  Data&  Data&
1662  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1663  {  {
1664  #if defined DOPROF    if (isProtected()) {
1665    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1666  #endif    }
1667    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1668    return (*this);    return (*this);
1669  }  }
# Line 1283  Data::operator+=(const Data& right) Line 1671  Data::operator+=(const Data& right)
1671  Data&  Data&
1672  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1673  {  {
1674  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1675    profData->binary++;    binaryOp(tmp,plus<double>());
 #endif  
   binaryOp(right,plus<double>());  
1676    return (*this);    return (*this);
1677  }  }
1678    
1679  Data&  Data&
1680  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1681  {  {
1682  #if defined DOPROF    if (isProtected()) {
1683    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1684  #endif    }
1685    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1686    return (*this);    return (*this);
1687  }  }
# Line 1303  Data::operator-=(const Data& right) Line 1689  Data::operator-=(const Data& right)
1689  Data&  Data&
1690  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1691  {  {
1692  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1693    profData->binary++;    binaryOp(tmp,minus<double>());
 #endif  
   binaryOp(right,minus<double>());  
1694    return (*this);    return (*this);
1695  }  }
1696    
1697  Data&  Data&
1698  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1699  {  {
1700  #if defined DOPROF    if (isProtected()) {
1701    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1702  #endif    }
1703    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1704    return (*this);    return (*this);
1705  }  }
# Line 1323  Data::operator*=(const Data& right) Line 1707  Data::operator*=(const Data& right)
1707  Data&  Data&
1708  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1709  {  {
1710  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1711    profData->binary++;    binaryOp(tmp,multiplies<double>());
 #endif  
   binaryOp(right,multiplies<double>());  
1712    return (*this);    return (*this);
1713  }  }
1714    
1715  Data&  Data&
1716  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1717  {  {
1718  #if defined DOPROF    if (isProtected()) {
1719    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1720  #endif    }
1721    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1722    return (*this);    return (*this);
1723  }  }
# Line 1343  Data::operator/=(const Data& right) Line 1725  Data::operator/=(const Data& right)
1725  Data&  Data&
1726  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1727  {  {
1728  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1729    profData->binary++;    binaryOp(tmp,divides<double>());
 #endif  
   binaryOp(right,divides<double>());  
1730    return (*this);    return (*this);
1731  }  }
1732    
1733  Data  Data
1734    Data::rpowO(const boost::python::object& left) const
1735    {
1736      Data left_d(left,*this);
1737      return left_d.powD(*this);
1738    }
1739    
1740    Data
1741  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1742  {  {
1743  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1744    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1745  }  }
1746    
1747  Data  Data
1748  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1749  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1750    Data result;    Data result;
1751    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1752    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1753         result.binaryOp(*this,escript::rpow);
1754      } else {
1755         result.copy(*this);
1756         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1757      }
1758    return result;    return result;
1759  }  }
1760    
1761    
1762  //  //
1763  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1764  Data  Data
# Line 1382  escript::operator+(const Data& left, con Line 1767  escript::operator+(const Data& left, con
1767    Data result;    Data result;
1768    //    //
1769    // perform a deep copy    // perform a deep copy
1770    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1771    result+=right;       result.copy(right);
1772         result+=left;
1773      } else {
1774         result.copy(left);
1775         result+=right;
1776      }
1777    return result;    return result;
1778  }  }
1779    
# Line 1395  escript::operator-(const Data& left, con Line 1785  escript::operator-(const Data& left, con
1785    Data result;    Data result;
1786    //    //
1787    // perform a deep copy    // perform a deep copy
1788    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1789    result-=right;       result=right.neg();
1790         result+=left;
1791      } else {
1792         result.copy(left);
1793         result-=right;
1794      }
1795    return result;    return result;
1796  }  }
1797    
# Line 1408  escript::operator*(const Data& left, con Line 1803  escript::operator*(const Data& left, con
1803    Data result;    Data result;
1804    //    //
1805    // perform a deep copy    // perform a deep copy
1806    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1807    result*=right;       result.copy(right);
1808         result*=left;
1809      } else {
1810         result.copy(left);
1811         result*=right;
1812      }
1813    return result;    return result;
1814  }  }
1815    
# Line 1421  escript::operator/(const Data& left, con Line 1821  escript::operator/(const Data& left, con
1821    Data result;    Data result;
1822    //    //
1823    // perform a deep copy    // perform a deep copy
1824    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1825    result/=right;       result=right.oneOver();
1826         result*=left;
1827      } else {
1828         result.copy(left);
1829         result/=right;
1830      }
1831    return result;    return result;
1832  }  }
1833    
# Line 1431  escript::operator/(const Data& left, con Line 1836  escript::operator/(const Data& left, con
1836  Data  Data
1837  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1838  {  {
1839    //    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;  
1840  }  }
1841    
1842  //  //
# Line 1447  escript::operator+(const Data& left, con Line 1844  escript::operator+(const Data& left, con
1844  Data  Data
1845  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1846  {  {
1847    //    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;  
1848  }  }
1849    
1850  //  //
# Line 1463  escript::operator-(const Data& left, con Line 1852  escript::operator-(const Data& left, con
1852  Data  Data
1853  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1854  {  {
1855    //    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;  
1856  }  }
1857    
1858  //  //
# Line 1479  escript::operator*(const Data& left, con Line 1860  escript::operator*(const Data& left, con
1860  Data  Data
1861  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1862  {  {
1863    //    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;  
1864  }  }
1865    
1866  //  //
# Line 1495  escript::operator/(const Data& left, con Line 1868  escript::operator/(const Data& left, con
1868  Data  Data
1869  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1870  {  {
1871    //    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;  
1872  }  }
1873    
1874  //  //
# Line 1508  escript::operator+(const boost::python:: Line 1876  escript::operator+(const boost::python::
1876  Data  Data
1877  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1878  {  {
1879    //    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;  
1880  }  }
1881    
1882  //  //
# Line 1521  escript::operator-(const boost::python:: Line 1884  escript::operator-(const boost::python::
1884  Data  Data
1885  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1886  {  {
1887    //    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;  
1888  }  }
1889    
1890  //  //
# Line 1534  escript::operator*(const boost::python:: Line 1892  escript::operator*(const boost::python::
1892  Data  Data
1893  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1894  {  {
1895    //    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;  
1896  }  }
1897    
1898  //  //
# Line 1587  escript::operator/(const boost::python:: Line 1940  escript::operator/(const boost::python::
1940  //  return ret;  //  return ret;
1941  //}  //}
1942    
1943    /* TODO */
1944    /* global reduction */
1945  Data  Data
1946  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1947  {  {
# Line 1601  Data::getItem(const boost::python::objec Line 1956  Data::getItem(const boost::python::objec
1956    return getSlice(slice_region);    return getSlice(slice_region);
1957  }  }
1958    
1959    /* TODO */
1960    /* global reduction */
1961  Data  Data
1962  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1963  {  {
# Line 1610  Data::getSlice(const DataArrayView::Regi Line 1967  Data::getSlice(const DataArrayView::Regi
1967    return Data(*this,region);    return Data(*this,region);
1968  }  }
1969    
1970    /* TODO */
1971    /* global reduction */
1972  void  void
1973  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1974                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1618  Data::setItemO(const boost::python::obje Line 1977  Data::setItemO(const boost::python::obje
1977    setItemD(key,tempData);    setItemD(key,tempData);
1978  }  }
1979    
1980    /* TODO */
1981    /* global reduction */
1982  void  void
1983  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1984                 const Data& value)                 const Data& value)
# Line 1635  Data::setItemD(const boost::python::obje Line 1996  Data::setItemD(const boost::python::obje
1996    }    }
1997  }  }
1998    
1999    /* TODO */
2000    /* global reduction */
2001  void  void
2002  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2003                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
2004  {  {
2005      if (isProtected()) {
2006            throw DataException("Error - attempt to update protected Data object.");
2007      }
2008  #if defined DOPROF  #if defined DOPROF
2009    profData->slicing++;    profData->slicing++;
2010  #endif  #endif
# Line 1676  Data::typeMatchRight(const Data& right) Line 2042  Data::typeMatchRight(const Data& right)
2042    }    }
2043  }  }
2044    
2045    /* TODO */
2046    /* global reduction */
2047  void  void
2048  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2049                       const boost::python::object& value)                       const boost::python::object& value)
2050  {  {
2051      if (isProtected()) {
2052            throw DataException("Error - attempt to update protected Data object.");
2053      }
2054    //    //
2055    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2056    tag();    tag();
# Line 1697  Data::setTaggedValue(int tagKey, Line 2068  Data::setTaggedValue(int tagKey,
2068    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2069  }  }
2070    
2071    /* TODO */
2072    /* global reduction */
2073  void  void
2074  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2075                              const DataArrayView& value)                              const DataArrayView& value)
2076  {  {
2077      if (isProtected()) {
2078            throw DataException("Error - attempt to update protected Data object.");
2079      }
2080    //    //
2081    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2082    tag();    tag();
# Line 1714  Data::setTaggedValueFromCPP(int tagKey, Line 2090  Data::setTaggedValueFromCPP(int tagKey,
2090    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
2091  }  }
2092    
2093    /* TODO */
2094    /* global reduction */
2095  int  int
2096  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2097  {  {
# Line 1721  Data::getTagNumber(int dpno) Line 2099  Data::getTagNumber(int dpno)
2099  }  }
2100    
2101  void  void
 Data::setRefValue(int ref,  
                   const boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
   
   //  
   // Load DataArray with values from data-points specified by ref  
   m_data->getRefValue(ref,valueDataArray);  
   
   //  
   // Load values from valueDataArray into return numarray  
   
   // extract the shape of the numarray  
   int rank = value.getrank();  
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
   
   if (rank==0) {  
       boost::python::numeric::array temp_numArray(valueView());  
       value = temp_numArray;  
   }  
   if (rank==1) {  
     for (int i=0; i < shape[0]; i++) {  
       value[i] = valueView(i);  
     }  
   }  
   if (rank==2) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         value[i][j] = valueView(i,j);  
       }  
     }  
   }  
   if (rank==3) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           value[i][j][k] = valueView(i,j,k);  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           for (int l=0; l < shape[3]; l++) {  
             value[i][j][k][l] = valueView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
 void  
2102  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
2103  {  {
2104    cout << "Archiving Data object to: " << fileName << endl;    cout << "Archiving Data object to: " << fileName << endl;
# Line 1836  Data::archiveData(const std::string file Line 2137  Data::archiveData(const std::string file
2137    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2138    int dataLength = getLength();    int dataLength = getLength();
2139    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2140    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2141    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2142      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
2143    }    }
2144    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2145    if (isTagged()) {    if (isTagged()) {
2146      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2147        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 1994  Data::extractData(const std::string file Line 2295  Data::extractData(const std::string file
2295        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2296      }      }
2297    }    }
2298    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2299    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2300      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2301    }    }
2302    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2303    if (dataType==2) {    if (dataType==2) {
2304      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2305        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 2047  Data::extractData(const std::string file Line 2348  Data::extractData(const std::string file
2348      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2349    }    }
2350    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2351      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2352        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2353      }      }
2354    }    }
# Line 2137  ostream& escript::operator<<(ostream& o, Line 2438  ostream& escript::operator<<(ostream& o,
2438    o << data.toString();    o << data.toString();
2439    return o;    return o;
2440  }  }
2441    
2442    Data
2443    escript::C_GeneralTensorProduct(Data& arg_0,
2444                         Data& arg_1,
2445                         int axis_offset,
2446                         int transpose)
2447    {
2448      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2449      // SM is the product of the last axis_offset entries in arg_0.getShape().
2450    
2451      #if defined DOPROF
2452        // profData->binary++;
2453      #endif
2454    
2455      // Interpolate if necessary and find an appropriate function space
2456      Data arg_0_Z, arg_1_Z;
2457      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2458        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2459          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2460          arg_1_Z = Data(arg_1);
2461        }
2462        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2463          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2464          arg_0_Z =Data(arg_0);
2465        }
2466        else {
2467          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2468        }
2469      } else {
2470          arg_0_Z = Data(arg_0);
2471          arg_1_Z = Data(arg_1);
2472      }
2473      // Get rank and shape of inputs
2474      int rank0 = arg_0_Z.getDataPointRank();
2475      int rank1 = arg_1_Z.getDataPointRank();
2476      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2477      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2478    
2479      // Prepare for the loops of the product and verify compatibility of shapes
2480      int start0=0, start1=0;
2481      if (transpose == 0)       {}
2482      else if (transpose == 1)  { start0 = axis_offset; }
2483      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2484      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2485    
2486      // Adjust the shapes for transpose
2487      DataArrayView::ShapeType tmpShape0;
2488      DataArrayView::ShapeType tmpShape1;
2489      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2490      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2491    
2492    #if 0
2493      // For debugging: show shape after transpose
2494      char tmp[100];
2495      std::string shapeStr;
2496      shapeStr = "(";
2497      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2498      shapeStr += ")";
2499      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2500      shapeStr = "(";
2501      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2502      shapeStr += ")";
2503      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2504    #endif
2505    
2506      // Prepare for the loops of the product
2507      int SL=1, SM=1, SR=1;
2508      for (int i=0; i<rank0-axis_offset; i++)   {
2509        SL *= tmpShape0[i];
2510      }
2511      for (int i=rank0-axis_offset; i<rank0; i++)   {
2512        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2513          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2514        }
2515        SM *= tmpShape0[i];
2516      }
2517      for (int i=axis_offset; i<rank1; i++)     {
2518        SR *= tmpShape1[i];
2519      }
2520    
2521      // Define the shape of the output
2522      DataArrayView::ShapeType shape2;
2523      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2524      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2525    
2526      // Declare output Data object
2527      Data res;
2528    
2529      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2530        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2531        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2532        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2533        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2534        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2535      }
2536      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2537    
2538        // Prepare the DataConstant input
2539        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2540        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2541    
2542        // Borrow DataTagged input from Data object
2543        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2544        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2545    
2546        // Prepare a DataTagged output 2
2547        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2548        res.tag();
2549        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2550        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2551    
2552        // Prepare offset into DataConstant
2553        int offset_0 = tmp_0->getPointOffset(0,0);
2554        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2555        // Get the views
2556        DataArrayView view_1 = tmp_1->getDefaultValue();
2557        DataArrayView view_2 = tmp_2->getDefaultValue();
2558        // Get the pointers to the actual data
2559        double *ptr_1 = &((view_1.getData())[0]);
2560        double *ptr_2 = &((view_2.getData())[0]);
2561        // Compute an MVP for the default
2562        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2563        // Compute an MVP for each tag
2564        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2565        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2566        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2567          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2568          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2569          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2570          double *ptr_1 = &view_1.getData(0);
2571          double *ptr_2 = &view_2.getData(0);
2572          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2573        }
2574    
2575      }
2576      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2577    
2578        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2579        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2580        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2581        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2582        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2583        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2584        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2585        int sampleNo_1,dataPointNo_1;
2586        int numSamples_1 = arg_1_Z.getNumSamples();
2587        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2588        int offset_0 = tmp_0->getPointOffset(0,0);
2589        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2590        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2591          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2592            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2593            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2594            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2595            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2596            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2597            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2598          }
2599        }
2600    
2601      }
2602      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2603    
2604        // Borrow DataTagged input from Data object
2605        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2606        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2607    
2608        // Prepare the DataConstant input
2609        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2610        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2611    
2612        // Prepare a DataTagged output 2
2613        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2614        res.tag();
2615        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2616        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2617    
2618        // Prepare offset into DataConstant
2619        int offset_1 = tmp_1->getPointOffset(0,0);
2620        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2621        // Get the views
2622        DataArrayView view_0 = tmp_0->getDefaultValue();
2623        DataArrayView view_2 = tmp_2->getDefaultValue();
2624        // Get the pointers to the actual data
2625        double *ptr_0 = &((view_0.getData())[0]);
2626        double *ptr_2 = &((view_2.getData())[0]);
2627        // Compute an MVP for the default
2628        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2629        // Compute an MVP for each tag
2630        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2631        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2632        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2633          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2634          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2635          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2636          double *ptr_0 = &view_0.getData(0);
2637          double *ptr_2 = &view_2.getData(0);
2638          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2639        }
2640    
2641      }
2642      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2643    
2644        // Borrow DataTagged input from Data object
2645        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2646        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2647    
2648        // Borrow DataTagged input from Data object
2649        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2650        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2651    
2652        // Prepare a DataTagged output 2
2653        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2654        res.tag();  // DataTagged output
2655        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2656        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2657    
2658        // Get the views
2659        DataArrayView view_0 = tmp_0->getDefaultValue();
2660        DataArrayView view_1 = tmp_1->getDefaultValue();
2661        DataArrayView view_2 = tmp_2->getDefaultValue();
2662        // Get the pointers to the actual data
2663        double *ptr_0 = &((view_0.getData())[0]);
2664        double *ptr_1 = &((view_1.getData())[0]);
2665        double *ptr_2 = &((view_2.getData())[0]);
2666        // Compute an MVP for the default
2667        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2668        // Merge the tags
2669        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2670        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2671        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2672        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2673          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2674        }
2675        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2676          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2677        }
2678        // Compute an MVP for each tag
2679        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2680        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2681          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2682          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2683          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2684          double *ptr_0 = &view_0.getData(0);
2685          double *ptr_1 = &view_1.getData(0);
2686          double *ptr_2 = &view_2.getData(0);
2687          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2688        }
2689    
2690      }
2691      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2692    
2693        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2694        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2695        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2696        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2697        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2698        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2699        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2700        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2701        int sampleNo_0,dataPointNo_0;
2702        int numSamples_0 = arg_0_Z.getNumSamples();
2703        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2704        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2705        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2706          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2707          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2708          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2709            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2710            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2711            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2712            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2713            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2714          }
2715        }
2716    
2717      }
2718      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2719    
2720        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2721        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2722        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2723        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2724        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2725        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2726        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2727        int sampleNo_0,dataPointNo_0;
2728        int numSamples_0 = arg_0_Z.getNumSamples();
2729        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2730        int offset_1 = tmp_1->getPointOffset(0,0);
2731        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2732        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2733          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2734            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2735            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2736            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2737            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2738            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2739            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2740          }
2741        }
2742    
2743    
2744      }
2745      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2746    
2747        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2748        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2749        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2750        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2751        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2752        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2753        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2754        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2755        int sampleNo_0,dataPointNo_0;
2756        int numSamples_0 = arg_0_Z.getNumSamples();
2757        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2758        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2759        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2760          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2761          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2762          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2763            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2764            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2765            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2766            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2767            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2768          }
2769        }
2770    
2771      }
2772      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2773    
2774        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2775        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2776        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2777        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2778        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2779        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2780        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2781        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2782        int sampleNo_0,dataPointNo_0;
2783        int numSamples_0 = arg_0_Z.getNumSamples();
2784        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2785        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2786        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2787          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2788            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2789            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2790            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2791            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2792            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2793            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2794            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2795          }
2796        }
2797    
2798      }
2799      else {
2800        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2801      }
2802    
2803      return res;
2804    }
2805    
2806    DataAbstract*
2807    Data::borrowData() const
2808    {
2809      return m_data.get();
2810    }
2811    
2812    /* Member functions specific to the MPI implementation */
2813    
2814    void
2815    Data::print()
2816    {
2817      int i,j;
2818      
2819      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2820      for( i=0; i<getNumSamples(); i++ )
2821      {
2822        printf( "[%6d]", i );
2823        for( j=0; j<getNumDataPointsPerSample(); j++ )
2824          printf( "\t%10.7g", (getSampleData(i))[j] );
2825        printf( "\n" );
2826      }
2827    }
2828    
2829    int
2830    Data::get_MPISize() const
2831    {
2832        int error, size;
2833    #ifdef PASO_MPI
2834        error = MPI_Comm_size( get_MPIComm(), &size );
2835    #else
2836        size = 1;
2837    #endif
2838        return size;
2839    }
2840    
2841    int
2842    Data::get_MPIRank() const
2843    {
2844        int error, rank;
2845    #ifdef PASO_MPI
2846        error = MPI_Comm_rank( get_MPIComm(), &rank );
2847    #else
2848        rank = 0;
2849    #endif
2850        return rank;
2851    }
2852    
2853    MPI_Comm
2854    Data::get_MPIComm() const
2855    {
2856    #ifdef PASO_MPI
2857        return MPI_COMM_WORLD;
2858    #else
2859        return -1;
2860    #endif
2861    }
2862    

Legend:
Removed from v.474  
changed lines
  Added in v.1031

  ViewVC Help
Powered by ViewVC 1.1.26