/[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 121 by jgs, Fri May 6 04:26:16 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/str.hpp>  #include <boost/python/dict.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/DataAlgorithm.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;
39  using namespace escript;  using namespace escript;
40    
41    #if defined DOPROF
42    //
43    // global table of profiling data for all Data objects
44    DataProf dataProfTable;
45    #endif
46    
47  Data::Data()  Data::Data()
48  {  {
49    //    //
# Line 52  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
56      // create entry in global profiling table for this object
57      profData = dataProfTable.newData();
58    #endif
59  }  }
60    
61  Data::Data(double value,  Data::Data(double value,
# Line 65  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
74      // create entry in global profiling table for this object
75      profData = dataProfTable.newData();
76    #endif
77  }  }
78    
79  Data::Data(double value,  Data::Data(double value,
# Line 75  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
89      // create entry in global profiling table for this object
90      profData = dataProfTable.newData();
91    #endif
92  }  }
93    
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
99      // create entry in global profiling table for this object
100      profData = dataProfTable.newData();
101    #endif
102  }  }
103    
104  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  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
114      // create entry in global profiling table for this object
115      profData = dataProfTable.newData();
116    #endif
117  }  }
118    
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 for Lutz, 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
135      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
136      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
137      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
138      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
139        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  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      m_protected=false;
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 122  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    }    }
161    #if defined DOPROF
162      // create entry in global profiling table for this object
163      profData = dataProfTable.newData();
164    #endif
165  }  }
166    
167  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  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
174      // create entry in global profiling table for this object
175      profData = dataProfTable.newData();
176    #endif
177  }  }
178    
179  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  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
186      // create entry in global profiling table for this object
187      profData = dataProfTable.newData();
188    #endif
189  }  }
190    
191  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  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
199      // create entry in global profiling table for this object
200      profData = dataProfTable.newData();
201    #endif
202  }  }
203    
204  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  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
224      // create entry in global profiling table for this object
225      profData = dataProfTable.newData();
226    #endif
227    }
228    
229    Data::~Data()
230    {
231    
232  }  }
233    
234  escriptDataC  escriptDataC
# Line 288  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 303  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 344  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
435  Data::wherePositive() const  Data::wherePositive() const
436  {  {
437    #if defined DOPROF
438      profData->where++;
439    #endif
440    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
441  }  }
442    
443  Data  Data
444  Data::whereNegative() const  Data::whereNegative() const
445  {  {
446    #if defined DOPROF
447      profData->where++;
448    #endif
449    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
450  }  }
451    
452  Data  Data
453  Data::whereNonNegative() const  Data::whereNonNegative() const
454  {  {
455    #if defined DOPROF
456      profData->where++;
457    #endif
458    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
459  }  }
460    
461  Data  Data
462  Data::whereNonPositive() const  Data::whereNonPositive() const
463  {  {
464    #if defined DOPROF
465      profData->where++;
466    #endif
467    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
468  }  }
469    
470  Data  Data
471  Data::whereZero() const  Data::whereZero(double tol) const
472  {  {
473    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));  #if defined DOPROF
474      profData->where++;
475    #endif
476      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    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));  #if defined DOPROF
484      profData->where++;
485    #endif
486      Data dataAbs=abs();
487      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
488  }  }
489    
490  Data  Data
491  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
492  {  {
493    #if defined DOPROF
494      profData->interpolate++;
495    #endif
496    return Data(*this,functionspace);    return Data(*this,functionspace);
497  }  }
498    
# Line 410  Data::probeInterpolation(const FunctionS Line 514  Data::probeInterpolation(const FunctionS
514  Data  Data
515  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
516  {  {
517    #if defined DOPROF
518      profData->grad++;
519    #endif
520    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
521      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
522    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 550  Data::getDataPointShape() const
550    return getPointDataView().getShape();    return getPointDataView().getShape();
551  }  }
552    
553    void
554    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
561      if (num_array.getrank()<getDataPointRank())
562          throw DataException("Rank of numarray does not match Data object rank");
563    
564      //
565      // check shape of num_array
566      for (int i=0; i<getDataPointRank(); i++) {
567        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
568           throw DataException("Shape of numarray does not match Data object rank");
569      }
570    
571      //
572      // make sure data is expanded:
573      if (!isExpanded()) {
574        expand();
575      }
576    
577      //
578      // and copy over
579      m_data->copyAll(num_array);
580    }
581    
582  const  const
583  boost::python::numeric::array  boost::python::numeric::array
584  Data::convertToNumArray()  Data::convertToNumArray()
# Line 649  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 697  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 735  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 748  Data::integrate() const Line 967  Data::integrate() const
967    int rank = getDataPointRank();    int rank = getDataPointRank();
968    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
969    
970    #if defined DOPROF
971      profData->integrate++;
972    #endif
973    
974    //    //
975    // calculate the integral values    // calculate the integral values
976    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 770  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 784  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 796  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 811  Data::integrate() const Line 1034  Data::integrate() const
1034  Data  Data
1035  Data::sin() const  Data::sin() const
1036  {  {
1037    #if defined DOPROF
1038      profData->unary++;
1039    #endif
1040    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1041  }  }
1042    
1043  Data  Data
1044  Data::cos() const  Data::cos() const
1045  {  {
1046    #if defined DOPROF
1047      profData->unary++;
1048    #endif
1049    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1050  }  }
1051    
1052  Data  Data
1053  Data::tan() const  Data::tan() const
1054  {  {
1055    #if defined DOPROF
1056      profData->unary++;
1057    #endif
1058    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1059  }  }
1060    
1061  Data  Data
1062  Data::log() const  Data::asin() const
1063    {
1064    #if defined DOPROF
1065      profData->unary++;
1066    #endif
1067      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1068    }
1069    
1070    Data
1071    Data::acos() const
1072    {
1073    #if defined DOPROF
1074      profData->unary++;
1075    #endif
1076      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1077    }
1078    
1079    Data
1080    Data::atan() const
1081    {
1082    #if defined DOPROF
1083      profData->unary++;
1084    #endif
1085      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1086    }
1087    
1088    Data
1089    Data::sinh() const
1090    {
1091    #if defined DOPROF
1092      profData->unary++;
1093    #endif
1094      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1095    }
1096    
1097    Data
1098    Data::cosh() const
1099    {
1100    #if defined DOPROF
1101      profData->unary++;
1102    #endif
1103      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1104    }
1105    
1106    Data
1107    Data::tanh() const
1108    {
1109    #if defined DOPROF
1110      profData->unary++;
1111    #endif
1112      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1113    }
1114    
1115    Data
1116    Data::asinh() const
1117    {
1118    #if defined DOPROF
1119      profData->unary++;
1120    #endif
1121      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1122    }
1123    
1124    Data
1125    Data::acosh() const
1126    {
1127    #if defined DOPROF
1128      profData->unary++;
1129    #endif
1130      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1131    }
1132    
1133    Data
1134    Data::atanh() const
1135  {  {
1136    #if defined DOPROF
1137      profData->unary++;
1138    #endif
1139      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1140    }
1141    
1142    Data
1143    Data::log10() const
1144    {
1145    #if defined DOPROF
1146      profData->unary++;
1147    #endif
1148    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1149  }  }
1150    
1151  Data  Data
1152  Data::ln() const  Data::log() const
1153  {  {
1154    #if defined DOPROF
1155      profData->unary++;
1156    #endif
1157    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1158  }  }
1159    
1160  Data  Data
1161  Data::sign() const  Data::sign() const
1162  {  {
1163    #if defined DOPROF
1164      profData->unary++;
1165    #endif
1166    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1167  }  }
1168    
1169  Data  Data
1170  Data::abs() const  Data::abs() const
1171  {  {
1172    #if defined DOPROF
1173      profData->unary++;
1174    #endif
1175    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1176  }  }
1177    
1178  Data  Data
1179  Data::neg() const  Data::neg() const
1180  {  {
1181    #if defined DOPROF
1182      profData->unary++;
1183    #endif
1184    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1185  }  }
1186    
1187  Data  Data
1188  Data::pos() const  Data::pos() const
1189  {  {
1190    return (*this);  #if defined DOPROF
1191      profData->unary++;
1192    #endif
1193      Data result;
1194      // perform a deep copy
1195      result.copy(*this);
1196      return result;
1197  }  }
1198    
1199  Data  Data
1200  Data::exp() const  Data::exp() const
1201  {  {
1202    #if defined DOPROF
1203      profData->unary++;
1204    #endif
1205    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1206  }  }
1207    
1208  Data  Data
1209  Data::sqrt() const  Data::sqrt() const
1210  {  {
1211    #if defined DOPROF
1212      profData->unary++;
1213    #endif
1214    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1215  }  }
1216    
1217  double  double
1218  Data::Lsup() const  Data::Lsup() const
1219  {  {
1220      double localValue, globalValue;
1221    #if defined DOPROF
1222      profData->reduction1++;
1223    #endif
1224    //    //
1225    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1226    return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
1227      AbsMax abs_max_func;
1228      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
1242      profData->reduction1++;
1243    #endif
1244    //    //
1245    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1246    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1247      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
1262      profData->reduction1++;
1263    #endif
1264    //    //
1265    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1266    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1267      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
1281      profData->reduction1++;
1282    #endif
1283    //    //
1284    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1285    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1286      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  {  {
1300    #if defined DOPROF
1301      profData->reduction2++;
1302    #endif
1303    //    //
1304    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1305    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1306      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1307  }  }
1308    
1309  Data  Data
1310  Data::minval() const  Data::minval() const
1311  {  {
1312    #if defined DOPROF
1313      profData->reduction2++;
1314    #endif
1315    //    //
1316    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1317    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1318      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1319    }
1320    
1321    Data
1322    Data::swapaxes(const int axis0, const int axis1) const
1323    {
1324         int axis0_tmp,axis1_tmp;
1325         #if defined DOPROF
1326         profData->unary++;
1327         #endif
1328         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
1484    Data::transpose(int axis_offset) const
1485    {
1486         #if defined DOPROF
1487         profData->unary++;
1488         #endif
1489         DataArrayView::ShapeType s=getDataPointShape();
1490         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
1508    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
1531            profData->unary++;
1532         #endif
1533         DataArrayView::ShapeType s=getDataPointShape();
1534         if (getDataPointRank()!=2)
1535            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
1550  Data::mindp() const  Data::mindp() const
1551  {  {
1552      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1553      // abort (for unknown reasons) if there are openmp directives with it in the
1554      // surrounding function
1555    
1556      int SampleNo;
1557      int DataPointNo;
1558        int ProcNo;
1559      calc_mindp(ProcNo,SampleNo,DataPointNo);
1560      return make_tuple(ProcNo,SampleNo,DataPointNo);
1561    }
1562    
1563    void
1564    Data::calc_mindp(   int& ProcNo,
1565                    int& SampleNo,
1566            int& DataPointNo) const
1567    {
1568      int i,j;
1569      int lowi=0,lowj=0;
1570      double min=numeric_limits<double>::max();
1571    
1572    Data temp=minval();    Data temp=minval();
1573    
1574    int numSamples=temp.getNumSamples();    int numSamples=temp.getNumSamples();
1575    int numDPPSample=temp.getNumDataPointsPerSample();    int numDPPSample=temp.getNumDataPointsPerSample();
1576    
1577    int i,j,lowi=0,lowj=0;    double next,local_min;
1578    double min=numeric_limits<double>::max();    int local_lowi,local_lowj;
1579    
1580    for (i=0; i<numSamples; i++) {    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1581      for (j=0; j<numDPPSample; j++) {    {
1582        double next=temp.getDataPoint(i,j)();      local_min=min;
1583        if (next<min) {      #pragma omp for private(i,j) schedule(static)
1584          min=next;      for (i=0; i<numSamples; i++) {
1585          lowi=i;        for (j=0; j<numDPPSample; j++) {
1586          lowj=j;          next=temp.getDataPoint(i,j)();
1587            if (next<local_min) {
1588              local_min=next;
1589              local_lowi=i;
1590              local_lowj=j;
1591            }
1592        }        }
1593      }      }
1594        #pragma omp critical
1595        if (local_min<min) {
1596          min=local_min;
1597          lowi=local_lowi;
1598          lowj=local_lowj;
1599        }
1600    }    }
1601    
1602    return make_tuple(lowi,lowj);  #ifdef PASO_MPI
1603  }      // determine the processor on which the minimum occurs
1604        next = temp.getDataPoint(lowi,lowj)();
1605  Data      int lowProc = 0;
1606  Data::length() const      double *globalMins = new double[get_MPISize()+1];
1607  {      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1608    return dp_algorithm(DataAlgorithmAdapter<Length>(0));      
1609  }      if( get_MPIRank()==0 ){
1610            next = globalMins[lowProc];
1611  Data          for( i=1; i<get_MPISize(); i++ )
1612  Data::trace() const              if( next>globalMins[i] ){
1613  {                  lowProc = i;
1614    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));                  next = globalMins[i];
1615  }              }
1616        }
1617  Data      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1618  Data::transpose(int axis) const  
1619  {      delete [] globalMins;
1620    // not implemented      ProcNo = lowProc;
1621    throw DataException("Error - Data::transpose not implemented yet.");  #else
1622    return Data();      ProcNo = 0;
1623    #endif
1624      SampleNo = lowi;
1625      DataPointNo = lowj;
1626  }  }
1627    
1628  void  void
1629  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1630  {  {
1631    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1632      args["data"]=boost::python::object(this);
1633      getDomain().saveDX(fileName,args);
1634    return;    return;
1635  }  }
1636    
1637  void  void
1638  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1639  {  {
1640    getDomain().saveVTK(fileName,*this);    boost::python::dict args;
1641      args["data"]=boost::python::object(this);
1642      getDomain().saveVTK(fileName,args);
1643    return;    return;
1644  }  }
1645    
1646  Data&  Data&
1647  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1648  {  {
1649      if (isProtected()) {
1650            throw DataException("Error - attempt to update protected Data object.");
1651      }
1652    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1653    return (*this);    return (*this);
1654  }  }
# Line 991  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    binaryOp(right,plus<double>());    Data tmp(right,getFunctionSpace(),false);
1660      binaryOp(tmp,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 (isProtected()) {
1668            throw DataException("Error - attempt to update protected Data object.");
1669      }
1670    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1671    return (*this);    return (*this);
1672  }  }
# Line 1005  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    binaryOp(right,minus<double>());    Data tmp(right,getFunctionSpace(),false);
1678      binaryOp(tmp,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 (isProtected()) {
1686            throw DataException("Error - attempt to update protected Data object.");
1687      }
1688    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1689    return (*this);    return (*this);
1690  }  }
# Line 1019  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    binaryOp(right,multiplies<double>());    Data tmp(right,getFunctionSpace(),false);
1696      binaryOp(tmp,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 (isProtected()) {
1704            throw DataException("Error - attempt to update protected Data object.");
1705      }
1706    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1707    return (*this);    return (*this);
1708  }  }
# Line 1033  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    binaryOp(right,divides<double>());    Data tmp(right,getFunctionSpace(),false);
1714      binaryOp(tmp,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    Data result;    Data tmp(right,getFunctionSpace(),false);
1729    result.copy(*this);    return powD(tmp);
   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  {  {
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 namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1749  Data  Data
1750  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1751  {  {
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    
1765  //  //
1766  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1767  Data  Data
1768  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1769  {  {
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    
1783  //  //
1784  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1785  Data  Data
1786  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1787  {  {
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    
1801  //  //
1802  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1803  Data  Data
1804  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1805  {  {
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    
1819  //  //
1820  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1828  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1836  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1844  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1852  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1860  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1868  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
1876  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
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  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1884  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1885  //{  //{
1886  //  /*  //  /*
# Line 1269  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 1283  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  {  {
1949    #if defined DOPROF
1950      profData->slicing++;
1951    #endif
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 1297  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)
1970  {  {
1971    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1972    
1973    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1974    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1975      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1313  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
1994      profData->slicing++;
1995    #endif
1996    Data tempValue(value);    Data tempValue(value);
1997    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1998    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1351  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 1372  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 1389  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
2081    Data::getTagNumber(int dpno)
2082    {
2083      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 1402  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 1428  Data::getRefValue(int ref, Line 2129  Data::getRefValue(int ref,
2129    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
2130    
2131    if (rank==0) {    if (rank==0) {
2132      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
2133          value = temp_numArray;
2134    }    }
2135    if (rank==1) {    if (rank==1) {
2136      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1436  Data::getRefValue(int ref, Line 2138  Data::getRefValue(int ref,
2138      }      }
2139    }    }
2140    if (rank==2) {    if (rank==2) {
2141      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2142          for (int j=0; j < shape[1]; j++) {
2143            value[i][j] = valueView(i,j);
2144          }
2145        }
2146    }    }
2147    if (rank==3) {    if (rank==3) {
2148      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2149          for (int j=0; j < shape[1]; j++) {
2150            for (int k=0; k < shape[2]; k++) {
2151              value[i][j][k] = valueView(i,j,k);
2152            }
2153          }
2154        }
2155    }    }
2156    if (rank==4) {    if (rank==4) {
2157      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2158          for (int j=0; j < shape[1]; j++) {
2159            for (int k=0; k < shape[2]; k++) {
2160              for (int l=0; l < shape[3]; l++) {
2161                value[i][j][k][l] = valueView(i,j,k,l);
2162              }
2163            }
2164          }
2165        }
2166    }    }
2167    
2168  }  }
# Line 1472  Data::archiveData(const std::string file Line 2192  Data::archiveData(const std::string file
2192      dataType = 3;      dataType = 3;
2193      cout << "\tdataType: DataExpanded" << endl;      cout << "\tdataType: DataExpanded" << endl;
2194    }    }
2195    
2196    if (dataType == -1) {    if (dataType == -1) {
2197      throw DataException("archiveData Error: undefined dataType");      throw DataException("archiveData Error: undefined dataType");
2198    }    }
# Line 1485  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 1511  Data::archiveData(const std::string file Line 2232  Data::archiveData(const std::string file
2232    cout << ">" << endl;    cout << ">" << endl;
2233    
2234    //    //
2235    // Write common data items to archive file    // Open archive file
2236    ofstream archiveFile;    ofstream archiveFile;
2237    archiveFile.open(fileName.data(), ios::out);    archiveFile.open(fileName.data(), ios::out);
2238    
# Line 1519  Data::archiveData(const std::string file Line 2240  Data::archiveData(const std::string file
2240      throw DataException("archiveData Error: problem opening archive file");      throw DataException("archiveData Error: problem opening archive file");
2241    }    }
2242    
2243      //
2244      // Write common data items to archive file
2245    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2246    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2247    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1542  Data::archiveData(const std::string file Line 2265  Data::archiveData(const std::string file
2265      throw DataException("archiveData Error: problem writing to archive file");      throw DataException("archiveData Error: problem writing to archive file");
2266    }    }
2267    
   archiveFile.close();  
   
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
   }  
   
2268    //    //
2269    // Collect and archive underlying data values for each Data type    // Archive underlying data values for each Data type
2270      int noValues;
2271    switch (dataType) {    switch (dataType) {
2272      case 0:      case 0:
2273        // DataEmpty        // DataEmpty
2274          noValues = 0;
2275          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2276          cout << "\tnoValues: " << noValues << endl;
2277        break;        break;
2278      case 1:      case 1:
2279        // DataConstant        // DataConstant
2280          noValues = m_data->getLength();
2281          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2282          cout << "\tnoValues: " << noValues << endl;
2283          if (m_data->archiveData(archiveFile,noValues)) {
2284            throw DataException("archiveData Error: problem writing data to archive file");
2285          }
2286        break;        break;
2287      case 2:      case 2:
2288        // DataTagged        // DataTagged
2289          noValues = m_data->getLength();
2290          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2291          cout << "\tnoValues: " << noValues << endl;
2292          if (m_data->archiveData(archiveFile,noValues)) {
2293            throw DataException("archiveData Error: problem writing data to archive file");
2294          }
2295        break;        break;
2296      case 3:      case 3:
2297        // DataExpanded        // DataExpanded
2298          noValues = m_data->getLength();
2299          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2300          cout << "\tnoValues: " << noValues << endl;
2301          if (m_data->archiveData(archiveFile,noValues)) {
2302            throw DataException("archiveData Error: problem writing data to archive file");
2303          }
2304        break;        break;
2305    }    }
2306    
2307      if (!archiveFile.good()) {
2308        throw DataException("archiveData Error: problem writing data to archive file");
2309      }
2310    
2311      //
2312      // Close archive file
2313      archiveFile.close();
2314    
2315      if (!archiveFile.good()) {
2316        throw DataException("archiveData Error: problem closing archive file");
2317      }
2318    
2319  }  }
2320    
2321  void  void
# Line 1590  Data::extractData(const std::string file Line 2341  Data::extractData(const std::string file
2341    int flatShape[4];    int flatShape[4];
2342    
2343    //    //
2344    // Open the archive file and read common data items    // Open the archive file
2345    ifstream archiveFile;    ifstream archiveFile;
2346    archiveFile.open(fileName.data(), ios::in);    archiveFile.open(fileName.data(), ios::in);
2347    
# Line 1598  Data::extractData(const std::string file Line 2349  Data::extractData(const std::string file
2349      throw DataException("extractData Error: problem opening archive file");      throw DataException("extractData Error: problem opening archive file");
2350    }    }
2351    
2352      //
2353      // Read common data items from archive file
2354    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2355    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2356    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1611  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 1626  Data::extractData(const std::string file Line 2379  Data::extractData(const std::string file
2379      throw DataException("extractData Error: problem reading from archive file");      throw DataException("extractData Error: problem reading from archive file");
2380    }    }
2381    
2382    archiveFile.close();    //
2383      // Verify the values just read from the archive file
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem closing archive file");  
   }  
   
