/[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 782 by bcumming, Tue Jul 18 00:47:47 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    #if defined DOPROF
55      // create entry in global profiling table for this object
56      profData = dataProfTable.newData();
57    #endif
58  }  }
59    
60  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 68  Data::Data(double value,
68    }    }
69    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
70    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
71    #if defined DOPROF
72      // create entry in global profiling table for this object
73      profData = dataProfTable.newData();
74    #endif
75  }  }
76    
77  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 82  Data::Data(double value,
82    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
83    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
84    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
85    #if defined DOPROF
86      // create entry in global profiling table for this object
87      profData = dataProfTable.newData();
88    #endif
89  }  }
90    
91  Data::Data(const Data& inData)  Data::Data(const Data& inData)
92  {  {
93    m_data=inData.m_data;    m_data=inData.m_data;
94    #if defined DOPROF
95      // create entry in global profiling table for this object
96      profData = dataProfTable.newData();
97    #endif
98  }  }
99    
100  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 105  Data::Data(const Data& inData,
105    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
106    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
107    m_data=temp_data;    m_data=temp_data;
108    #if defined DOPROF
109      // create entry in global profiling table for this object
110      profData = dataProfTable.newData();
111    #endif
112  }  }
113    
114  Data::Data(const Data& inData,  Data::Data(const Data& inData,
115             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
116  {  {
117    #if defined DOPROF
118      // create entry in global profiling table for this object
119      profData = dataProfTable.newData();
120    #endif
121    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
122      m_data=inData.m_data;      m_data=inData.m_data;
123    } else {    } else {
124        #if defined DOPROF
125        profData->interpolate++;
126        #endif
127      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
128      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
129      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
130      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
131      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
132      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
133      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
134        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 125  Data::Data(const DataTagged::TagListType Line 151  Data::Data(const DataTagged::TagListType
151    if (expanded) {    if (expanded) {
152      expand();      expand();
153    }    }
154    #if defined DOPROF
155      // create entry in global profiling table for this object
156      profData = dataProfTable.newData();
157    #endif
158  }  }
159    
160  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 162  Data::Data(const numeric::array& value,
162             bool expanded)             bool expanded)
163  {  {
164    initialise(value,what,expanded);    initialise(value,what,expanded);
165    #if defined DOPROF
166      // create entry in global profiling table for this object
167      profData = dataProfTable.newData();
168    #endif
169  }  }
170    
171  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 173  Data::Data(const DataArrayView& value,
173             bool expanded)             bool expanded)
174  {  {
175    initialise(value,what,expanded);    initialise(value,what,expanded);
176    #if defined DOPROF
177      // create entry in global profiling table for this object
178      profData = dataProfTable.newData();
179    #endif
180  }  }
181    
182  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 185  Data::Data(const object& value,
185  {  {
186    numeric::array asNumArray(value);    numeric::array asNumArray(value);
187    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
188    #if defined DOPROF
189      // create entry in global profiling table for this object
190      profData = dataProfTable.newData();
191    #endif
192  }  }
193    
194  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 209  Data::Data(const object& value,
209      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
210      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
211    }    }
212    #if defined DOPROF
213      // create entry in global profiling table for this object
214      profData = dataProfTable.newData();
215    #endif
216    }
217    
218    Data::~Data()
219    {
220    
221  }  }
222    
223  escriptDataC  escriptDataC
# Line 288  Data::isTagged() const Line 339  Data::isTagged() const
339    return (temp!=0);    return (temp!=0);
340  }  }
341    
342    /* TODO */
343    /* global reduction -- the local data being empty does not imply that it is empty on other processers*/
344  bool  bool
345  Data::isEmpty() const  Data::isEmpty() const
346  {  {
# Line 353  Data::reshapeDataPoint(const DataArrayVi Line 406  Data::reshapeDataPoint(const DataArrayVi
406  Data  Data
407  Data::wherePositive() const  Data::wherePositive() const
408  {  {
409    #if defined DOPROF
410      profData->where++;
411    #endif
412    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
413  }  }
414    
415  Data  Data
416  Data::whereNegative() const  Data::whereNegative() const
417  {  {
418    #if defined DOPROF
419      profData->where++;
420    #endif
421    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
422  }  }
423    
424  Data  Data
425  Data::whereNonNegative() const  Data::whereNonNegative() const
426  {  {
427    #if defined DOPROF
428      profData->where++;
429    #endif
430    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
431  }  }
432    
433  Data  Data
434  Data::whereNonPositive() const  Data::whereNonPositive() const
435  {  {
436    #if defined DOPROF
437      profData->where++;
438    #endif
439    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
440  }  }
441    
442  Data  Data
443  Data::whereZero() const  Data::whereZero(double tol) const
444  {  {
445    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));  #if defined DOPROF
446      profData->where++;
447    #endif
448      Data dataAbs=abs();
449      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
450  }  }
451    
452  Data  Data
453  Data::whereNonZero() const  Data::whereNonZero(double tol) const
454  {  {
455    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));  #if defined DOPROF
456      profData->where++;
457    #endif
458      Data dataAbs=abs();
459      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
460  }  }
461    
462  Data  Data
463  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
464  {  {
465    #if defined DOPROF
466      profData->interpolate++;
467    #endif
468    return Data(*this,functionspace);    return Data(*this,functionspace);
469  }  }
470    
# Line 410  Data::probeInterpolation(const FunctionS Line 486  Data::probeInterpolation(const FunctionS
486  Data  Data
487  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
488  {  {
489    #if defined DOPROF
490      profData->grad++;
491    #endif
492    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
493      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
494    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 522  Data::getDataPointShape() const
522    return getPointDataView().getShape();    return getPointDataView().getShape();
523  }  }
524    
525    void
526    Data::fillFromNumArray(const boost::python::numeric::array num_array)
527    {
528      //
529      // check rank
530      if (num_array.getrank()<getDataPointRank())
531          throw DataException("Rank of numarray does not match Data object rank");
532    
533      //
534      // check shape of num_array
535      for (int i=0; i<getDataPointRank(); i++) {
536        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
537           throw DataException("Shape of numarray does not match Data object rank");
538      }
539    
540      //
541      // make sure data is expanded:
542      if (!isExpanded()) {
543        expand();
544      }
545    
546      //
547      // and copy over
548      m_data->copyAll(num_array);
549    }
550    
551  const  const
552  boost::python::numeric::array  boost::python::numeric::array
553  Data::convertToNumArray()  Data::convertToNumArray()
# Line 649  Data::convertToNumArrayFromSampleNo(int Line 754  Data::convertToNumArrayFromSampleNo(int
754    
755  const  const
756  boost::python::numeric::array  boost::python::numeric::array
757    #ifdef PASO_MPI
758    Data::convertToNumArrayFromDPNo(int procNo,
759                                    int sampleNo,
760                                                                    int dataPointNo)
761    
762    #else
763  Data::convertToNumArrayFromDPNo(int sampleNo,  Data::convertToNumArrayFromDPNo(int sampleNo,
764                                  int dataPointNo)                                  int dataPointNo)
765    #endif
766  {  {
767    //    //
768    // Check a valid sample number has been supplied    // Check a valid sample number has been supplied
# Line 748  Data::integrate() const Line 860  Data::integrate() const
860    int rank = getDataPointRank();    int rank = getDataPointRank();
861    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
862    
863    #if defined DOPROF
864      profData->integrate++;
865    #endif
866    
867    //    //
868    // calculate the integral values    // calculate the integral values
869    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 770  Data::integrate() const Line 886  Data::integrate() const
886      }      }
887    }    }
888    if (rank==2) {    if (rank==2) {
889      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
890      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
891        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
892          index = i + shape[0] * j;             index = i + shape[0] * j;
893          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
894        }           }
895      }         }
896    }    }
897    if (rank==3) {    if (rank==3) {
898      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 784  Data::integrate() const Line 900  Data::integrate() const
900        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
901          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
902            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
903            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
904          }          }
905        }        }
906      }      }
# Line 796  Data::integrate() const Line 912  Data::integrate() const
912          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
913            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
914              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
915              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
916            }            }
917          }          }
918        }        }
# Line 811  Data::integrate() const Line 927  Data::integrate() const
927  Data  Data
928  Data::sin() const  Data::sin() const
929  {  {
930    #if defined DOPROF
931      profData->unary++;
932    #endif
933    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
934  }  }
935    
936  Data  Data
937  Data::cos() const  Data::cos() const
938  {  {
939    #if defined DOPROF
940      profData->unary++;
941    #endif
942    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
943  }  }
944    
945  Data  Data
946  Data::tan() const  Data::tan() const
947  {  {
948    #if defined DOPROF
949      profData->unary++;
950    #endif
951    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
952  }  }
953    
954  Data  Data
955  Data::log() const  Data::asin() const
956    {
957    #if defined DOPROF
958      profData->unary++;
959    #endif
960      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
961    }
962    
963    Data
964    Data::acos() const
965    {
966    #if defined DOPROF
967      profData->unary++;
968    #endif
969      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
970    }
971    
972    Data
973    Data::atan() const
974    {
975    #if defined DOPROF
976      profData->unary++;
977    #endif
978      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
979    }
980    
981    Data
982    Data::sinh() const
983    {
984    #if defined DOPROF
985      profData->unary++;
986    #endif
987      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
988    }
989    
990    Data
991    Data::cosh() const
992  {  {
993    #if defined DOPROF
994      profData->unary++;
995    #endif
996      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
997    }
998    
999    Data
1000    Data::tanh() const
1001    {
1002    #if defined DOPROF
1003      profData->unary++;
1004    #endif
1005      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1006    }
1007    
1008    Data
1009    Data::asinh() const
1010    {
1011    #if defined DOPROF
1012      profData->unary++;
1013    #endif
1014      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1015    }
1016    
1017    Data
1018    Data::acosh() const
1019    {
1020    #if defined DOPROF
1021      profData->unary++;
1022    #endif
1023      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1024    }
1025    
1026    Data
1027    Data::atanh() const
1028    {
1029    #if defined DOPROF
1030      profData->unary++;
1031    #endif
1032      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1033    }
1034    
1035    Data
1036    Data::log10() const
1037    {
1038    #if defined DOPROF
1039      profData->unary++;
1040    #endif
1041    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1042  }  }
1043    
1044  Data  Data
1045  Data::ln() const  Data::log() const
1046  {  {
1047    #if defined DOPROF
1048      profData->unary++;
1049    #endif
1050    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1051  }  }
1052    
1053  Data  Data
1054  Data::sign() const  Data::sign() const
1055  {  {
1056    #if defined DOPROF
1057      profData->unary++;
1058    #endif
1059    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1060  }  }
1061    
1062  Data  Data
1063  Data::abs() const  Data::abs() const
1064  {  {
1065    #if defined DOPROF
1066      profData->unary++;
1067    #endif
1068    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1069  }  }
1070    
1071  Data  Data
1072  Data::neg() const  Data::neg() const
1073  {  {
1074    #if defined DOPROF
1075      profData->unary++;
1076    #endif
1077    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1078  }  }
1079    
1080  Data  Data
1081  Data::pos() const  Data::pos() const
1082  {  {
1083    return (*this);  #if defined DOPROF
1084      profData->unary++;
1085    #endif
1086      Data result;
1087      // perform a deep copy
1088      result.copy(*this);
1089      return result;
1090  }  }
1091    
1092  Data  Data
1093  Data::exp() const  Data::exp() const
1094  {  {
1095    #if defined DOPROF
1096      profData->unary++;
1097    #endif
1098    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1099  }  }
1100    
1101  Data  Data
1102  Data::sqrt() const  Data::sqrt() const
1103  {  {
1104    #if defined DOPROF
1105      profData->unary++;
1106    #endif
1107    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1108  }  }
1109    
1110  double  double
1111  Data::Lsup() const  Data::Lsup() const
1112  {  {
1113      double localValue, globalValue;
1114    #if defined DOPROF
1115      profData->reduction1++;
1116    #endif
1117    //    //
1118    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1119    return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
1120      AbsMax abs_max_func;
1121      localValue = algorithm(abs_max_func,0);
1122    #ifdef PASO_MPI
1123      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1124      return globalValue;
1125    #else
1126      return localValue;
1127    #endif
1128  }  }
1129    
1130  double  double
1131  Data::Linf() const  Data::Linf() const
1132  {  {
1133      double localValue, globalValue;
1134    #if defined DOPROF
1135      profData->reduction1++;
1136    #endif
1137    //    //
1138    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1139    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1140      localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1141    
1142    #ifdef PASO_MPI
1143      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1144      return globalValue;
1145    #else
1146      return localValue;
1147    #endif
1148  }  }
1149    
1150  double  double
1151  Data::sup() const  Data::sup() const
1152  {  {
1153      double localValue, globalValue;
1154    #if defined DOPROF
1155      profData->reduction1++;
1156    #endif
1157    //    //
1158    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1159    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1160      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1161    #ifdef PASO_MPI
1162      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1163      return globalValue;
1164    #else
1165      return localValue;
1166    #endif
1167  }  }
1168    
1169  double  double
1170  Data::inf() const  Data::inf() const
1171  {  {
1172      double localValue, globalValue;
1173    #if defined DOPROF
1174      profData->reduction1++;
1175    #endif
1176    //    //
1177    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1178    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1179      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1180    #ifdef PASO_MPI
1181      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1182      return globalValue;
1183    #else
1184      return localValue;
1185    #endif
1186  }  }
1187    
1188    /* TODO */
1189    /* global reduction */
1190  Data  Data
1191  Data::maxval() const  Data::maxval() const
1192  {  {
1193    #if defined DOPROF
1194      profData->reduction2++;
1195    #endif
1196    //    //
1197    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1198    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1199      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1200  }  }
1201    
1202  Data  Data
1203  Data::minval() const  Data::minval() const
1204  {  {
1205    #if defined DOPROF
1206      profData->reduction2++;
1207    #endif
1208    //    //
1209    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1210    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1211      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1212    }
1213    
1214    Data
1215    Data::trace() const
1216    {
1217    #if defined DOPROF
1218      profData->reduction2++;
1219    #endif
1220      Trace trace_func;
1221      return dp_algorithm(trace_func,0);
1222    }
1223    
1224    Data
1225    Data::symmetric() const
1226    {
1227         #if defined DOPROF
1228            profData->unary++;
1229         #endif
1230         // check input
1231         DataArrayView::ShapeType s=getDataPointShape();
1232         if (getDataPointRank()==2) {
1233            if(s[0] != s[1])
1234               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1235         }
1236         else if (getDataPointRank()==4) {
1237            if(!(s[0] == s[2] && s[1] == s[3]))
1238               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1239         }
1240         else {
1241            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1242         }
1243         Data ev(0.,getDataPointShape(),getFunctionSpace());
1244         ev.typeMatchRight(*this);
1245         m_data->symmetric(ev.m_data.get());
1246         return ev;
1247    }
1248    
1249    Data
1250    Data::nonsymmetric() const
1251    {
1252         #if defined DOPROF
1253            profData->unary++;
1254         #endif
1255         // check input
1256         DataArrayView::ShapeType s=getDataPointShape();
1257         if (getDataPointRank()==2) {
1258            if(s[0] != s[1])
1259               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1260            DataArrayView::ShapeType ev_shape;
1261            ev_shape.push_back(s[0]);
1262            ev_shape.push_back(s[1]);
1263            Data ev(0.,ev_shape,getFunctionSpace());
1264            ev.typeMatchRight(*this);
1265            m_data->nonsymmetric(ev.m_data.get());
1266            return ev;
1267         }
1268         else if (getDataPointRank()==4) {
1269            if(!(s[0] == s[2] && s[1] == s[3]))
1270               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1271            DataArrayView::ShapeType ev_shape;
1272            ev_shape.push_back(s[0]);
1273            ev_shape.push_back(s[1]);
1274            ev_shape.push_back(s[2]);
1275            ev_shape.push_back(s[3]);
1276            Data ev(0.,ev_shape,getFunctionSpace());
1277            ev.typeMatchRight(*this);
1278            m_data->nonsymmetric(ev.m_data.get());
1279            return ev;
1280         }
1281         else {
1282            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1283         }
1284    }
1285    
1286    Data
1287    Data::matrixtrace(int axis_offset) const
1288    {
1289         #if defined DOPROF
1290            profData->unary++;
1291         #endif
1292         DataArrayView::ShapeType s=getDataPointShape();
1293         if (getDataPointRank()==2) {
1294            DataArrayView::ShapeType ev_shape;
1295            Data ev(0.,ev_shape,getFunctionSpace());
1296            ev.typeMatchRight(*this);
1297            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1298            return ev;
1299         }
1300         if (getDataPointRank()==3) {
1301            DataArrayView::ShapeType ev_shape;
1302            if (axis_offset==0) {
1303              int s2=s[2];
1304              ev_shape.push_back(s2);
1305            }
1306            else if (axis_offset==1) {
1307              int s0=s[0];
1308              ev_shape.push_back(s0);
1309            }
1310            Data ev(0.,ev_shape,getFunctionSpace());
1311            ev.typeMatchRight(*this);
1312            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1313            return ev;
1314         }
1315         if (getDataPointRank()==4) {
1316            DataArrayView::ShapeType ev_shape;
1317            if (axis_offset==0) {
1318              ev_shape.push_back(s[2]);
1319              ev_shape.push_back(s[3]);
1320            }
1321            else if (axis_offset==1) {
1322              ev_shape.push_back(s[0]);
1323              ev_shape.push_back(s[3]);
1324            }
1325        else if (axis_offset==2) {
1326          ev_shape.push_back(s[0]);
1327          ev_shape.push_back(s[1]);
1328        }
1329            Data ev(0.,ev_shape,getFunctionSpace());
1330            ev.typeMatchRight(*this);
1331        m_data->matrixtrace(ev.m_data.get(), axis_offset);
1332            return ev;
1333         }
1334         else {
1335            throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1336         }
1337    }
1338    
1339    Data
1340    Data::transpose(int axis_offset) const
1341    {
1342    #if defined DOPROF
1343         profData->reduction2++;
1344    #endif
1345         DataArrayView::ShapeType s=getDataPointShape();
1346         DataArrayView::ShapeType ev_shape;
1347         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1348         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1349         int rank=getDataPointRank();
1350         if (axis_offset<0 || axis_offset>rank) {
1351            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1352         }
1353         for (int i=0; i<rank; i++) {
1354           int index = (axis_offset+i)%rank;
1355           ev_shape.push_back(s[index]); // Append to new shape
1356         }
1357         Data ev(0.,ev_shape,getFunctionSpace());
1358         ev.typeMatchRight(*this);
1359         m_data->transpose(ev.m_data.get(), axis_offset);
1360         return ev;
1361    }
1362    
1363    Data
1364    Data::eigenvalues() const
1365    {
1366         #if defined DOPROF
1367            profData->unary++;
1368         #endif
1369         // check input
1370         DataArrayView::ShapeType s=getDataPointShape();
1371         if (getDataPointRank()!=2)
1372            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1373         if(s[0] != s[1])
1374            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1375         // create return
1376         DataArrayView::ShapeType ev_shape(1,s[0]);
1377         Data ev(0.,ev_shape,getFunctionSpace());
1378         ev.typeMatchRight(*this);
1379         m_data->eigenvalues(ev.m_data.get());
1380         return ev;
1381    }
1382    
1383    const boost::python::tuple
1384    Data::eigenvalues_and_eigenvectors(const double tol) const
1385    {
1386         #if defined DOPROF
1387            profData->unary++;
1388         #endif
1389         DataArrayView::ShapeType s=getDataPointShape();
1390         if (getDataPointRank()!=2)
1391            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1392         if(s[0] != s[1])
1393            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1394         // create return
1395         DataArrayView::ShapeType ev_shape(1,s[0]);
1396         Data ev(0.,ev_shape,getFunctionSpace());
1397         ev.typeMatchRight(*this);
1398         DataArrayView::ShapeType V_shape(2,s[0]);
1399         Data V(0.,V_shape,getFunctionSpace());
1400         V.typeMatchRight(*this);
1401         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1402         return make_tuple(boost::python::object(ev),boost::python::object(V));
1403  }  }
1404    
1405  const boost::python::tuple  const boost::python::tuple
1406  Data::mindp() const  Data::mindp() const
1407  {  {
1408      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1409      // abort (for unknown reasons) if there are openmp directives with it in the
1410      // surrounding function
1411    
1412      int SampleNo;
1413      int DataPointNo;
1414    #ifdef PASO_MPI
1415        int ProcNo;
1416      calc_mindp(ProcNo,SampleNo,DataPointNo);
1417      return make_tuple(ProcNo,SampleNo,DataPointNo);
1418    #else
1419      calc_mindp(SampleNo,DataPointNo);
1420      return make_tuple(SampleNo,DataPointNo);
1421    #endif
1422    }
1423    
1424    void
1425    #ifdef PASO_MPI
1426    Data::calc_mindp(   int& ProcNo,
1427                    int& SampleNo,
1428            int& DataPointNo) const
1429    #else
1430    Data::calc_mindp(   int& SampleNo,
1431            int& DataPointNo) const
1432    #endif
1433    {
1434      int i,j;
1435      int lowi=0,lowj=0;
1436      double min=numeric_limits<double>::max();
1437    
1438    Data temp=minval();    Data temp=minval();
1439    
1440    int numSamples=temp.getNumSamples();    int numSamples=temp.getNumSamples();
1441    int numDPPSample=temp.getNumDataPointsPerSample();    int numDPPSample=temp.getNumDataPointsPerSample();
1442    
1443    int i,j,lowi=0,lowj=0;    double next,local_min;
1444    double min=numeric_limits<double>::max();    int local_lowi,local_lowj;
1445    
1446    for (i=0; i<numSamples; i++) {    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1447      for (j=0; j<numDPPSample; j++) {    {
1448        double next=temp.getDataPoint(i,j)();      local_min=min;
1449        if (next<min) {      #pragma omp for private(i,j) schedule(static)
1450          min=next;      for (i=0; i<numSamples; i++) {
1451          lowi=i;        for (j=0; j<numDPPSample; j++) {
1452          lowj=j;          next=temp.getDataPoint(i,j)();
1453            if (next<local_min) {
1454              local_min=next;
1455              local_lowi=i;
1456              local_lowj=j;
1457            }
1458        }        }
1459      }      }
1460        #pragma omp critical
1461        if (local_min<min) {
1462          min=local_min;
1463          lowi=local_lowi;
1464          lowj=local_lowj;
1465        }
1466    }    }
1467    
1468    return make_tuple(lowi,lowj);  #ifdef PASO_MPI
1469  }      // determine the processor on which the minimum occurs
1470        next = temp.getDataPoint(lowi,lowj)();
1471  Data      int lowProc = 0;
1472  Data::length() const      double *globalMins = new double[get_MPISize()+1];
1473  {      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1474    return dp_algorithm(DataAlgorithmAdapter<Length>(0));      
1475  }      if( get_MPIRank()==0 ){
1476            next = globalMins[lowProc];
1477  Data          for( i=1; i<get_MPISize(); i++ )
1478  Data::trace() const              if( next>globalMins[i] ){
1479  {                  lowProc = i;
1480    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));                  next = globalMins[i];
1481  }              }
1482        }
1483  Data      MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1484  Data::transpose(int axis) const  
1485  {      delete [] globalMins;
1486    // not implemented      ProcNo = lowProc;
1487    throw DataException("Error - Data::transpose not implemented yet.");  #endif
1488    return Data();    SampleNo = lowi;
1489      DataPointNo = lowj;
1490  }  }
1491    
1492  void  void
1493  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1494  {  {
1495    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1496      args["data"]=boost::python::object(this);
1497      getDomain().saveDX(fileName,args);
1498    return;    return;
1499  }  }
1500    
1501  void  void
1502  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1503  {  {
1504    getDomain().saveVTK(fileName,*this);    boost::python::dict args;
1505      args["data"]=boost::python::object(this);
1506      getDomain().saveVTK(fileName,args);
1507    return;    return;
1508  }  }
1509    
1510  Data&  Data&
1511  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1512  {  {
1513    #if defined DOPROF
1514      profData->binary++;
1515    #endif
1516    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1517    return (*this);    return (*this);
1518  }  }
# Line 991  Data::operator+=(const Data& right) Line 1520  Data::operator+=(const Data& right)
1520  Data&  Data&
1521  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1522  {  {
1523    #if defined DOPROF
1524      profData->binary++;
1525    #endif
1526    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1527    return (*this);    return (*this);
1528  }  }
# Line 998  Data::operator+=(const boost::python::ob Line 1530  Data::operator+=(const boost::python::ob
1530  Data&  Data&
1531  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1532  {  {
1533    #if defined DOPROF
1534      profData->binary++;
1535    #endif
1536    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1537    return (*this);    return (*this);
1538  }  }
# Line 1005  Data::operator-=(const Data& right) Line 1540  Data::operator-=(const Data& right)
1540  Data&  Data&
1541  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1542  {  {
1543    #if defined DOPROF
1544      profData->binary++;
1545    #endif
1546    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1547    return (*this);    return (*this);
1548  }  }
# Line 1012  Data::operator-=(const boost::python::ob Line 1550  Data::operator-=(const boost::python::ob
1550  Data&  Data&
1551  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1552  {  {
1553    #if defined DOPROF
1554      profData->binary++;
1555    #endif
1556    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1557    return (*this);    return (*this);
1558  }  }
# Line 1019  Data::operator*=(const Data& right) Line 1560  Data::operator*=(const Data& right)
1560  Data&  Data&
1561  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1562  {  {
1563    #if defined DOPROF
1564      profData->binary++;
1565    #endif
1566    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1567    return (*this);    return (*this);
1568  }  }
# Line 1026  Data::operator*=(const boost::python::ob Line 1570  Data::operator*=(const boost::python::ob
1570  Data&  Data&
1571  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1572  {  {
1573    #if defined DOPROF
1574      profData->binary++;
1575    #endif
1576    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1577    return (*this);    return (*this);
1578  }  }
# Line 1033  Data::operator/=(const Data& right) Line 1580  Data::operator/=(const Data& right)
1580  Data&  Data&
1581  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1582  {  {
1583    #if defined DOPROF
1584      profData->binary++;
1585    #endif
1586    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1587    return (*this);    return (*this);
1588  }  }
1589    
1590  Data  Data
1591    Data::rpowO(const boost::python::object& left) const
1592    {
1593    #if defined DOPROF
1594      profData->binary++;
1595    #endif
1596      Data left_d(left,*this);
1597      return left_d.powD(*this);
1598    }
1599    
1600    Data
1601  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1602  {  {
1603    #if defined DOPROF
1604      profData->binary++;
1605    #endif
1606    Data result;    Data result;
1607    result.copy(*this);    result.copy(*this);
1608    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 1049  Data::powO(const boost::python::object& Line 1612  Data::powO(const boost::python::object&
1612  Data  Data
1613  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1614  {  {
1615    #if defined DOPROF
1616      profData->binary++;
1617    #endif
1618    Data result;    Data result;
1619    result.copy(*this);    result.copy(*this);
1620    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1621    return result;    return result;
1622  }  }
1623    
1624    
1625  //  //
1626  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1627  Data  Data
1628  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1629  {  {
# Line 1069  escript::operator+(const Data& left, con Line 1636  escript::operator+(const Data& left, con
1636  }  }
1637    
1638  //  //
1639  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1640  Data  Data
1641  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1642  {  {
# Line 1082  escript::operator-(const Data& left, con Line 1649  escript::operator-(const Data& left, con
1649  }  }
1650    
1651  //  //
1652  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1653  Data  Data
1654  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1655  {  {
# Line 1095  escript::operator*(const Data& left, con Line 1662  escript::operator*(const Data& left, con
1662  }  }
1663    
1664  //  //
1665  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1666  Data  Data
1667  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1668  {  {
# Line 1108  escript::operator/(const Data& left, con Line 1675  escript::operator/(const Data& left, con
1675  }  }
1676    
1677  //  //
1678  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1679  Data  Data
1680  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1681  {  {
# Line 1124  escript::operator+(const Data& left, con Line 1691  escript::operator+(const Data& left, con
1691  }  }
1692    
1693  //  //
1694  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1695  Data  Data
1696  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1697  {  {
# Line 1140  escript::operator-(const Data& left, con Line 1707  escript::operator-(const Data& left, con
1707  }  }
1708    
1709  //  //
1710  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1711  Data  Data
1712  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1713  {  {
# Line 1156  escript::operator*(const Data& left, con Line 1723  escript::operator*(const Data& left, con
1723  }  }
1724    
1725  //  //
1726  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1727  Data  Data
1728  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1729  {  {
# Line 1172  escript::operator/(const Data& left, con Line 1739  escript::operator/(const Data& left, con
1739  }  }
1740    
1741  //  //
1742  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1743  Data  Data
1744  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1745  {  {
# Line 1185  escript::operator+(const boost::python:: Line 1752  escript::operator+(const boost::python::
1752  }  }
1753    
1754  //  //
1755  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1756  Data  Data
1757  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1758  {  {
# Line 1198  escript::operator-(const boost::python:: Line 1765  escript::operator-(const boost::python::
1765  }  }
1766    
1767  //  //
1768  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1769  Data  Data
1770  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1771  {  {
# Line 1211  escript::operator*(const boost::python:: Line 1778  escript::operator*(const boost::python::
1778  }  }
1779    
1780  //  //
1781  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1782  Data  Data
1783  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1784  {  {
# Line 1224  escript::operator/(const boost::python:: Line 1791  escript::operator/(const boost::python::
1791  }  }
1792    
1793  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1794  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1795  //{  //{
1796  //  /*  //  /*
# Line 1269  escript::operator/(const boost::python:: Line 1835  escript::operator/(const boost::python::
1835  //  return ret;  //  return ret;
1836  //}  //}
1837    
1838    /* TODO */
1839    /* global reduction */
1840  Data  Data
1841  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1842  {  {
# Line 1283  Data::getItem(const boost::python::objec Line 1851  Data::getItem(const boost::python::objec
1851    return getSlice(slice_region);    return getSlice(slice_region);
1852  }  }
1853    
1854    /* TODO */
1855    /* global reduction */
1856  Data  Data
1857  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1858  {  {
1859    #if defined DOPROF
1860      profData->slicing++;
1861    #endif
1862    return Data(*this,region);    return Data(*this,region);
1863  }  }
1864    
1865    /* TODO */
1866    /* global reduction */
1867  void  void
1868  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1869                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1297  Data::setItemO(const boost::python::obje Line 1872  Data::setItemO(const boost::python::obje
1872    setItemD(key,tempData);    setItemD(key,tempData);
1873  }  }
1874    
1875    /* TODO */
1876    /* global reduction */
1877  void  void
1878  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1879                 const Data& value)                 const Data& value)
1880  {  {
1881    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1882    
1883    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1884    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1885      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 1891  Data::setItemD(const boost::python::obje
1891    }    }
1892  }  }
1893    
1894    /* TODO */
1895    /* global reduction */
1896  void  void
1897  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1898                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1899  {  {
1900    #if defined DOPROF
1901      profData->slicing++;
1902    #endif
1903    Data tempValue(value);    Data tempValue(value);
1904    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1905    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1351  Data::typeMatchRight(const Data& right) Line 1934  Data::typeMatchRight(const Data& right)
1934    }    }
1935  }  }
1936    
1937    /* TODO */
1938    /* global reduction */
1939  void  void
1940  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1941                       const boost::python::object& value)                       const boost::python::object& value)
# Line 1372  Data::setTaggedValue(int tagKey, Line 1957  Data::setTaggedValue(int tagKey,
1957    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1958  }  }
1959    
1960    /* TODO */
1961    /* global reduction */
1962  void  void
1963  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1964                              const DataArrayView& value)                              const DataArrayView& value)
# Line 1389  Data::setTaggedValueFromCPP(int tagKey, Line 1976  Data::setTaggedValueFromCPP(int tagKey,
1976    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1977  }  }
1978    
1979    /* TODO */
1980    /* global reduction */
1981    int
1982    Data::getTagNumber(int dpno)
1983    {
1984      return m_data->getTagNumber(dpno);
1985    }
1986    
1987    /* TODO */
1988    /* global reduction */
1989  void  void
1990  Data::setRefValue(int ref,  Data::setRefValue(int ref,
1991                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
# Line 1402  Data::setRefValue(int ref, Line 1999  Data::setRefValue(int ref,
1999    m_data->setRefValue(ref,valueDataArray);    m_data->setRefValue(ref,valueDataArray);
2000  }  }
2001    
2002    /* TODO */
2003    /* global reduction */
2004  void  void
2005  Data::getRefValue(int ref,  Data::getRefValue(int ref,
2006                    boost::python::numeric::array& value)                    boost::python::numeric::array& value)
# Line 1428  Data::getRefValue(int ref, Line 2027  Data::getRefValue(int ref,
2027    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
2028    
2029    if (rank==0) {    if (rank==0) {
2030      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
2031          value = temp_numArray;
2032    }    }
2033    if (rank==1) {    if (rank==1) {
2034      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1436  Data::getRefValue(int ref, Line 2036  Data::getRefValue(int ref,
2036      }      }
2037    }    }
2038    if (rank==2) {    if (rank==2) {
2039      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2040          for (int j=0; j < shape[1]; j++) {
2041            value[i][j] = valueView(i,j);
2042          }
2043        }
2044    }    }
2045    if (rank==3) {    if (rank==3) {
2046      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2047          for (int j=0; j < shape[1]; j++) {
2048            for (int k=0; k < shape[2]; k++) {
2049              value[i][j][k] = valueView(i,j,k);
2050            }
2051          }
2052        }
2053    }    }
2054    if (rank==4) {    if (rank==4) {
2055      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2056          for (int j=0; j < shape[1]; j++) {
2057            for (int k=0; k < shape[2]; k++) {
2058              for (int l=0; l < shape[3]; l++) {
2059                value[i][j][k][l] = valueView(i,j,k,l);
2060              }
2061            }
2062          }
2063        }
2064    }    }
2065    
2066  }  }
# Line 1472  Data::archiveData(const std::string file Line 2090  Data::archiveData(const std::string file
2090      dataType = 3;      dataType = 3;
2091      cout << "\tdataType: DataExpanded" << endl;      cout << "\tdataType: DataExpanded" << endl;
2092    }    }
2093    
2094    if (dataType == -1) {    if (dataType == -1) {
2095      throw DataException("archiveData Error: undefined dataType");      throw DataException("archiveData Error: undefined dataType");
2096    }    }
# Line 1485  Data::archiveData(const std::string file Line 2104  Data::archiveData(const std::string file
2104    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2105    int dataLength = getLength();    int dataLength = getLength();
2106    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2107    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2108    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2109      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2110    }    }
2111    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2112    if (isTagged()) {    if (isTagged()) {
2113      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2114        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 1511  Data::archiveData(const std::string file Line 2130  Data::archiveData(const std::string file
2130    cout << ">" << endl;    cout << ">" << endl;
2131    
2132    //    //
2133    // Write common data items to archive file    // Open archive file
2134    ofstream archiveFile;    ofstream archiveFile;
2135    archiveFile.open(fileName.data(), ios::out);    archiveFile.open(fileName.data(), ios::out);
2136    
# Line 1519  Data::archiveData(const std::string file Line 2138  Data::archiveData(const std::string file
2138      throw DataException("archiveData Error: problem opening archive file");      throw DataException("archiveData Error: problem opening archive file");
2139    }    }
2140    
2141      //
2142      // Write common data items to archive file
2143    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2144    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2145    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 2163  Data::archiveData(const std::string file
2163      throw DataException("archiveData Error: problem writing to archive file");      throw DataException("archiveData Error: problem writing to archive file");
2164    }    }
2165    
   archiveFile.close();  
   
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
   }  
   
