/[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/escript/src/Data/Data.cpp revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC trunk/escript/src/Data.cpp revision 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
2    
3   ******************************************************************************  /*
4   *                                                                            *   ************************************************************
5   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *   *          Copyright 2006 by ACcESS MNRF                   *
6   *                                                                            *   *                                                          *
7   * This software is the property of ACcESS.  No part of this code             *   *              http://www.access.edu.au                    *
8   * may be copied in any form or by any means without the expressed written    *   *       Primary Business: Queensland, Australia            *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *  Licensed under the Open Software License version 3.0    *
10   * by any unauthorised person is illegal unless that                          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
11   * person has a software license agreement with ACcESS.                       *   *                                                          *
12   *                                                                            *   ************************************************************
13   ******************************************************************************  */
14    #include "Data.h"
15    
16    #include "DataExpanded.h"
17    #include "DataConstant.h"
18    #include "DataTagged.h"
19    #include "DataEmpty.h"
20    #include "DataArray.h"
21    #include "DataArrayView.h"
22    #include "DataProf.h"
23    #include "FunctionSpaceFactory.h"
24    #include "AbstractContinuousDomain.h"
25    #include "UnaryFuncs.h"
26    
 ******************************************************************************/  
   
 #include "escript/Data/Data.h"  
   
 #include <iostream>  
27  #include <fstream>  #include <fstream>
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
 #include <exception>  
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
 #include <boost/python/str.hpp>  
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
34  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
35    
 #include "escript/Data/DataException.h"  
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataEmpty.h"  
 #include "escript/Data/DataArray.h"  
 #include "escript/Data/DataProf.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/AbstractContinuousDomain.h"  
 #include "escript/Data/UnaryFuncs.h"  
   
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
# Line 59  Data::Data() Line 51  Data::Data()
51    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
52    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
53    m_data=temp_data;    m_data=temp_data;
54      m_protected=false;
55  #if defined DOPROF  #if defined DOPROF
56    // create entry in global profiling table for this object    // create entry in global profiling table for this object
57    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 76  Data::Data(double value, Line 69  Data::Data(double value,
69    }    }
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
72      m_protected=false;
73  #if defined DOPROF  #if defined DOPROF
74    // create entry in global profiling table for this object    // create entry in global profiling table for this object
75    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 90  Data::Data(double value, Line 84  Data::Data(double value,
84    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
85    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
86    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
87      m_protected=false;
88  #if defined DOPROF  #if defined DOPROF
89    // create entry in global profiling table for this object    // create entry in global profiling table for this object
90    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 99  Data::Data(double value, Line 94  Data::Data(double value,
94  Data::Data(const Data& inData)  Data::Data(const Data& inData)
95  {  {
96    m_data=inData.m_data;    m_data=inData.m_data;
97      m_protected=inData.isProtected();
98  #if defined DOPROF  #if defined DOPROF
99    // create entry in global profiling table for this object    // create entry in global profiling table for this object
100    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 113  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
110    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
111    m_data=temp_data;    m_data=temp_data;
112      m_protected=false;
113  #if defined DOPROF  #if defined DOPROF
114    // create entry in global profiling table for this object    // create entry in global profiling table for this object
115    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 122  Data::Data(const Data& inData, Line 119  Data::Data(const Data& inData,
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122    #if defined DOPROF
123      // create entry in global profiling table for this object
124      profData = dataProfTable.newData();
125    #endif
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    } else {
129        #if defined DOPROF
130        profData->interpolate++;
131        #endif
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
133      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
134      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 138  Data::Data(const Data& inData, Line 142  Data::Data(const Data& inData,
142      }      }
143      m_data=tmp.m_data;      m_data=tmp.m_data;
144    }    }
145  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 153  Data::Data(const DataTagged::TagListType Line 154  Data::Data(const DataTagged::TagListType
154    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
155    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
156    m_data=temp_data;    m_data=temp_data;
157      m_protected=false;
158    if (expanded) {    if (expanded) {
159      expand();      expand();
160    }    }
# Line 167  Data::Data(const numeric::array& value, Line 169  Data::Data(const numeric::array& value,
169             bool expanded)             bool expanded)
170  {  {
171    initialise(value,what,expanded);    initialise(value,what,expanded);
172      m_protected=false;
173  #if defined DOPROF  #if defined DOPROF
174    // create entry in global profiling table for this object    // create entry in global profiling table for this object
175    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 178  Data::Data(const DataArrayView& value, Line 181  Data::Data(const DataArrayView& value,
181             bool expanded)             bool expanded)
182  {  {
183    initialise(value,what,expanded);    initialise(value,what,expanded);
184      m_protected=false;
185  #if defined DOPROF  #if defined DOPROF
186    // create entry in global profiling table for this object    // create entry in global profiling table for this object
187    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 190  Data::Data(const object& value, Line 194  Data::Data(const object& value,
194  {  {
195    numeric::array asNumArray(value);    numeric::array asNumArray(value);
196    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
197      m_protected=false;
198  #if defined DOPROF  #if defined DOPROF
199    // create entry in global profiling table for this object    // create entry in global profiling table for this object
200    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 214  Data::Data(const object& value, Line 219  Data::Data(const object& value,
219      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
220      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
221    }    }
222      m_protected=false;
223  #if defined DOPROF  #if defined DOPROF
224    // create entry in global profiling table for this object    // create entry in global profiling table for this object
225    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 344  Data::isTagged() const Line 350  Data::isTagged() const
350    return (temp!=0);    return (temp!=0);
351  }  }
352    
353    /* TODO */
354    /* global reduction -- the local data being empty does not imply that it is empty on other processers*/
355  bool  bool
356  Data::isEmpty() const  Data::isEmpty() const
357  {  {
# Line 359  Data::isConstant() const Line 367  Data::isConstant() const
367  }  }
368    
369  void  void
370    Data::setProtection()
371    {
372       m_protected=true;
373    }
374    
375    bool
376    Data::isProtected() const
377    {
378       return m_protected;
379    }
380    
381    
382    
383    void
384  Data::expand()  Data::expand()
385  {  {
386    if (isConstant()) {    if (isConstant()) {
# Line 443  Data::whereNonPositive() const Line 465  Data::whereNonPositive() const
465  }  }
466    
467  Data  Data
468  Data::whereZero() const  Data::whereZero(double tol) const
469  {  {
470  #if defined DOPROF  #if defined DOPROF
471    profData->where++;    profData->where++;
472  #endif  #endif
473    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
474      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
475  }  }
476    
477  Data  Data
478  Data::whereNonZero() const  Data::whereNonZero(double tol) const
479  {  {
480  #if defined DOPROF  #if defined DOPROF
481    profData->where++;    profData->where++;
482  #endif  #endif
483    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
484      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
485  }  }
486    
487  Data  Data
# Line 526  Data::getDataPointShape() const Line 550  Data::getDataPointShape() const
550  void  void
551  Data::fillFromNumArray(const boost::python::numeric::array num_array)  Data::fillFromNumArray(const boost::python::numeric::array num_array)
552  {  {
553      if (isProtected()) {
554            throw DataException("Error - attempt to update protected Data object.");
555      }
556    //    //
557    // check rank    // check rank
558    if (num_array.getrank()<getDataPointRank())    if (num_array.getrank()<getDataPointRank())
# Line 755  Data::convertToNumArrayFromSampleNo(int Line 782  Data::convertToNumArrayFromSampleNo(int
782    
783  const  const
784  boost::python::numeric::array  boost::python::numeric::array
785  Data::convertToNumArrayFromDPNo(int sampleNo,  Data::convertToNumArrayFromDPNo(int procNo,
786                                  int dataPointNo)                                  int sampleNo,
787                                                                    int dataPointNo)
788    
789  {  {
790        size_t length=0;
791        int i, j, k, l, pos;
792    
793    //    //
794    // Check a valid sample number has been supplied    // Check a valid sample number has been supplied
795    if (sampleNo >= getNumSamples()) {    if (sampleNo >= getNumSamples()) {
# Line 803  Data::convertToNumArrayFromDPNo(int samp Line 835  Data::convertToNumArrayFromDPNo(int samp
835      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
836    }    }
837    
838        // added for the MPI communication
839        length=1;
840        for( i=0; i<arrayRank; i++ )
841            length *= arrayShape[i];
842        double *tmpData = new double[length];
843    
844    //    //
845    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
846    DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
847        // updated for the MPI case
848        if( get_MPIRank()==procNo ){
849            // create a view of the data if it is stored locally
850            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
851            
852            // pack the data from the view into tmpData for MPI communication
853            pos=0;
854            switch( dataPointRank ){
855                case 0 :
856                    tmpData[0] = dataPointView();
857                    break;
858                case 1 :        
859                    for( i=0; i<dataPointShape[0]; i++ )
860                        tmpData[i]=dataPointView(i);
861                    break;
862                case 2 :        
863                    for( i=0; i<dataPointShape[0]; i++ )
864                        for( j=0; j<dataPointShape[1]; j++, pos++ )
865                            tmpData[pos]=dataPointView(i,j);
866                    break;
867                case 3 :        
868                    for( i=0; i<dataPointShape[0]; i++ )
869                        for( j=0; j<dataPointShape[1]; j++ )
870                            for( k=0; k<dataPointShape[2]; k++, pos++ )
871                                tmpData[pos]=dataPointView(i,j,k);
872                    break;
873                case 4 :
874                    for( i=0; i<dataPointShape[0]; i++ )
875                        for( j=0; j<dataPointShape[1]; j++ )
876                            for( k=0; k<dataPointShape[2]; k++ )
877                                for( l=0; l<dataPointShape[3]; l++, pos++ )
878                                    tmpData[pos]=dataPointView(i,j,k,l);
879                    break;
880            }
881        }
882    #ifdef PASO_MPI
883            // broadcast the data to all other processes
884            MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
885    #endif
886    
887        // unpack the data
888        switch( dataPointRank ){
889            case 0 :
890                numArray[i]=tmpData[0];
891                break;
892            case 1 :        
893                for( i=0; i<dataPointShape[0]; i++ )
894                    numArray[i]=tmpData[i];
895                break;
896            case 2 :        
897                for( i=0; i<dataPointShape[0]; i++ )
898                    for( j=0; j<dataPointShape[1]; j++ )
899                        tmpData[i+j*dataPointShape[0]];
900                break;
901            case 3 :        
902                for( i=0; i<dataPointShape[0]; i++ )
903                    for( j=0; j<dataPointShape[1]; j++ )
904                        for( k=0; k<dataPointShape[2]; k++ )
905                            tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
906                break;
907            case 4 :
908                for( i=0; i<dataPointShape[0]; i++ )
909                    for( j=0; j<dataPointShape[1]; j++ )
910                        for( k=0; k<dataPointShape[2]; k++ )
911                            for( l=0; l<dataPointShape[3]; l++ )
912                                tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
913                break;
914        }
915    
916        delete [] tmpData;  
917    /*
918    if (dataPointRank==0) {    if (dataPointRank==0) {
919      numArray[0]=dataPointView();      numArray[0]=dataPointView();
920    }    }
# Line 841  Data::convertToNumArrayFromDPNo(int samp Line 950  Data::convertToNumArrayFromDPNo(int samp
950        }        }
951      }      }
952    }    }
953    */
954    
955    //    //
956    // return the loaded array    // return the loaded array
# Line 880  Data::integrate() const Line 990  Data::integrate() const
990      }      }
991    }    }
992    if (rank==2) {    if (rank==2) {
993      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
994      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
995        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
996          index = i + shape[0] * j;             index = i + shape[0] * j;
997          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
998        }           }
999      }         }
1000    }    }
1001    if (rank==3) {    if (rank==3) {
1002      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 894  Data::integrate() const Line 1004  Data::integrate() const
1004        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
1005          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1006            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
1007            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
1008          }          }
1009        }        }
1010      }      }
# Line 906  Data::integrate() const Line 1016  Data::integrate() const
1016          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1017            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
1018              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1019              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
1020            }            }
1021          }          }
1022        }        }
# Line 1027  Data::atanh() const Line 1137  Data::atanh() const
1137  }  }
1138    
1139  Data  Data
1140  Data::log() const  Data::log10() const
1141  {  {
1142  #if defined DOPROF  #if defined DOPROF
1143    profData->unary++;    profData->unary++;
# Line 1036  Data::log() const Line 1146  Data::log() const
1146  }  }
1147    
1148  Data  Data
1149  Data::ln() const  Data::log() const
1150  {  {
1151  #if defined DOPROF  #if defined DOPROF
1152    profData->unary++;    profData->unary++;
# Line 1104  Data::sqrt() const Line 1214  Data::sqrt() const
1214  double  double
1215  Data::Lsup() const  Data::Lsup() const
1216  {  {
1217      double localValue, globalValue;
1218  #if defined DOPROF  #if defined DOPROF
1219    profData->reduction1++;    profData->reduction1++;
1220  #endif  #endif
1221    //    //
1222    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1223    
1224    AbsMax abs_max_func;    AbsMax abs_max_func;
1225    return algorithm(abs_max_func,0);    localValue = algorithm(abs_max_func,0);
1226    #ifdef PASO_MPI
1227      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1228      return globalValue;
1229    #else
1230      return localValue;
1231    #endif
1232  }  }
1233    
1234  double  double
1235  Data::Linf() const  Data::Linf() const
1236  {  {
1237      double localValue, globalValue;
1238  #if defined DOPROF  #if defined DOPROF
1239    profData->reduction1++;    profData->reduction1++;
1240  #endif  #endif
1241    //    //
1242    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1243    AbsMin abs_min_func;    AbsMin abs_min_func;
1244    return algorithm(abs_min_func,numeric_limits<double>::max());    localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1245    
1246    #ifdef PASO_MPI
1247      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1248      return globalValue;
1249    #else
1250      return localValue;
1251    #endif
1252  }  }
1253    
1254  double  double
1255  Data::sup() const  Data::sup() const
1256  {  {
1257      double localValue, globalValue;
1258  #if defined DOPROF  #if defined DOPROF
1259    profData->reduction1++;    profData->reduction1++;
1260  #endif  #endif
1261    //    //
1262    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1263    FMax fmax_func;    FMax fmax_func;
1264    return algorithm(fmax_func,numeric_limits<double>::max()*-1);    localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1265    #ifdef PASO_MPI
1266      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1267      return globalValue;
1268    #else
1269      return localValue;
1270    #endif
1271  }  }
1272    
1273  double  double
1274  Data::inf() const  Data::inf() const
1275  {  {
1276      double localValue, globalValue;
1277  #if defined DOPROF  #if defined DOPROF
1278    profData->reduction1++;    profData->reduction1++;
1279  #endif  #endif
1280    //    //
1281    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1282    FMin fmin_func;    FMin fmin_func;
1283    return algorithm(fmin_func,numeric_limits<double>::max());    localValue = algorithm(fmin_func,numeric_limits<double>::max());
1284    #ifdef PASO_MPI
1285      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1286      return globalValue;
1287    #else
1288      return localValue;
1289    #endif
1290  }  }
1291    
1292    /* TODO */
1293    /* global reduction */
1294  Data  Data
1295  Data::maxval() const  Data::maxval() const
1296  {  {
# Line 1174  Data::minval() const Line 1316  Data::minval() const
1316  }  }
1317    
1318  Data  Data
1319  Data::length() const  Data::trace() const
1320  {  {
1321  #if defined DOPROF  #if defined DOPROF
1322    profData->reduction2++;    profData->reduction2++;
1323  #endif  #endif
1324    Length len_func;    Trace trace_func;
1325    return dp_algorithm(len_func,0);    return dp_algorithm(trace_func,0);
1326  }  }
1327    
1328  Data  Data
1329  Data::trace() const  Data::symmetric() const
1330  {  {
1331  #if defined DOPROF       #if defined DOPROF
1332    profData->reduction2++;          profData->unary++;
1333  #endif       #endif
1334    Trace trace_func;       // check input
1335    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1336         if (getDataPointRank()==2) {
1337            if(s[0] != s[1])
1338               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1339         }
1340         else if (getDataPointRank()==4) {
1341            if(!(s[0] == s[2] && s[1] == s[3]))
1342               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1343         }
1344         else {
1345            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1346         }
1347         Data ev(0.,getDataPointShape(),getFunctionSpace());
1348         ev.typeMatchRight(*this);
1349         m_data->symmetric(ev.m_data.get());
1350         return ev;
1351    }
1352    
1353    Data
1354    Data::nonsymmetric() const
1355    {
1356         #if defined DOPROF
1357            profData->unary++;
1358         #endif
1359         // check input
1360         DataArrayView::ShapeType s=getDataPointShape();
1361         if (getDataPointRank()==2) {
1362            if(s[0] != s[1])
1363               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1364            DataArrayView::ShapeType ev_shape;
1365            ev_shape.push_back(s[0]);
1366            ev_shape.push_back(s[1]);
1367            Data ev(0.,ev_shape,getFunctionSpace());
1368            ev.typeMatchRight(*this);
1369            m_data->nonsymmetric(ev.m_data.get());
1370            return ev;
1371         }
1372         else if (getDataPointRank()==4) {
1373            if(!(s[0] == s[2] && s[1] == s[3]))
1374               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1375            DataArrayView::ShapeType ev_shape;
1376            ev_shape.push_back(s[0]);
1377            ev_shape.push_back(s[1]);
1378            ev_shape.push_back(s[2]);
1379            ev_shape.push_back(s[3]);
1380            Data ev(0.,ev_shape,getFunctionSpace());
1381            ev.typeMatchRight(*this);
1382            m_data->nonsymmetric(ev.m_data.get());
1383            return ev;
1384         }
1385         else {
1386            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1387         }
1388    }
1389    
1390    Data
1391    Data::matrixtrace(int axis_offset) const
1392    {
1393         #if defined DOPROF
1394            profData->unary++;
1395         #endif
1396         DataArrayView::ShapeType s=getDataPointShape();
1397         if (getDataPointRank()==2) {
1398            DataArrayView::ShapeType ev_shape;
1399            Data ev(0.,ev_shape,getFunctionSpace());
1400            ev.typeMatchRight(*this);
1401            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1402            return ev;
1403         }
1404         if (getDataPointRank()==3) {
1405            DataArrayView::ShapeType ev_shape;
1406            if (axis_offset==0) {
1407              int s2=s[2];
1408              ev_shape.push_back(s2);
1409            }
1410            else if (axis_offset==1) {
1411              int s0=s[0];
1412              ev_shape.push_back(s0);
1413            }
1414            Data ev(0.,ev_shape,getFunctionSpace());
1415            ev.typeMatchRight(*this);
1416            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1417            return ev;
1418         }
1419         if (getDataPointRank()==4) {
1420            DataArrayView::ShapeType ev_shape;
1421            if (axis_offset==0) {
1422              ev_shape.push_back(s[2]);
1423              ev_shape.push_back(s[3]);
1424            }
1425            else if (axis_offset==1) {
1426              ev_shape.push_back(s[0]);
1427              ev_shape.push_back(s[3]);
1428            }
1429        else if (axis_offset==2) {
1430          ev_shape.push_back(s[0]);
1431          ev_shape.push_back(s[1]);
1432        }
1433            Data ev(0.,ev_shape,getFunctionSpace());
1434            ev.typeMatchRight(*this);
1435        m_data->matrixtrace(ev.m_data.get(), axis_offset);
1436            return ev;
1437         }
1438         else {
1439            throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1440         }
1441    }
1442    
1443    Data
1444    Data::transpose(int axis_offset) const
1445    {
1446    #if defined DOPROF
1447         profData->reduction2++;
1448    #endif
1449         DataArrayView::ShapeType s=getDataPointShape();
1450         DataArrayView::ShapeType ev_shape;
1451         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1452         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1453         int rank=getDataPointRank();
1454         if (axis_offset<0 || axis_offset>rank) {
1455            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1456         }
1457         for (int i=0; i<rank; i++) {
1458           int index = (axis_offset+i)%rank;
1459           ev_shape.push_back(s[index]); // Append to new shape
1460         }
1461         Data ev(0.,ev_shape,getFunctionSpace());
1462         ev.typeMatchRight(*this);
1463         m_data->transpose(ev.m_data.get(), axis_offset);
1464         return ev;
1465  }  }
1466    
1467  Data  Data
1468  Data::transpose(int axis) const  Data::eigenvalues() const
1469    {
1470         #if defined DOPROF
1471            profData->unary++;
1472         #endif
1473         // check input
1474         DataArrayView::ShapeType s=getDataPointShape();
1475         if (getDataPointRank()!=2)
1476            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1477         if(s[0] != s[1])
1478            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1479         // create return
1480         DataArrayView::ShapeType ev_shape(1,s[0]);
1481         Data ev(0.,ev_shape,getFunctionSpace());
1482         ev.typeMatchRight(*this);
1483         m_data->eigenvalues(ev.m_data.get());
1484         return ev;
1485    }
1486    
1487    const boost::python::tuple
1488    Data::eigenvalues_and_eigenvectors(const double tol) const
1489  {  {
1490  #if defined DOPROF       #if defined DOPROF
1491    profData->reduction2++;          profData->unary++;
1492  #endif       #endif
1493    // not implemented       DataArrayView::ShapeType s=getDataPointShape();
1494    throw DataException("Error - Data::transpose not implemented yet.");       if (getDataPointRank()!=2)
1495    return Data();          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1496         if(s[0] != s[1])
1497            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1498         // create return
1499         DataArrayView::ShapeType ev_shape(1,s[0]);
1500         Data ev(0.,ev_shape,getFunctionSpace());
1501         ev.typeMatchRight(*this);
1502         DataArrayView::ShapeType V_shape(2,s[0]);
1503         Data V(0.,V_shape,getFunctionSpace());
1504         V.typeMatchRight(*this);
1505         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1506         return make_tuple(boost::python::object(ev),boost::python::object(V));
1507  }  }
1508    
1509  const boost::python::tuple  const boost::python::tuple
# Line 1213  Data::mindp() const Line 1515  Data::mindp() const
1515    
1516    int SampleNo;    int SampleNo;
1517    int DataPointNo;    int DataPointNo;
1518        int ProcNo;
1519    calc_mindp(SampleNo,DataPointNo);    calc_mindp(ProcNo,SampleNo,DataPointNo);
1520      return make_tuple(ProcNo,SampleNo,DataPointNo);
   return make_tuple(SampleNo,DataPointNo);  
1521  }  }
1522    
1523  void  void
1524  Data::calc_mindp(int& SampleNo,  Data::calc_mindp(   int& ProcNo,
1525                   int& DataPointNo) const                  int& SampleNo,
1526            int& DataPointNo) const
1527  {  {
1528    int i,j;    int i,j;
1529    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1257  Data::calc_mindp(int& SampleNo, Line 1559  Data::calc_mindp(int& SampleNo,
1559      }      }
1560    }    }
1561    
1562    #ifdef PASO_MPI
1563        // determine the processor on which the minimum occurs
1564        next = temp.getDataPoint(lowi,lowj)();
1565        int lowProc = 0;
1566        double *globalMins = new double[get_MPISize()+1];
1567        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1568        
1569        if( get_MPIRank()==0 ){
1570            next = globalMins[lowProc];
1571            for( i=1; i<get_MPISize(); i++ )
1572                if( next>globalMins[i] ){
1573                    lowProc = i;
1574                    next = globalMins[i];
1575                }
1576        }
1577        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1578    
1579        delete [] globalMins;
1580        ProcNo = lowProc;
1581    #else
1582        ProcNo = 0;
1583    #endif
1584    SampleNo = lowi;    SampleNo = lowi;
1585    DataPointNo = lowj;    DataPointNo = lowj;
1586  }  }
# Line 1282  Data::saveVTK(std::string fileName) cons Line 1606  Data::saveVTK(std::string fileName) cons
1606  Data&  Data&
1607  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1608  {  {
1609      if (isProtected()) {
1610            throw DataException("Error - attempt to update protected Data object.");
1611      }
1612  #if defined DOPROF  #if defined DOPROF
1613    profData->binary++;    profData->binary++;
1614  #endif  #endif
# Line 1292  Data::operator+=(const Data& right) Line 1619  Data::operator+=(const Data& right)
1619  Data&  Data&
1620  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1621  {  {
1622      if (isProtected()) {
1623            throw DataException("Error - attempt to update protected Data object.");
1624      }
1625  #if defined DOPROF  #if defined DOPROF
1626    profData->binary++;    profData->binary++;
1627  #endif  #endif
# Line 1302  Data::operator+=(const boost::python::ob Line 1632  Data::operator+=(const boost::python::ob
1632  Data&  Data&
1633  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1634  {  {
1635      if (isProtected()) {
1636            throw DataException("Error - attempt to update protected Data object.");
1637      }
1638  #if defined DOPROF  #if defined DOPROF
1639    profData->binary++;    profData->binary++;
1640  #endif  #endif
# Line 1312  Data::operator-=(const Data& right) Line 1645  Data::operator-=(const Data& right)
1645  Data&  Data&
1646  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1647  {  {
1648      if (isProtected()) {
1649            throw DataException("Error - attempt to update protected Data object.");
1650      }
1651  #if defined DOPROF  #if defined DOPROF
1652    profData->binary++;    profData->binary++;
1653  #endif  #endif
# Line 1322  Data::operator-=(const boost::python::ob Line 1658  Data::operator-=(const boost::python::ob
1658  Data&  Data&
1659  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1660  {  {
1661      if (isProtected()) {
1662            throw DataException("Error - attempt to update protected Data object.");
1663      }
1664  #if defined DOPROF  #if defined DOPROF
1665    profData->binary++;    profData->binary++;
1666  #endif  #endif
# Line 1332  Data::operator*=(const Data& right) Line 1671  Data::operator*=(const Data& right)
1671  Data&  Data&
1672  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1673  {  {
1674      if (isProtected()) {
1675            throw DataException("Error - attempt to update protected Data object.");
1676      }
1677  #if defined DOPROF  #if defined DOPROF
1678    profData->binary++;    profData->binary++;
1679  #endif  #endif
# Line 1342  Data::operator*=(const boost::python::ob Line 1684  Data::operator*=(const boost::python::ob
1684  Data&  Data&
1685  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1686  {  {
1687      if (isProtected()) {
1688            throw DataException("Error - attempt to update protected Data object.");
1689      }
1690  #if defined DOPROF  #if defined DOPROF
1691    profData->binary++;    profData->binary++;
1692  #endif  #endif
# Line 1352  Data::operator/=(const Data& right) Line 1697  Data::operator/=(const Data& right)
1697  Data&  Data&
1698  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1699  {  {
1700      if (isProtected()) {
1701            throw DataException("Error - attempt to update protected Data object.");
1702      }
1703  #if defined DOPROF  #if defined DOPROF
1704    profData->binary++;    profData->binary++;
1705  #endif  #endif
# Line 1360  Data::operator/=(const boost::python::ob Line 1708  Data::operator/=(const boost::python::ob
1708  }  }
1709    
1710  Data  Data
1711    Data::rpowO(const boost::python::object& left) const
1712    {
1713      if (isProtected()) {
1714            throw DataException("Error - attempt to update protected Data object.");
1715      }
1716    #if defined DOPROF
1717      profData->binary++;
1718    #endif
1719      Data left_d(left,*this);
1720      return left_d.powD(*this);
1721    }
1722    
1723    Data
1724  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1725  {  {
1726  #if defined DOPROF  #if defined DOPROF
# Line 1383  Data::powD(const Data& right) const Line 1744  Data::powD(const Data& right) const
1744    return result;    return result;
1745  }  }
1746    
1747    
1748  //  //
1749  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1750  Data  Data
# Line 1596  escript::operator/(const boost::python:: Line 1958  escript::operator/(const boost::python::
1958  //  return ret;  //  return ret;
1959  //}  //}
1960    
1961    /* TODO */
1962    /* global reduction */
1963  Data  Data
1964  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1965  {  {
# Line 1610  Data::getItem(const boost::python::objec Line 1974  Data::getItem(const boost::python::objec
1974    return getSlice(slice_region);    return getSlice(slice_region);
1975  }  }
1976    
1977    /* TODO */
1978    /* global reduction */
1979  Data  Data
1980  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1981  {  {
# Line 1619  Data::getSlice(const DataArrayView::Regi Line 1985  Data::getSlice(const DataArrayView::Regi
1985    return Data(*this,region);    return Data(*this,region);
1986  }  }
1987    
1988    /* TODO */
1989    /* global reduction */
1990  void  void
1991  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1992                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1627  Data::setItemO(const boost::python::obje Line 1995  Data::setItemO(const boost::python::obje
1995    setItemD(key,tempData);    setItemD(key,tempData);
1996  }  }
1997    
1998    /* TODO */
1999    /* global reduction */
2000  void  void
2001  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2002                 const Data& value)                 const Data& value)
# Line 1644  Data::setItemD(const boost::python::obje Line 2014  Data::setItemD(const boost::python::obje
2014    }    }
2015  }  }
2016    
2017    /* TODO */
2018    /* global reduction */
2019  void  void
2020  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2021                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
2022  {  {
2023      if (isProtected()) {
2024            throw DataException("Error - attempt to update protected Data object.");
2025      }
2026  #if defined DOPROF  #if defined DOPROF
2027    profData->slicing++;    profData->slicing++;
2028  #endif  #endif
# Line 1685  Data::typeMatchRight(const Data& right) Line 2060  Data::typeMatchRight(const Data& right)
2060    }    }
2061  }  }
2062    
2063    /* TODO */
2064    /* global reduction */
2065  void  void
2066  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2067                       const boost::python::object& value)                       const boost::python::object& value)
2068  {  {
2069      if (isProtected()) {
2070            throw DataException("Error - attempt to update protected Data object.");
2071      }
2072    //    //
2073    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2074    tag();    tag();
# Line 1706  Data::setTaggedValue(int tagKey, Line 2086  Data::setTaggedValue(int tagKey,
2086    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2087  }  }
2088    
2089    /* TODO */
2090    /* global reduction */
2091  void  void
2092  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2093                              const DataArrayView& value)                              const DataArrayView& value)
2094  {  {
2095      if (isProtected()) {
2096            throw DataException("Error - attempt to update protected Data object.");
2097      }
2098    //    //
2099    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2100    tag();    tag();
# Line 1723  Data::setTaggedValueFromCPP(int tagKey, Line 2108  Data::setTaggedValueFromCPP(int tagKey,
2108    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
2109  }  }
2110    
2111    /* TODO */
2112    /* global reduction */
2113  int  int
2114  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2115  {  {
2116    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
2117  }  }
2118    
2119    /* TODO */
2120    /* global reduction */
2121  void  void
2122  Data::setRefValue(int ref,  Data::setRefValue(int ref,
2123                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
2124  {  {
2125      if (isProtected()) {
2126            throw DataException("Error - attempt to update protected Data object.");
2127      }
2128    //    //
2129    // Construct DataArray from boost::python::object input value    // Construct DataArray from boost::python::object input value
2130    DataArray valueDataArray(value);    DataArray valueDataArray(value);
# Line 1742  Data::setRefValue(int ref, Line 2134  Data::setRefValue(int ref,
2134    m_data->setRefValue(ref,valueDataArray);    m_data->setRefValue(ref,valueDataArray);
2135  }  }
2136    
2137    /* TODO */
2138    /* global reduction */
2139  void  void
2140  Data::getRefValue(int ref,  Data::getRefValue(int ref,
2141                    boost::python::numeric::array& value)                    boost::python::numeric::array& value)
# Line 1845  Data::archiveData(const std::string file Line 2239  Data::archiveData(const std::string file
2239    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2240    int dataLength = getLength();    int dataLength = getLength();
2241    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2242    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2243    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2245    }    }
2246    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2247    if (isTagged()) {    if (isTagged()) {
2248      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 2003  Data::extractData(const std::string file Line 2397  Data::extractData(const std::string file
2397        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2398      }      }
2399    }    }
2400    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2401    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2402      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2403    }    }
2404    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2405    if (dataType==2) {    if (dataType==2) {
2406      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2407        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 2146  ostream& escript::operator<<(ostream& o, Line 2540  ostream& escript::operator<<(ostream& o,
2540    o << data.toString();    o << data.toString();
2541    return o;    return o;
2542  }  }
2543    
2544    /* Member functions specific to the MPI implementation */
2545    
2546    void
2547    Data::print()
2548    {
2549      int i,j;
2550      
2551      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2552      for( i=0; i<getNumSamples(); i++ )
2553      {
2554        printf( "[%6d]", i );
2555        for( j=0; j<getNumDataPointsPerSample(); j++ )
2556          printf( "\t%10.7g", (getSampleData(i))[j] );
2557        printf( "\n" );
2558      }
2559    }
2560    
2561    int
2562    Data::get_MPISize() const
2563    {
2564        int error, size;
2565    #ifdef PASO_MPI
2566        error = MPI_Comm_size( get_MPIComm(), &size );
2567    #else
2568        size = 1;
2569    #endif
2570        return size;
2571    }
2572    
2573    int
2574    Data::get_MPIRank() const
2575    {
2576        int error, rank;
2577    #ifdef PASO_MPI
2578        error = MPI_Comm_rank( get_MPIComm(), &rank );
2579    #else
2580        rank = 0;
2581    #endif
2582        return rank;
2583    }
2584    
2585    MPI_Comm
2586    Data::get_MPIComm() const
2587    {
2588    #ifdef PASO_MPI
2589        return MPI_COMM_WORLD;
2590    #else
2591        return -1;
2592    #endif
2593    }
2594    

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

  ViewVC Help
Powered by ViewVC 1.1.26