2384    switch (dataType) {    switch (dataType) {
2385      case 0:      case 0:
2386        cout << "\tdataType: DataEmpty" << endl;        cout << "\tdataType: DataEmpty" << endl;
# Line 1686  Data::extractData(const std::string file Line 2435  Data::extractData(const std::string file
2435    
2436    //    //
2437    // Load this DataVector with the appropriate values    // Load this DataVector with the appropriate values
2438      int noValues;
2439      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2440      cout << "\tnoValues: " << noValues << endl;
2441    switch (dataType) {    switch (dataType) {
2442      case 0:      case 0:
2443        // DataEmpty        // DataEmpty
2444          if (noValues != 0) {
2445            throw DataException("extractData Error: problem reading data from archive file");
2446          }
2447        break;        break;
2448      case 1:      case 1:
2449        // DataConstant        // DataConstant
2450          if (dataVec.extractData(archiveFile,noValues)) {
2451            throw DataException("extractData Error: problem reading data from archive file");
2452          }
2453        break;        break;
2454      case 2:      case 2:
2455        // DataTagged        // DataTagged
2456          if (dataVec.extractData(archiveFile,noValues)) {
2457            throw DataException("extractData Error: problem reading data from archive file");
2458          }
2459        break;        break;
2460      case 3:      case 3:
2461        // DataExpanded        // DataExpanded
2462          if (dataVec.extractData(archiveFile,noValues)) {
2463            throw DataException("extractData Error: problem reading data from archive file");
2464          }
2465        break;        break;
2466    }    }
2467    
2468      if (!archiveFile.good()) {
2469        throw DataException("extractData Error: problem reading from archive file");
2470      }
2471    
2472      //
2473      // Close archive file
2474      archiveFile.close();
2475    
2476      if (!archiveFile.good()) {
2477        throw DataException("extractData Error: problem closing archive file");
2478      }
2479    
2480    //    //
2481    // Construct an appropriate Data object    // Construct an appropriate Data object
2482    DataAbstract* tempData;    DataAbstract* tempData;
# Line 1724  Data::extractData(const std::string file Line 2500  Data::extractData(const std::string file
2500    }    }
2501    shared_ptr<DataAbstract> temp_data(tempData);    shared_ptr<DataAbstract> temp_data(tempData);
2502    m_data=temp_data;    m_data=temp_data;
   
2503  }  }
2504    
2505  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
# Line 1732  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.121  
changed lines
  Added in v.854

  ViewVC Help
Powered by ViewVC 1.1.26