2166    //    //
2167    // Collect and archive underlying data values for each Data type    // Archive underlying data values for each Data type
2168      int noValues;
2169    switch (dataType) {    switch (dataType) {
2170      case 0:      case 0:
2171        // DataEmpty        // DataEmpty
2172          noValues = 0;
2173          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2174          cout << "\tnoValues: " << noValues << endl;
2175        break;        break;
2176      case 1:      case 1:
2177        // DataConstant        // DataConstant
2178          noValues = m_data->getLength();
2179          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2180          cout << "\tnoValues: " << noValues << endl;
2181          if (m_data->archiveData(archiveFile,noValues)) {
2182            throw DataException("archiveData Error: problem writing data to archive file");
2183          }
2184        break;        break;
2185      case 2:      case 2:
2186        // DataTagged        // DataTagged
2187          noValues = m_data->getLength();
2188          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2189          cout << "\tnoValues: " << noValues << endl;
2190          if (m_data->archiveData(archiveFile,noValues)) {
2191            throw DataException("archiveData Error: problem writing data to archive file");
2192          }
2193        break;        break;
2194      case 3:      case 3:
2195        // DataExpanded        // DataExpanded
2196          noValues = m_data->getLength();
2197          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2198          cout << "\tnoValues: " << noValues << endl;
2199          if (m_data->archiveData(archiveFile,noValues)) {
2200            throw DataException("archiveData Error: problem writing data to archive file");
2201          }
2202        break;        break;
2203    }    }
2204    
2205      if (!archiveFile.good()) {
2206        throw DataException("archiveData Error: problem writing data to archive file");
2207      }
2208    
2209      //
2210      // Close archive file
2211      archiveFile.close();
2212    
2213      if (!archiveFile.good()) {
2214        throw DataException("archiveData Error: problem closing archive file");
2215      }
2216    
2217  }  }
2218    
2219  void  void
# Line 1590  Data::extractData(const std::string file Line 2239  Data::extractData(const std::string file
2239    int flatShape[4];    int flatShape[4];
2240    
2241    //    //
2242    // Open the archive file and read common data items    // Open the archive file
2243    ifstream archiveFile;    ifstream archiveFile;
2244    archiveFile.open(fileName.data(), ios::in);    archiveFile.open(fileName.data(), ios::in);
2245    
# Line 1598  Data::extractData(const std::string file Line 2247  Data::extractData(const std::string file
2247      throw DataException("extractData Error: problem opening archive file");      throw DataException("extractData Error: problem opening archive file");
2248    }    }
2249    
2250      //
2251      // Read common data items from archive file
2252    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2253    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2254    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 2262  Data::extractData(const std::string file
2262        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2263      }      }
2264    }    }
2265    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2266    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2267      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2268    }    }
2269    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2270    if (dataType==2) {    if (dataType==2) {
2271      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2272        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 2277  Data::extractData(const std::string file
2277      throw DataException("extractData Error: problem reading from archive file");      throw DataException("extractData Error: problem reading from archive file");
2278    }    }
2279    
2280    archiveFile.close();    //
2281      // Verify the values just read from the archive file
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem closing archive file");  
   }  
   
