/[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

trunk/esys2/escript/src/Data/Data.cpp revision 153 by jgs, Tue Oct 25 01:51:20 2005 UTC trunk/escript/src/Data.cpp revision 854 by gross, Thu Sep 21 05:29:42 2006 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
2    
3   ******************************************************************************  /*
4   *                                                                            *   ************************************************************
5   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *   *          Copyright 2006 by ACcESS MNRF                   *
6   *                                                                            *   *                                                          *
7   * This software is the property of ACcESS.  No part of this code             *   *              http://www.access.edu.au                    *
8   * may be copied in any form or by any means without the expressed written    *   *       Primary Business: Queensland, Australia            *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *  Licensed under the Open Software License version 3.0    *
10   * by any unauthorised person is illegal unless that                          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
11   * person has a software license agreement with ACcESS.                       *   *                                                          *
12   *                                                                            *   ************************************************************
13   ******************************************************************************  */
14    #include "Data.h"
15    
16    #include "DataExpanded.h"
17    #include "DataConstant.h"
18    #include "DataTagged.h"
19    #include "DataEmpty.h"
20    #include "DataArray.h"
21    #include "DataArrayView.h"
22    #include "DataProf.h"
23    #include "FunctionSpaceFactory.h"
24    #include "AbstractContinuousDomain.h"
25    #include "UnaryFuncs.h"
26    
 ******************************************************************************/  
   
 #include "escript/Data/Data.h"  
   
 #include <iostream>  
27  #include <fstream>  #include <fstream>
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
 #include <exception>  
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
 #include <boost/python/str.hpp>  
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
34  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
35    
 #include "escript/Data/DataException.h"  
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataEmpty.h"  
 #include "escript/Data/DataArray.h"  
 #include "escript/Data/DataProf.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/AbstractContinuousDomain.h"  
 #include "escript/Data/UnaryFuncs.h"  
   
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
# Line 59  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 76  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 90  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 99  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 113  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 122  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 138  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 153  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 167  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 178  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 190  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 214  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 344  Data::isTagged() const Line 350  Data::isTagged() const
350    return (temp!=0);    return (temp!=0);
351  }  }
352    
353    /* TODO */
354    /* global reduction -- the local data being empty does not imply that it is empty on other processers*/
355  bool  bool
356  Data::isEmpty() const  Data::isEmpty() const
357  {  {
# Line 359  Data::isConstant() const Line 367  Data::isConstant() const
367  }  }
368    
369  void  void
370    Data::setProtection()
371    {
372       m_protected=true;
373    }
374    
375    bool
376    Data::isProtected() const
377    {
378       return m_protected;
379    }
380    
381    
382    
383    void
384  Data::expand()  Data::expand()
385  {  {
386    if (isConstant()) {    if (isConstant()) {
# Line 400  Data::tag() Line 422  Data::tag()
422    }    }
423  }  }
424    
425  void  Data
426  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
427  {  {
428    m_data->reshapeDataPoint(shape);  #if defined DOPROF
429      profData->where++;
430    #endif
431      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
432  }  }
433    
434  Data  Data
# Line 443  Data::whereNonPositive() const Line 468  Data::whereNonPositive() const
468  }  }
469    
470  Data  Data
471  Data::whereZero() const  Data::whereZero(double tol) const
472  {  {
473  #if defined DOPROF  #if defined DOPROF
474    profData->where++;    profData->where++;
475  #endif  #endif
476    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
477      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
478  }  }
479    
480  Data  Data
481  Data::whereNonZero() const  Data::whereNonZero(double tol) const
482  {  {
483  #if defined DOPROF  #if defined DOPROF
484    profData->where++;    profData->where++;
485  #endif  #endif
486    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
487      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
488  }  }
489    
490  Data  Data
# Line 526  Data::getDataPointShape() const Line 553  Data::getDataPointShape() const
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 755  Data::convertToNumArrayFromSampleNo(int Line 785  Data::convertToNumArrayFromSampleNo(int
785    
786  const  const
787  boost::python::numeric::array  boost::python::numeric::array
788  Data::convertToNumArrayFromDPNo(int sampleNo,  Data::convertToNumArrayFromDPNo(int procNo,
789                                  int dataPointNo)                                  int sampleNo,
790                                                                    int dataPointNo)
791    
792  {  {
793        size_t length=0;
794        int i, j, k, l, pos;
795    
796    //    //
797    // Check a valid sample number has been supplied    // Check a valid sample number has been supplied
798    if (sampleNo >= getNumSamples()) {    if (sampleNo >= getNumSamples()) {
# Line 803  Data::convertToNumArrayFromDPNo(int samp Line 838  Data::convertToNumArrayFromDPNo(int samp
838      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
839    }    }
840    
841        // added for the MPI communication
842        length=1;
843        for( i=0; i<arrayRank; i++ )
844            length *= arrayShape[i];
845        double *tmpData = new double[length];
846    
847    //    //
848    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
849    DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
850        // updated for the MPI case
851        if( get_MPIRank()==procNo ){
852            // create a view of the data if it is stored locally
853            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
854            
855            // pack the data from the view into tmpData for MPI communication
856            pos=0;
857            switch( dataPointRank ){
858                case 0 :
859                    tmpData[0] = dataPointView();
860                    break;
861                case 1 :        
862                    for( i=0; i<dataPointShape[0]; i++ )
863                        tmpData[i]=dataPointView(i);
864                    break;
865                case 2 :        
866                    for( i=0; i<dataPointShape[0]; i++ )
867                        for( j=0; j<dataPointShape[1]; j++, pos++ )
868                            tmpData[pos]=dataPointView(i,j);
869                    break;
870                case 3 :        
871                    for( i=0; i<dataPointShape[0]; i++ )
872                        for( j=0; j<dataPointShape[1]; j++ )
873                            for( k=0; k<dataPointShape[2]; k++, pos++ )
874                                tmpData[pos]=dataPointView(i,j,k);
875                    break;
876                case 4 :
877                    for( i=0; i<dataPointShape[0]; i++ )
878                        for( j=0; j<dataPointShape[1]; j++ )
879                            for( k=0; k<dataPointShape[2]; k++ )
880                                for( l=0; l<dataPointShape[3]; l++, pos++ )
881                                    tmpData[pos]=dataPointView(i,j,k,l);
882                    break;
883            }
884        }
885    #ifdef PASO_MPI
886            // broadcast the data to all other processes
887            MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
888    #endif
889    
890        // unpack the data
891        switch( dataPointRank ){
892            case 0 :
893                numArray[i]=tmpData[0];
894                break;
895            case 1 :        
896                for( i=0; i<dataPointShape[0]; i++ )
897                    numArray[i]=tmpData[i];
898                break;
899            case 2 :        
900                for( i=0; i<dataPointShape[0]; i++ )
901                    for( j=0; j<dataPointShape[1]; j++ )
902                        tmpData[i+j*dataPointShape[0]];
903                break;
904            case 3 :        
905                for( i=0; i<dataPointShape[0]; i++ )
906                    for( j=0; j<dataPointShape[1]; j++ )
907                        for( k=0; k<dataPointShape[2]; k++ )
908                            tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
909                break;
910            case 4 :
911                for( i=0; i<dataPointShape[0]; i++ )
912                    for( j=0; j<dataPointShape[1]; j++ )
913                        for( k=0; k<dataPointShape[2]; k++ )
914                            for( l=0; l<dataPointShape[3]; l++ )
915                                tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
916                break;
917        }
918    
919        delete [] tmpData;  
920    /*
921    if (dataPointRank==0) {    if (dataPointRank==0) {
922      numArray[0]=dataPointView();      numArray[0]=dataPointView();
923    }    }
# Line 841  Data::convertToNumArrayFromDPNo(int samp Line 953  Data::convertToNumArrayFromDPNo(int samp
953        }        }
954      }      }
955    }    }
956    */
957    
958    //    //
959    // return the loaded array    // return the loaded array
# Line 880  Data::integrate() const Line 993  Data::integrate() const
993      }      }
994    }    }
995    if (rank==2) {    if (rank==2) {
996      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
997      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
998        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
999          index = i + shape[0] * j;             index = i + shape[0] * j;
1000          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
1001        }           }
1002      }         }
1003    }    }
1004    if (rank==3) {    if (rank==3) {
1005      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 894  Data::integrate() const Line 1007  Data::integrate() const
1007        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
1008          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1009            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
1010            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
1011          }          }
1012        }        }
1013      }      }
# Line 906  Data::integrate() const Line 1019  Data::integrate() const
1019          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1020            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
1021              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1022              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
1023            }            }
1024          }          }
1025        }        }
# Line 1027  Data::atanh() const Line 1140  Data::atanh() const
1140  }  }
1141    
1142  Data  Data
1143  Data::log() const  Data::log10() const
1144  {  {
1145  #if defined DOPROF  #if defined DOPROF
1146    profData->unary++;    profData->unary++;
# Line 1036  Data::log() const Line 1149  Data::log() const
1149  }  }
1150    
1151  Data  Data
1152  Data::ln() const  Data::log() const
1153  {  {
1154  #if defined DOPROF  #if defined DOPROF
1155    profData->unary++;    profData->unary++;
# Line 1104  Data::sqrt() const Line 1217  Data::sqrt() const
1217  double  double
1218  Data::Lsup() const  Data::Lsup() const
1219  {  {
1220      double localValue, globalValue;
1221  #if defined DOPROF  #if defined DOPROF
1222    profData->reduction1++;    profData->reduction1++;
1223  #endif  #endif
1224    //    //
1225    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1226    
1227    AbsMax abs_max_func;    AbsMax abs_max_func;
1228    return algorithm(abs_max_func,0);    localValue = algorithm(abs_max_func,0);
1229    #ifdef PASO_MPI
1230      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1231      return globalValue;
1232    #else
1233      return localValue;
1234    #endif
1235  }  }
1236    
1237  double  double
1238  Data::Linf() const  Data::Linf() const
1239  {  {
1240      double localValue, globalValue;
1241  #if defined DOPROF  #if defined DOPROF
1242    profData->reduction1++;    profData->reduction1++;
1243  #endif  #endif
1244    //    //
1245    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1246    AbsMin abs_min_func;    AbsMin abs_min_func;
1247    return algorithm(abs_min_func,numeric_limits<double>::max());    localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1248    
1249    #ifdef PASO_MPI
1250      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1251      return globalValue;
1252    #else
1253      return localValue;
1254    #endif
1255  }  }
1256    
1257  double  double
1258  Data::sup() const  Data::sup() const
1259  {  {
1260      double localValue, globalValue;
1261  #if defined DOPROF  #if defined DOPROF
1262    profData->reduction1++;    profData->reduction1++;
1263  #endif  #endif
1264    //    //
1265    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1266    FMax fmax_func;    FMax fmax_func;
1267    return algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1268    #ifdef PASO_MPI
1269      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1270      return globalValue;
1271    #else
1272      return localValue;
1273    #endif
1274  }  }
1275    
1276  double  double
1277  Data::inf() const  Data::inf() const
1278  {  {
1279      double localValue, globalValue;
1280  #if defined DOPROF  #if defined DOPROF
1281    profData->reduction1++;    profData->reduction1++;
1282  #endif  #endif
1283    //    //
1284    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1285    FMin fmin_func;    FMin fmin_func;
1286    return algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1287    #ifdef PASO_MPI
1288      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1289      return globalValue;
1290    #else
1291      return localValue;
1292    #endif
1293  }  }
1294    
1295    /* TODO */
1296    /* global reduction */
1297  Data  Data
1298  Data::maxval() const  Data::maxval() const
1299  {  {
# Line 1174  Data::minval() const Line 1319  Data::minval() const
1319  }  }
1320    
1321  Data  Data
1322  Data::length() const  Data::swapaxes(const int axis0, const int axis1) const
1323  {  {
1324  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1325    profData->reduction2++;       #if defined DOPROF
1326  #endif       profData->unary++;
1327    Length len_func;       #endif
1328    return dp_algorithm(len_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1329         DataArrayView::ShapeType ev_shape;
1330         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1331         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1332         int rank=getDataPointRank();
1333         if (rank<2) {
1334            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1335         }
1336         if (axis0<0 || axis0>rank-1) {
1337            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1338         }
1339         if (axis1<0 || axis1>rank-1) {
1340             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1341         }
1342         if (axis0 == axis1) {
1343             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1344         }
1345         if (axis0 > axis1) {
1346             axis0_tmp=axis1;
1347             axis1_tmp=axis0;
1348         } else {
1349             axis0_tmp=axis0;
1350             axis1_tmp=axis1;
1351         }
1352         for (int i=0; i<rank; i++) {
1353           if (i == axis0_tmp) {
1354              ev_shape.push_back(s[axis1_tmp]);
1355           } else if (i == axis1_tmp) {
1356              ev_shape.push_back(s[axis0_tmp]);
1357           } else {
1358              ev_shape.push_back(s[i]);
1359           }
1360         }
1361         Data ev(0.,ev_shape,getFunctionSpace());
1362         ev.typeMatchRight(*this);
1363         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1364         return ev;
1365    
1366    }
1367    
1368    Data
1369    Data::symmetric() const
1370    {
1371         #if defined DOPROF
1372            profData->unary++;
1373         #endif
1374         // check input
1375         DataArrayView::ShapeType s=getDataPointShape();
1376         if (getDataPointRank()==2) {
1377            if(s[0] != s[1])
1378               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1379         }
1380         else if (getDataPointRank()==4) {
1381            if(!(s[0] == s[2] && s[1] == s[3]))
1382               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1383         }
1384         else {
1385            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1386         }
1387         Data ev(0.,getDataPointShape(),getFunctionSpace());
1388         ev.typeMatchRight(*this);
1389         m_data->symmetric(ev.m_data.get());
1390         return ev;
1391    }
1392    
1393    Data
1394    Data::nonsymmetric() const
1395    {
1396         #if defined DOPROF
1397            profData->unary++;
1398         #endif
1399         // check input
1400         DataArrayView::ShapeType s=getDataPointShape();
1401         if (getDataPointRank()==2) {
1402            if(s[0] != s[1])
1403               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1404            DataArrayView::ShapeType ev_shape;
1405            ev_shape.push_back(s[0]);
1406            ev_shape.push_back(s[1]);
1407            Data ev(0.,ev_shape,getFunctionSpace());
1408            ev.typeMatchRight(*this);
1409            m_data->nonsymmetric(ev.m_data.get());
1410            return ev;
1411         }
1412         else if (getDataPointRank()==4) {
1413            if(!(s[0] == s[2] && s[1] == s[3]))
1414               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1415            DataArrayView::ShapeType ev_shape;
1416            ev_shape.push_back(s[0]);
1417            ev_shape.push_back(s[1]);
1418            ev_shape.push_back(s[2]);
1419            ev_shape.push_back(s[3]);
1420            Data ev(0.,ev_shape,getFunctionSpace());
1421            ev.typeMatchRight(*this);
1422            m_data->nonsymmetric(ev.m_data.get());
1423            return ev;
1424         }
1425         else {
1426            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1427         }
1428    }
1429    
1430    Data
1431    Data::trace(int axis_offset) const
1432    {
1433         #if defined DOPROF
1434            profData->unary++;
1435         #endif
1436         DataArrayView::ShapeType s=getDataPointShape();
1437         if (getDataPointRank()==2) {
1438            DataArrayView::ShapeType ev_shape;
1439            Data ev(0.,ev_shape,getFunctionSpace());
1440            ev.typeMatchRight(*this);
1441            m_data->trace(ev.m_data.get(), axis_offset);
1442            return ev;
1443         }
1444         if (getDataPointRank()==3) {
1445            DataArrayView::ShapeType ev_shape;
1446            if (axis_offset==0) {
1447              int s2=s[2];
1448              ev_shape.push_back(s2);
1449            }
1450            else if (axis_offset==1) {
1451              int s0=s[0];
1452              ev_shape.push_back(s0);
1453            }
1454            Data ev(0.,ev_shape,getFunctionSpace());
1455            ev.typeMatchRight(*this);
1456            m_data->trace(ev.m_data.get(), axis_offset);
1457            return ev;
1458         }
1459         if (getDataPointRank()==4) {
1460            DataArrayView::ShapeType ev_shape;
1461            if (axis_offset==0) {
1462              ev_shape.push_back(s[2]);
1463              ev_shape.push_back(s[3]);
1464            }
1465            else if (axis_offset==1) {
1466              ev_shape.push_back(s[0]);
1467              ev_shape.push_back(s[3]);
1468            }
1469        else if (axis_offset==2) {
1470          ev_shape.push_back(s[0]);
1471          ev_shape.push_back(s[1]);
1472        }
1473            Data ev(0.,ev_shape,getFunctionSpace());
1474            ev.typeMatchRight(*this);
1475        m_data->trace(ev.m_data.get(), axis_offset);
1476            return ev;
1477         }
1478         else {
1479            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1480         }
1481  }  }
1482    
1483  Data  Data
1484  Data::trace() const  Data::transpose(int axis_offset) const
1485  {  {
1486  #if defined DOPROF       #if defined DOPROF
1487    profData->reduction2++;       profData->unary++;
1488  #endif       #endif
1489    Trace trace_func;       DataArrayView::ShapeType s=getDataPointShape();
1490    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType ev_shape;
1491         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1492         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1493         int rank=getDataPointRank();
1494         if (axis_offset<0 || axis_offset>rank) {
1495            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1496         }
1497         for (int i=0; i<rank; i++) {
1498           int index = (axis_offset+i)%rank;
1499           ev_shape.push_back(s[index]); // Append to new shape
1500         }
1501         Data ev(0.,ev_shape,getFunctionSpace());
1502         ev.typeMatchRight(*this);
1503         m_data->transpose(ev.m_data.get(), axis_offset);
1504         return ev;
1505  }  }
1506    
1507  Data  Data
1508  Data::transpose(int axis) const  Data::eigenvalues() const
1509    {
1510         #if defined DOPROF
1511            profData->unary++;
1512         #endif
1513         // check input
1514         DataArrayView::ShapeType s=getDataPointShape();
1515         if (getDataPointRank()!=2)
1516            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1517         if(s[0] != s[1])
1518            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1519         // create return
1520         DataArrayView::ShapeType ev_shape(1,s[0]);
1521         Data ev(0.,ev_shape,getFunctionSpace());
1522         ev.typeMatchRight(*this);
1523         m_data->eigenvalues(ev.m_data.get());
1524         return ev;
1525    }
1526    
1527    const boost::python::tuple
1528    Data::eigenvalues_and_eigenvectors(const double tol) const
1529  {  {
1530  #if defined DOPROF       #if defined DOPROF
1531    profData->reduction2++;          profData->unary++;
1532  #endif       #endif
1533    // not implemented       DataArrayView::ShapeType s=getDataPointShape();
1534    throw DataException("Error - Data::transpose not implemented yet.");       if (getDataPointRank()!=2)
1535    return Data();          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1536         if(s[0] != s[1])
1537            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1538         // create return
1539         DataArrayView::ShapeType ev_shape(1,s[0]);
1540         Data ev(0.,ev_shape,getFunctionSpace());
1541         ev.typeMatchRight(*this);
1542         DataArrayView::ShapeType V_shape(2,s[0]);
1543         Data V(0.,V_shape,getFunctionSpace());
1544         V.typeMatchRight(*this);
1545         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1546         return make_tuple(boost::python::object(ev),boost::python::object(V));
1547  }  }
1548    
1549  const boost::python::tuple  const boost::python::tuple
# Line 1213  Data::mindp() const Line 1555  Data::mindp() const
1555    
1556    int SampleNo;    int SampleNo;
1557    int DataPointNo;    int DataPointNo;
1558        int ProcNo;
1559    calc_mindp(SampleNo,DataPointNo);    calc_mindp(ProcNo,SampleNo,DataPointNo);
1560      return make_tuple(ProcNo,SampleNo,DataPointNo);
   return make_tuple(SampleNo,DataPointNo);  
1561  }  }
1562    
1563  void  void
1564  Data::calc_mindp(int& SampleNo,  Data::calc_mindp(   int& ProcNo,
1565                   int& DataPointNo) const                  int& SampleNo,
1566            int& DataPointNo) const
1567  {  {
1568    int i,j;    int i,j;
1569    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1257  Data::calc_mindp(int& SampleNo, Line 1599  Data::calc_mindp(int& SampleNo,
1599      }      }
1600    }    }
1601    
1602    #ifdef PASO_MPI
1603        // determine the processor on which the minimum occurs
1604        next = temp.getDataPoint(lowi,lowj)();
1605        int lowProc = 0;
1606        double *globalMins = new double[get_MPISize()+1];
1607        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1608        
1609        if( get_MPIRank()==0 ){
1610            next = globalMins[lowProc];
1611            for( i=1; i<get_MPISize(); i++ )
1612                if( next>globalMins[i] ){
1613                    lowProc = i;
1614                    next = globalMins[i];
1615                }
1616        }
1617        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1618    
1619        delete [] globalMins;
1620        ProcNo = lowProc;
1621    #else
1622        ProcNo = 0;
1623    #endif
1624    SampleNo = lowi;    SampleNo = lowi;
1625    DataPointNo = lowj;    DataPointNo = lowj;
1626  }  }
# Line 1282  Data::saveVTK(std::string fileName) cons Line 1646  Data::saveVTK(std::string fileName) cons
1646  Data&  Data&
1647  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1648  {  {
1649  #if defined DOPROF    if (isProtected()) {
1650    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1651  #endif    }
1652    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1653    return (*this);    return (*this);
1654  }  }
# Line 1292  Data::operator+=(const Data& right) Line 1656  Data::operator+=(const Data& right)
1656  Data&  Data&
1657  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1658  {  {
1659  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1660    profData->binary++;    binaryOp(tmp,plus<double>());
 #endif  
   binaryOp(right,plus<double>());  
1661    return (*this);    return (*this);
1662  }  }
1663    
1664  Data&  Data&
1665  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1666  {  {
1667  #if defined DOPROF    if (isProtected()) {
1668    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1669  #endif    }
1670    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1671    return (*this);    return (*this);
1672  }  }
# Line 1312  Data::operator-=(const Data& right) Line 1674  Data::operator-=(const Data& right)
1674  Data&  Data&
1675  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1676  {  {
1677  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1678    profData->binary++;    binaryOp(tmp,minus<double>());
 #endif  
   binaryOp(right,minus<double>());  
1679    return (*this);    return (*this);
1680  }  }
1681    
1682  Data&  Data&
1683  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1684  {  {
1685  #if defined DOPROF    if (isProtected()) {
1686    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1687  #endif    }
1688    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1689    return (*this);    return (*this);
1690  }  }
# Line 1332  Data::operator*=(const Data& right) Line 1692  Data::operator*=(const Data& right)
1692  Data&  Data&
1693  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1694  {  {
1695  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1696    profData->binary++;    binaryOp(tmp,multiplies<double>());
 #endif  
   binaryOp(right,multiplies<double>());  
1697    return (*this);    return (*this);
1698  }  }
1699    
1700  Data&  Data&
1701  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1702  {  {
1703  #if defined DOPROF    if (isProtected()) {
1704    profData->binary++;          throw DataException("Error - attempt to update protected Data object.");
1705  #endif    }
1706    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1707    return (*this);    return (*this);
1708  }  }
# Line 1352  Data::operator/=(const Data& right) Line 1710  Data::operator/=(const Data& right)
1710  Data&  Data&
1711  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1712  {  {
1713  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1714    profData->binary++;    binaryOp(tmp,divides<double>());
 #endif  
   binaryOp(right,divides<double>());  
1715    return (*this);    return (*this);
1716  }  }
1717    
1718  Data  Data
1719    Data::rpowO(const boost::python::object& left) const
1720    {
1721      Data left_d(left,*this);
1722      return left_d.powD(*this);
1723    }
1724    
1725    Data
1726  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1727  {  {
1728  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1729    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1730  }  }
1731    
1732  Data  Data
1733  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1734  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1735    Data result;    Data result;
1736    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1737    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1738         result.binaryOp(*this,escript::rpow);
1739      } else {
1740         result.copy(*this);
1741         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1742      }
1743    return result;    return result;
1744  }  }
1745    
1746    
1747  //  //
1748  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1749  Data  Data
# Line 1391  escript::operator+(const Data& left, con Line 1752  escript::operator+(const Data& left, con
1752    Data result;    Data result;
1753    //    //
1754    // perform a deep copy    // perform a deep copy
1755    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1756    result+=right;       result.copy(right);
1757         result+=left;
1758      } else {
1759         result.copy(left);
1760         result+=right;
1761      }
1762    return result;    return result;
1763  }  }
1764    
# Line 1404  escript::operator-(const Data& left, con Line 1770  escript::operator-(const Data& left, con
1770    Data result;    Data result;
1771    //    //
1772    // perform a deep copy    // perform a deep copy
1773    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1774    result-=right;       result=right.neg();
1775         result+=left;
1776      } else {
1777         result.copy(left);
1778         result-=right;
1779      }
1780    return result;    return result;
1781  }  }
1782    
# Line 1417  escript::operator*(const Data& left, con Line 1788  escript::operator*(const Data& left, con
1788    Data result;    Data result;
1789    //    //
1790    // perform a deep copy    // perform a deep copy
1791    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1792    result*=right;       result.copy(right);
1793         result*=left;
1794      } else {
1795         result.copy(left);
1796         result*=right;
1797      }
1798    return result;    return result;
1799  }  }
1800    
# Line 1430  escript::operator/(const Data& left, con Line 1806  escript::operator/(const Data& left, con
1806    Data result;    Data result;
1807    //    //
1808    // perform a deep copy    // perform a deep copy
1809    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1810    result/=right;       result=right.oneOver();
1811         result*=left;
1812      } else {
1813         result.copy(left);
1814         result/=right;
1815      }
1816    return result;    return result;
1817  }  }
1818    
# Line 1440  escript::operator/(const Data& left, con Line 1821  escript::operator/(const Data& left, con
1821  Data  Data
1822  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1823  {  {
1824    //    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;  
1825  }  }
1826    
1827  //  //
# Line 1456  escript::operator+(const Data& left, con Line 1829  escript::operator+(const Data& left, con
1829  Data  Data
1830  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1831  {  {
1832    //    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;  
1833  }  }
1834    
1835  //  //
# Line 1472  escript::operator-(const Data& left, con Line 1837  escript::operator-(const Data& left, con
1837  Data  Data
1838  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1839  {  {
1840    //    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;  
1841  }  }
1842    
1843  //  //
# Line 1488  escript::operator*(const Data& left, con Line 1845  escript::operator*(const Data& left, con
1845  Data  Data
1846  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1847  {  {
1848    //    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;  
1849  }  }
1850    
1851  //  //
# Line 1504  escript::operator/(const Data& left, con Line 1853  escript::operator/(const Data& left, con
1853  Data  Data
1854  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1855  {  {
1856    //    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;  
1857  }  }
1858    
1859  //  //
# Line 1517  escript::operator+(const boost::python:: Line 1861  escript::operator+(const boost::python::
1861  Data  Data
1862  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1863  {  {
1864    //    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;  
1865  }  }
1866    
1867  //  //
# Line 1530  escript::operator-(const boost::python:: Line 1869  escript::operator-(const boost::python::
1869  Data  Data
1870  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1871  {  {
1872    //    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;  
1873  }  }
1874    
1875  //  //
# Line 1543  escript::operator*(const boost::python:: Line 1877  escript::operator*(const boost::python::
1877  Data  Data
1878  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1879  {  {
1880    //    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;  
1881  }  }
1882    
1883  //  //
# Line 1596  escript::operator/(const boost::python:: Line 1925  escript::operator/(const boost::python::
1925  //  return ret;  //  return ret;
1926  //}  //}
1927    
1928    /* TODO */
1929    /* global reduction */
1930  Data  Data
1931  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1932  {  {
# Line 1610  Data::getItem(const boost::python::objec Line 1941  Data::getItem(const boost::python::objec
1941    return getSlice(slice_region);    return getSlice(slice_region);
1942  }  }
1943    
1944    /* TODO */
1945    /* global reduction */
1946  Data  Data
1947  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1948  {  {
# Line 1619  Data::getSlice(const DataArrayView::Regi Line 1952  Data::getSlice(const DataArrayView::Regi
1952    return Data(*this,region);    return Data(*this,region);
1953  }  }
1954    
1955    /* TODO */
1956    /* global reduction */
1957  void  void
1958  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1959                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1627  Data::setItemO(const boost::python::obje Line 1962  Data::setItemO(const boost::python::obje
1962    setItemD(key,tempData);    setItemD(key,tempData);
1963  }  }
1964    
1965    /* TODO */
1966    /* global reduction */
1967  void  void
1968  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1969                 const Data& value)                 const Data& value)
# Line 1644  Data::setItemD(const boost::python::obje Line 1981  Data::setItemD(const boost::python::obje
1981    }    }
1982  }  }
1983    
1984    /* TODO */
1985    /* global reduction */
1986  void  void
1987  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1988                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1989  {  {
1990      if (isProtected()) {
1991            throw DataException("Error - attempt to update protected Data object.");
1992      }
1993  #if defined DOPROF  #if defined DOPROF
1994    profData->slicing++;    profData->slicing++;
1995  #endif  #endif
# Line 1685  Data::typeMatchRight(const Data& right) Line 2027  Data::typeMatchRight(const Data& right)
2027    }    }
2028  }  }
2029    
2030    /* TODO */
2031    /* global reduction */
2032  void  void
2033  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2034                       const boost::python::object& value)                       const boost::python::object& value)
2035  {  {
2036      if (isProtected()) {
2037            throw DataException("Error - attempt to update protected Data object.");
2038      }
2039    //    //
2040    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2041    tag();    tag();
# Line 1706  Data::setTaggedValue(int tagKey, Line 2053  Data::setTaggedValue(int tagKey,
2053    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2054  }  }
2055    
2056    /* TODO */
2057    /* global reduction */
2058  void  void
2059  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2060                              const DataArrayView& value)                              const DataArrayView& value)
2061  {  {
2062      if (isProtected()) {
2063            throw DataException("Error - attempt to update protected Data object.");
2064      }
2065    //    //
2066    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2067    tag();    tag();
# Line 1723  Data::setTaggedValueFromCPP(int tagKey, Line 2075  Data::setTaggedValueFromCPP(int tagKey,
2075    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
2076  }  }
2077    
2078    /* TODO */
2079    /* global reduction */
2080  int  int
2081  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2082  {  {
2083    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
2084  }  }
2085    
2086    /* TODO */
2087    /* global reduction */
2088  void  void
2089  Data::setRefValue(int ref,  Data::setRefValue(int ref,
2090                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
2091  {  {
2092      if (isProtected()) {
2093            throw DataException("Error - attempt to update protected Data object.");
2094      }
2095    //    //
2096    // Construct DataArray from boost::python::object input value    // Construct DataArray from boost::python::object input value
2097    DataArray valueDataArray(value);    DataArray valueDataArray(value);
# Line 1742  Data::setRefValue(int ref, Line 2101  Data::setRefValue(int ref,
2101    m_data->setRefValue(ref,valueDataArray);    m_data->setRefValue(ref,valueDataArray);
2102  }  }
2103    
2104    /* TODO */
2105    /* global reduction */
2106  void  void
2107  Data::getRefValue(int ref,  Data::getRefValue(int ref,
2108                    boost::python::numeric::array& value)                    boost::python::numeric::array& value)
# Line 1845  Data::archiveData(const std::string file Line 2206  Data::archiveData(const std::string file
2206    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2207    int dataLength = getLength();    int dataLength = getLength();
2208    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2209    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2210    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2211      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2212    }    }
2213    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2214    if (isTagged()) {    if (isTagged()) {
2215      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2216        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 2003  Data::extractData(const std::string file Line 2364  Data::extractData(const std::string file
2364        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2365      }      }
2366    }    }
2367    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2368    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2369      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2370    }    }
2371    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2372    if (dataType==2) {    if (dataType==2) {
2373      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2374        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 2146  ostream& escript::operator<<(ostream& o, Line 2507  ostream& escript::operator<<(ostream& o,
2507    o << data.toString();    o << data.toString();
2508    return o;    return o;
2509  }  }
2510    
2511    Data
2512    escript::C_GeneralTensorProduct(Data& arg_0,
2513                         Data& arg_1,
2514                         int axis_offset,
2515                         int transpose)
2516    {
2517      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2518      // SM is the product of the last axis_offset entries in arg_0.getShape().
2519    
2520      #if defined DOPROF
2521        // profData->binary++;
2522      #endif
2523    
2524      // Interpolate if necessary and find an appropriate function space
2525      Data arg_0_Z, arg_1_Z;
2526      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2527        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2528          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2529          arg_1_Z = Data(arg_1);
2530        }
2531        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2532          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2533          arg_0_Z =Data(arg_0);
2534        }
2535        else {
2536          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2537        }
2538      } else {
2539          arg_0_Z = Data(arg_0);
2540          arg_1_Z = Data(arg_1);
2541      }
2542      // Get rank and shape of inputs
2543      int rank0 = arg_0_Z.getDataPointRank();
2544      int rank1 = arg_1_Z.getDataPointRank();
2545      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2546      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2547    
2548      // Prepare for the loops of the product and verify compatibility of shapes
2549      int start0=0, start1=0;
2550      if (transpose == 0)       {}
2551      else if (transpose == 1)  { start0 = axis_offset; }
2552      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2553      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2554    
2555      // Adjust the shapes for transpose
2556      DataArrayView::ShapeType tmpShape0;
2557      DataArrayView::ShapeType tmpShape1;
2558      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2559      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2560    
2561    #if 0
2562      // For debugging: show shape after transpose
2563      char tmp[100];
2564      std::string shapeStr;
2565      shapeStr = "(";
2566      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2567      shapeStr += ")";
2568      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2569      shapeStr = "(";
2570      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2571      shapeStr += ")";
2572      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2573    #endif
2574    
2575      // Prepare for the loops of the product
2576      int SL=1, SM=1, SR=1;
2577      for (int i=0; i<rank0-axis_offset; i++)   {
2578        SL *= tmpShape0[i];
2579      }
2580      for (int i=rank0-axis_offset; i<rank0; i++)   {
2581        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2582          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2583        }
2584        SM *= tmpShape0[i];
2585      }
2586      for (int i=axis_offset; i<rank1; i++)     {
2587        SR *= tmpShape1[i];
2588      }
2589    
2590      // Define the shape of the output
2591      DataArrayView::ShapeType shape2;
2592      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2593      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2594    
2595      // Declare output Data object
2596      Data res;
2597    
2598      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2599        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2600        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2601        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2602        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2603        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2604      }
2605      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2606    
2607        // Prepare the DataConstant input
2608        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2609        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2610    
2611        // Borrow DataTagged input from Data object
2612        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2613        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2614    
2615        // Prepare a DataTagged output 2
2616        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2617        res.tag();
2618        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2619        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2620    
2621        // Prepare offset into DataConstant
2622        int offset_0 = tmp_0->getPointOffset(0,0);
2623        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2624        // Get the views
2625        DataArrayView view_1 = tmp_1->getDefaultValue();
2626        DataArrayView view_2 = tmp_2->getDefaultValue();
2627        // Get the pointers to the actual data
2628        double *ptr_1 = &((view_1.getData())[0]);
2629        double *ptr_2 = &((view_2.getData())[0]);
2630        // Compute an MVP for the default
2631        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2632        // Compute an MVP for each tag
2633        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2634        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2635        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2636          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2637          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2638          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2639          double *ptr_1 = &view_1.getData(0);
2640          double *ptr_2 = &view_2.getData(0);
2641          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2642        }
2643    
2644      }
2645      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2646    
2647        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2648        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2649        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2650        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2651        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2652        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2653        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2654        int sampleNo_1,dataPointNo_1;
2655        int numSamples_1 = arg_1_Z.getNumSamples();
2656        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2657        int offset_0 = tmp_0->getPointOffset(0,0);
2658        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2659        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2660          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2661            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2662            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2663            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2664            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2665            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2666            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2667          }
2668        }
2669    
2670      }
2671      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2672    
2673        // Borrow DataTagged input from Data object
2674        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2675        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2676    
2677        // Prepare the DataConstant input
2678        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2679        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2680    
2681        // Prepare a DataTagged output 2
2682        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2683        res.tag();
2684        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2685        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2686    
2687        // Prepare offset into DataConstant
2688        int offset_1 = tmp_1->getPointOffset(0,0);
2689        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2690        // Get the views
2691        DataArrayView view_0 = tmp_0->getDefaultValue();
2692        DataArrayView view_2 = tmp_2->getDefaultValue();
2693        // Get the pointers to the actual data
2694        double *ptr_0 = &((view_0.getData())[0]);
2695        double *ptr_2 = &((view_2.getData())[0]);
2696        // Compute an MVP for the default
2697        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2698        // Compute an MVP for each tag
2699        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2700        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2701        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2702          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2703          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2704          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2705          double *ptr_0 = &view_0.getData(0);
2706          double *ptr_2 = &view_2.getData(0);
2707          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2708        }
2709    
2710      }
2711      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2712    
2713        // Borrow DataTagged input from Data object
2714        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2715        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2716    
2717        // Borrow DataTagged input from Data object
2718        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2719        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2720    
2721        // Prepare a DataTagged output 2
2722        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2723        res.tag();  // DataTagged output
2724        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2725        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2726    
2727        // Get the views
2728        DataArrayView view_0 = tmp_0->getDefaultValue();
2729        DataArrayView view_1 = tmp_1->getDefaultValue();
2730        DataArrayView view_2 = tmp_2->getDefaultValue();
2731        // Get the pointers to the actual data
2732        double *ptr_0 = &((view_0.getData())[0]);
2733        double *ptr_1 = &((view_1.getData())[0]);
2734        double *ptr_2 = &((view_2.getData())[0]);
2735        // Compute an MVP for the default
2736        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2737        // Merge the tags
2738        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2739        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2740        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2741        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2742          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2743        }
2744        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2745          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2746        }
2747        // Compute an MVP for each tag
2748        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2749        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2750          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2751          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2752          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2753          double *ptr_0 = &view_0.getData(0);
2754          double *ptr_1 = &view_1.getData(0);
2755          double *ptr_2 = &view_2.getData(0);
2756          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2757        }
2758    
2759      }
2760      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2761    
2762        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2763        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2764        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2765        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2766        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2767        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2768        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2769        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2770        int sampleNo_0,dataPointNo_0;
2771        int numSamples_0 = arg_0_Z.getNumSamples();
2772        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2773        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2774        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2775          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2776          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2777          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2778            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2779            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2780            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2781            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2782            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2783          }
2784        }
2785    
2786      }
2787      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2788    
2789        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2790        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2791        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2792        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2793        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2794        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2795        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2796        int sampleNo_0,dataPointNo_0;
2797        int numSamples_0 = arg_0_Z.getNumSamples();
2798        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2799        int offset_1 = tmp_1->getPointOffset(0,0);
2800        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2801        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2802          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2803            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2804            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2805            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2806            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2807            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2808            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2809          }
2810        }
2811    
2812    
2813      }
2814      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2815    
2816        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2817        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2818        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2819        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2820        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2821        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2822        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2823        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2824        int sampleNo_0,dataPointNo_0;
2825        int numSamples_0 = arg_0_Z.getNumSamples();
2826        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2827        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2828        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2829          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2830          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2831          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2832            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2833            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2834            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2835            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2836            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2837          }
2838        }
2839    
2840      }
2841      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2842    
2843        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2844        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2845        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2846        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2847        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2848        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2849        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2850        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2851        int sampleNo_0,dataPointNo_0;
2852        int numSamples_0 = arg_0_Z.getNumSamples();
2853        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2854        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2855        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2856          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2857            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2858            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2859            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2860            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2861            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2862            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2863            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2864          }
2865        }
2866    
2867      }
2868      else {
2869        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2870      }
2871    
2872      return res;
2873    }
2874    
2875    DataAbstract*
2876    Data::borrowData() const
2877    {
2878      return m_data.get();
2879    }
2880    
2881    /* Member functions specific to the MPI implementation */
2882    
2883    void
2884    Data::print()
2885    {
2886      int i,j;
2887      
2888      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2889      for( i=0; i<getNumSamples(); i++ )
2890      {
2891        printf( "[%6d]", i );
2892        for( j=0; j<getNumDataPointsPerSample(); j++ )
2893          printf( "\t%10.7g", (getSampleData(i))[j] );
2894        printf( "\n" );
2895      }
2896    }
2897    
2898    int
2899    Data::get_MPISize() const
2900    {
2901        int error, size;
2902    #ifdef PASO_MPI
2903        error = MPI_Comm_size( get_MPIComm(), &size );
2904    #else
2905        size = 1;
2906    #endif
2907        return size;
2908    }
2909    
2910    int
2911    Data::get_MPIRank() const
2912    {
2913        int error, rank;
2914    #ifdef PASO_MPI
2915        error = MPI_Comm_rank( get_MPIComm(), &rank );
2916    #else
2917        rank = 0;
2918    #endif
2919        return rank;
2920    }
2921    
2922    MPI_Comm
2923    Data::get_MPIComm() const
2924    {
2925    #ifdef PASO_MPI
2926        return MPI_COMM_WORLD;
2927    #else
2928        return -1;
2929    #endif
2930    }
2931    

Legend:
Removed from v.153  
changed lines
  Added in v.854

  ViewVC Help
Powered by ViewVC 1.1.26