2282    switch (dataType) {    switch (dataType) {
2283      case 0:      case 0:
2284        cout << "\tdataType: DataEmpty" << endl;        cout << "\tdataType: DataEmpty" << endl;
# Line 1686  Data::extractData(const std::string file Line 2333  Data::extractData(const std::string file
2333    
2334    //    //
2335    // Load this DataVector with the appropriate values    // Load this DataVector with the appropriate values
2336      int noValues;
2337      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2338      cout << "\tnoValues: " << noValues << endl;
2339    switch (dataType) {    switch (dataType) {
2340      case 0:      case 0:
2341        // DataEmpty        // DataEmpty
2342          if (noValues != 0) {
2343            throw DataException("extractData Error: problem reading data from archive file");
2344          }
2345        break;        break;
2346      case 1:      case 1:
2347        // DataConstant        // DataConstant
2348          if (dataVec.extractData(archiveFile,noValues)) {
2349            throw DataException("extractData Error: problem reading data from archive file");
2350          }
2351        break;        break;
2352      case 2:      case 2:
2353        // DataTagged        // DataTagged
2354          if (dataVec.extractData(archiveFile,noValues)) {
2355            throw DataException("extractData Error: problem reading data from archive file");
2356          }
2357        break;        break;
2358      case 3:      case 3:
2359        // DataExpanded        // DataExpanded
2360          if (dataVec.extractData(archiveFile,noValues)) {
2361            throw DataException("extractData Error: problem reading data from archive file");
2362          }
2363        break;        break;
2364    }    }
2365    
2366      if (!archiveFile.good()) {
2367        throw DataException("extractData Error: problem reading from archive file");
2368      }
2369    
2370      //
2371      // Close archive file
2372      archiveFile.close();
2373    
2374      if (!archiveFile.good()) {
2375        throw DataException("extractData Error: problem closing archive file");
2376      }
2377    
2378    //    //
2379    // Construct an appropriate Data object    // Construct an appropriate Data object
2380    DataAbstract* tempData;    DataAbstract* tempData;
# Line 1724  Data::extractData(const std::string file Line 2398  Data::extractData(const std::string file
2398    }    }
2399    shared_ptr<DataAbstract> temp_data(tempData);    shared_ptr<DataAbstract> temp_data(tempData);
2400    m_data=temp_data;    m_data=temp_data;
   
2401  }  }
2402    
2403  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
# Line 1732  ostream& escript::operator<<(ostream& o, Line 2405  ostream& escript::operator<<(ostream& o,
2405    o << data.toString();    o << data.toString();
2406    return o;    return o;
2407  }  }
2408    
2409    /* Member functions specific to the MPI implementation */
2410    #ifdef PASO_MPI
2411    
2412    void
2413    Data::print()
2414    {
2415      int i,j;
2416      
2417      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2418      for( i=0; i<getNumSamples(); i++ )
2419      {
2420        printf( "[%6d]", i );
2421        for( j=0; j<getNumDataPointsPerSample(); j++ )
2422          printf( "\t%10.7g", (getSampleData(i))[j] );
2423        printf( "\n" );
2424      }
2425    }
2426    
2427    int
2428    Data::get_MPISize() const
2429    {
2430        int error, size;
2431        error = MPI_Comm_size( get_MPIComm(), &size );
2432    
2433        return size;
2434    }
2435    
2436    int
2437    Data::get_MPIRank() const
2438    {
2439        int error, rank;
2440        error = MPI_Comm_rank( get_MPIComm(), &rank );
2441    
2442        return rank;
2443    }
2444    
2445    MPI_Comm
2446    Data::get_MPIComm() const
2447    {
2448        return MPI_COMM_WORLD;
2449    }
2450    #endif

Legend:
Removed from v.121  
changed lines
  Added in v.782

  ViewVC Help
Powered by ViewVC 1.1.26