/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 474 by jgs, Mon Jan 30 04:23:44 2006 UTC revision 804 by gross, Thu Aug 10 01:12:16 2006 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
   
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *  
  *                                                                            *  
  * This software is the property of ACcESS.  No part of this code             *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that                          *  
  * person has a software license agreement with ACcESS.                       *  
  *                                                                            *  
  ******************************************************************************  
   
 ******************************************************************************/  
2    
3    /*
4     ************************************************************
5     *          Copyright 2006 by ACcESS MNRF                   *
6     *                                                          *
7     *              http://www.access.edu.au                    *
8     *       Primary Business: Queensland, Australia            *
9     *  Licensed under the Open Software License version 3.0    *
10     *     http://www.opensource.org/licenses/osl-3.0.php       *
11     *                                                          *
12     ************************************************************
13    */
14  #include "Data.h"  #include "Data.h"
15    
 #include <iostream>  
 #include <fstream>  
 #include <algorithm>  
 #include <vector>  
 #include <exception>  
 #include <functional>  
 #include <math.h>  
   
 #include <boost/python/dict.hpp>  
 #include <boost/python/str.hpp>  
 #include <boost/python/extract.hpp>  
 #include <boost/python/long.hpp>  
 #include <boost/python/tuple.hpp>  
   
 #include "DataException.h"  
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataConstant.h"  #include "DataConstant.h"
18  #include "DataTagged.h"  #include "DataTagged.h"
19  #include "DataEmpty.h"  #include "DataEmpty.h"
20  #include "DataArray.h"  #include "DataArray.h"
21    #include "DataArrayView.h"
22  #include "DataProf.h"  #include "DataProf.h"
23  #include "FunctionSpaceFactory.h"  #include "FunctionSpaceFactory.h"
24  #include "AbstractContinuousDomain.h"  #include "AbstractContinuousDomain.h"
25  #include "UnaryFuncs.h"  #include "UnaryFuncs.h"
26    
27    #include <fstream>
28    #include <algorithm>
29    #include <vector>
30    #include <functional>
31    
32    #include <boost/python/dict.hpp>
33    #include <boost/python/extract.hpp>
34    #include <boost/python/long.hpp>
35    
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
# Line 60  Data::Data() Line 51  Data::Data()
51    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
52    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
53    m_data=temp_data;    m_data=temp_data;
54      m_protected=false;
55  #if defined DOPROF  #if defined DOPROF
56    // create entry in global profiling table for this object    // create entry in global profiling table for this object
57    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 77  Data::Data(double value, Line 69  Data::Data(double value,
69    }    }
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
72      m_protected=false;
73  #if defined DOPROF  #if defined DOPROF
74    // create entry in global profiling table for this object    // create entry in global profiling table for this object
75    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 91  Data::Data(double value, Line 84  Data::Data(double value,
84    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
85    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
86    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
87      m_protected=false;
88  #if defined DOPROF  #if defined DOPROF
89    // create entry in global profiling table for this object    // create entry in global profiling table for this object
90    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 100  Data::Data(double value, Line 94  Data::Data(double value,
94  Data::Data(const Data& inData)  Data::Data(const Data& inData)
95  {  {
96    m_data=inData.m_data;    m_data=inData.m_data;
97      m_protected=inData.isProtected();
98  #if defined DOPROF  #if defined DOPROF
99    // create entry in global profiling table for this object    // create entry in global profiling table for this object
100    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 114  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
110    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
111    m_data=temp_data;    m_data=temp_data;
112      m_protected=false;
113  #if defined DOPROF  #if defined DOPROF
114    // create entry in global profiling table for this object    // create entry in global profiling table for this object
115    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 123  Data::Data(const Data& inData, Line 119  Data::Data(const Data& inData,
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122    #if defined DOPROF
123      // create entry in global profiling table for this object
124      profData = dataProfTable.newData();
125    #endif
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    } else {
129        #if defined DOPROF
130        profData->interpolate++;
131        #endif
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
133      // Note: Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
134      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
# Line 139  Data::Data(const Data& inData, Line 142  Data::Data(const Data& inData,
142      }      }
143      m_data=tmp.m_data;      m_data=tmp.m_data;
144    }    }
145  #if defined DOPROF    m_protected=false;
   // create entry in global profiling table for this object  
   profData = dataProfTable.newData();  
 #endif  
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 154  Data::Data(const DataTagged::TagListType Line 154  Data::Data(const DataTagged::TagListType
154    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
155    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
156    m_data=temp_data;    m_data=temp_data;
157      m_protected=false;
158    if (expanded) {    if (expanded) {
159      expand();      expand();
160    }    }
# Line 168  Data::Data(const numeric::array& value, Line 169  Data::Data(const numeric::array& value,
169             bool expanded)             bool expanded)
170  {  {
171    initialise(value,what,expanded);    initialise(value,what,expanded);
172      m_protected=false;
173  #if defined DOPROF  #if defined DOPROF
174    // create entry in global profiling table for this object    // create entry in global profiling table for this object
175    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 179  Data::Data(const DataArrayView& value, Line 181  Data::Data(const DataArrayView& value,
181             bool expanded)             bool expanded)
182  {  {
183    initialise(value,what,expanded);    initialise(value,what,expanded);
184      m_protected=false;
185  #if defined DOPROF  #if defined DOPROF
186    // create entry in global profiling table for this object    // create entry in global profiling table for this object
187    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 191  Data::Data(const object& value, Line 194  Data::Data(const object& value,
194  {  {
195    numeric::array asNumArray(value);    numeric::array asNumArray(value);
196    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
197      m_protected=false;
198  #if defined DOPROF  #if defined DOPROF
199    // create entry in global profiling table for this object    // create entry in global profiling table for this object
200    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 215  Data::Data(const object& value, Line 219  Data::Data(const object& value,
219      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
220      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
221    }    }
222      m_protected=false;
223  #if defined DOPROF  #if defined DOPROF
224    // create entry in global profiling table for this object    // create entry in global profiling table for this object
225    profData = dataProfTable.newData();    profData = dataProfTable.newData();
# Line 345  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 360  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 444  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 527  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 756  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 804  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 842  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 1105  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 1175  Data::minval() const Line 1316  Data::minval() const
1316  }  }
1317    
1318  Data  Data
1319  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1320  {  {
1321  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1322    profData->reduction2++;       #if defined DOPROF
1323  #endif       profData->unary++;
1324    Trace trace_func;       #endif
1325    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1326         DataArrayView::ShapeType ev_shape;
1327         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1328         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1329         int rank=getDataPointRank();
1330         if (rank<2) {
1331            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1332         }
1333         if (axis0<0 || axis0>rank-1) {
1334            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1335         }
1336         if (axis1<0 || axis1>rank-1) {
1337             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1338         }
1339         if (axis0 == axis1) {
1340             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1341         }
1342         if (axis0 > axis1) {
1343             axis0_tmp=axis1;
1344             axis1_tmp=axis0;
1345         } else {
1346             axis0_tmp=axis0;
1347             axis1_tmp=axis1;
1348         }
1349         for (int i=0; i<rank; i++) {
1350           if (i == axis0_tmp) {
1351              ev_shape.push_back(s[axis1_tmp]);
1352           } else if (i == axis1_tmp) {
1353              ev_shape.push_back(s[axis0_tmp]);
1354           } else {
1355              ev_shape.push_back(s[i]);
1356           }
1357         }
1358         Data ev(0.,ev_shape,getFunctionSpace());
1359         ev.typeMatchRight(*this);
1360         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1361         return ev;
1362    
1363    }
1364    
1365    Data
1366    Data::symmetric() const
1367    {
1368         #if defined DOPROF
1369            profData->unary++;
1370         #endif
1371         // check input
1372         DataArrayView::ShapeType s=getDataPointShape();
1373         if (getDataPointRank()==2) {
1374            if(s[0] != s[1])
1375               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1376         }
1377         else if (getDataPointRank()==4) {
1378            if(!(s[0] == s[2] && s[1] == s[3]))
1379               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1380         }
1381         else {
1382            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1383         }
1384         Data ev(0.,getDataPointShape(),getFunctionSpace());
1385         ev.typeMatchRight(*this);
1386         m_data->symmetric(ev.m_data.get());
1387         return ev;
1388    }
1389    
1390    Data
1391    Data::nonsymmetric() const
1392    {
1393         #if defined DOPROF
1394            profData->unary++;
1395         #endif
1396         // check input
1397         DataArrayView::ShapeType s=getDataPointShape();
1398         if (getDataPointRank()==2) {
1399            if(s[0] != s[1])
1400               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1401            DataArrayView::ShapeType ev_shape;
1402            ev_shape.push_back(s[0]);
1403            ev_shape.push_back(s[1]);
1404            Data ev(0.,ev_shape,getFunctionSpace());
1405            ev.typeMatchRight(*this);
1406            m_data->nonsymmetric(ev.m_data.get());
1407            return ev;
1408         }
1409         else if (getDataPointRank()==4) {
1410            if(!(s[0] == s[2] && s[1] == s[3]))
1411               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1412            DataArrayView::ShapeType ev_shape;
1413            ev_shape.push_back(s[0]);
1414            ev_shape.push_back(s[1]);
1415            ev_shape.push_back(s[2]);
1416            ev_shape.push_back(s[3]);
1417            Data ev(0.,ev_shape,getFunctionSpace());
1418            ev.typeMatchRight(*this);
1419            m_data->nonsymmetric(ev.m_data.get());
1420            return ev;
1421         }
1422         else {
1423            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1424         }
1425    }
1426    
1427    Data
1428    Data::trace(int axis_offset) const
1429    {
1430         #if defined DOPROF
1431            profData->unary++;
1432         #endif
1433         DataArrayView::ShapeType s=getDataPointShape();
1434         if (getDataPointRank()==2) {
1435            DataArrayView::ShapeType ev_shape;
1436            Data ev(0.,ev_shape,getFunctionSpace());
1437            ev.typeMatchRight(*this);
1438            m_data->trace(ev.m_data.get(), axis_offset);
1439            return ev;
1440         }
1441         if (getDataPointRank()==3) {
1442            DataArrayView::ShapeType ev_shape;
1443            if (axis_offset==0) {
1444              int s2=s[2];
1445              ev_shape.push_back(s2);
1446            }
1447            else if (axis_offset==1) {
1448              int s0=s[0];
1449              ev_shape.push_back(s0);
1450            }
1451            Data ev(0.,ev_shape,getFunctionSpace());
1452            ev.typeMatchRight(*this);
1453            m_data->trace(ev.m_data.get(), axis_offset);
1454            return ev;
1455         }
1456         if (getDataPointRank()==4) {
1457            DataArrayView::ShapeType ev_shape;
1458            if (axis_offset==0) {
1459              ev_shape.push_back(s[2]);
1460              ev_shape.push_back(s[3]);
1461            }
1462            else if (axis_offset==1) {
1463              ev_shape.push_back(s[0]);
1464              ev_shape.push_back(s[3]);
1465            }
1466        else if (axis_offset==2) {
1467          ev_shape.push_back(s[0]);
1468          ev_shape.push_back(s[1]);
1469        }
1470            Data ev(0.,ev_shape,getFunctionSpace());
1471            ev.typeMatchRight(*this);
1472        m_data->trace(ev.m_data.get(), axis_offset);
1473            return ev;
1474         }
1475         else {
1476            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1477         }
1478  }  }
1479    
1480  Data  Data
1481  Data::transpose(int axis) const  Data::transpose(int axis_offset) const
1482    {
1483         #if defined DOPROF
1484         profData->unary++;
1485         #endif
1486         DataArrayView::ShapeType s=getDataPointShape();
1487         DataArrayView::ShapeType ev_shape;
1488         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1489         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1490         int rank=getDataPointRank();
1491         if (axis_offset<0 || axis_offset>rank) {
1492            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1493         }
1494         for (int i=0; i<rank; i++) {
1495           int index = (axis_offset+i)%rank;
1496           ev_shape.push_back(s[index]); // Append to new shape
1497         }
1498         Data ev(0.,ev_shape,getFunctionSpace());
1499         ev.typeMatchRight(*this);
1500         m_data->transpose(ev.m_data.get(), axis_offset);
1501         return ev;
1502    }
1503    
1504    Data
1505    Data::eigenvalues() const
1506    {
1507         #if defined DOPROF
1508            profData->unary++;
1509         #endif
1510         // check input
1511         DataArrayView::ShapeType s=getDataPointShape();
1512         if (getDataPointRank()!=2)
1513            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1514         if(s[0] != s[1])
1515            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1516         // create return
1517         DataArrayView::ShapeType ev_shape(1,s[0]);
1518         Data ev(0.,ev_shape,getFunctionSpace());
1519         ev.typeMatchRight(*this);
1520         m_data->eigenvalues(ev.m_data.get());
1521         return ev;
1522    }
1523    
1524    const boost::python::tuple
1525    Data::eigenvalues_and_eigenvectors(const double tol) const
1526  {  {
1527  #if defined DOPROF       #if defined DOPROF
1528    profData->reduction2++;          profData->unary++;
1529  #endif       #endif
1530    // not implemented       DataArrayView::ShapeType s=getDataPointShape();
1531    throw DataException("Error - Data::transpose not implemented yet.");       if (getDataPointRank()!=2)
1532    return Data();          throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1533         if(s[0] != s[1])
1534            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1535         // create return
1536         DataArrayView::ShapeType ev_shape(1,s[0]);
1537         Data ev(0.,ev_shape,getFunctionSpace());
1538         ev.typeMatchRight(*this);
1539         DataArrayView::ShapeType V_shape(2,s[0]);
1540         Data V(0.,V_shape,getFunctionSpace());
1541         V.typeMatchRight(*this);
1542         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1543         return make_tuple(boost::python::object(ev),boost::python::object(V));
1544  }  }
1545    
1546  const boost::python::tuple  const boost::python::tuple
# Line 1204  Data::mindp() const Line 1552  Data::mindp() const
1552    
1553    int SampleNo;    int SampleNo;
1554    int DataPointNo;    int DataPointNo;
1555        int ProcNo;
1556    calc_mindp(SampleNo,DataPointNo);    calc_mindp(ProcNo,SampleNo,DataPointNo);
1557      return make_tuple(ProcNo,SampleNo,DataPointNo);
   return make_tuple(SampleNo,DataPointNo);  
1558  }  }
1559    
1560  void  void
1561  Data::calc_mindp(int& SampleNo,  Data::calc_mindp(   int& ProcNo,
1562                   int& DataPointNo) const                  int& SampleNo,
1563            int& DataPointNo) const
1564  {  {
1565    int i,j;    int i,j;
1566    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1248  Data::calc_mindp(int& SampleNo, Line 1596  Data::calc_mindp(int& SampleNo,
1596      }      }
1597    }    }
1598    
1599    #ifdef PASO_MPI
1600        // determine the processor on which the minimum occurs
1601        next = temp.getDataPoint(lowi,lowj)();
1602        int lowProc = 0;
1603        double *globalMins = new double[get_MPISize()+1];
1604        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1605        
1606        if( get_MPIRank()==0 ){
1607            next = globalMins[lowProc];
1608            for( i=1; i<get_MPISize(); i++ )
1609                if( next>globalMins[i] ){
1610                    lowProc = i;
1611                    next = globalMins[i];
1612                }
1613        }
1614        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1615    
1616        delete [] globalMins;
1617        ProcNo = lowProc;
1618    #else
1619        ProcNo = 0;
1620    #endif
1621    SampleNo = lowi;    SampleNo = lowi;
1622    DataPointNo = lowj;    DataPointNo = lowj;
1623  }  }
# Line 1273  Data::saveVTK(std::string fileName) cons Line 1643  Data::saveVTK(std::string fileName) cons
1643  Data&  Data&
1644  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1645  {  {
1646      if (isProtected()) {
1647            throw DataException("Error - attempt to update protected Data object.");
1648      }
1649  #if defined DOPROF  #if defined DOPROF
1650    profData->binary++;    profData->binary++;
1651  #endif  #endif
# Line 1283  Data::operator+=(const Data& right) Line 1656  Data::operator+=(const Data& right)
1656  Data&  Data&
1657  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1658  {  {
1659      if (isProtected()) {
1660            throw DataException("Error - attempt to update protected Data object.");
1661      }
1662  #if defined DOPROF  #if defined DOPROF
1663    profData->binary++;    profData->binary++;
1664  #endif  #endif
# Line 1293  Data::operator+=(const boost::python::ob Line 1669  Data::operator+=(const boost::python::ob
1669  Data&  Data&
1670  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1671  {  {
1672      if (isProtected()) {
1673            throw DataException("Error - attempt to update protected Data object.");
1674      }
1675  #if defined DOPROF  #if defined DOPROF
1676    profData->binary++;    profData->binary++;
1677  #endif  #endif
# Line 1303  Data::operator-=(const Data& right) Line 1682  Data::operator-=(const Data& right)
1682  Data&  Data&
1683  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1684  {  {
1685      if (isProtected()) {
1686            throw DataException("Error - attempt to update protected Data object.");
1687      }
1688  #if defined DOPROF  #if defined DOPROF
1689    profData->binary++;    profData->binary++;
1690  #endif  #endif
# Line 1313  Data::operator-=(const boost::python::ob Line 1695  Data::operator-=(const boost::python::ob
1695  Data&  Data&
1696  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1697  {  {
1698      if (isProtected()) {
1699            throw DataException("Error - attempt to update protected Data object.");
1700      }
1701  #if defined DOPROF  #if defined DOPROF
1702    profData->binary++;    profData->binary++;
1703  #endif  #endif
# Line 1323  Data::operator*=(const Data& right) Line 1708  Data::operator*=(const Data& right)
1708  Data&  Data&
1709  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1710  {  {
1711      if (isProtected()) {
1712            throw DataException("Error - attempt to update protected Data object.");
1713      }
1714  #if defined DOPROF  #if defined DOPROF
1715    profData->binary++;    profData->binary++;
1716  #endif  #endif
# Line 1333  Data::operator*=(const boost::python::ob Line 1721  Data::operator*=(const boost::python::ob
1721  Data&  Data&
1722  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1723  {  {
1724      if (isProtected()) {
1725            throw DataException("Error - attempt to update protected Data object.");
1726      }
1727  #if defined DOPROF  #if defined DOPROF
1728    profData->binary++;    profData->binary++;
1729  #endif  #endif
# Line 1343  Data::operator/=(const Data& right) Line 1734  Data::operator/=(const Data& right)
1734  Data&  Data&
1735  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1736  {  {
1737      if (isProtected()) {
1738            throw DataException("Error - attempt to update protected Data object.");
1739      }
1740  #if defined DOPROF  #if defined DOPROF
1741    profData->binary++;    profData->binary++;
1742  #endif  #endif
# Line 1351  Data::operator/=(const boost::python::ob Line 1745  Data::operator/=(const boost::python::ob
1745  }  }
1746    
1747  Data  Data
1748    Data::rpowO(const boost::python::object& left) const
1749    {
1750      if (isProtected()) {
1751            throw DataException("Error - attempt to update protected Data object.");
1752      }
1753    #if defined DOPROF
1754      profData->binary++;
1755    #endif
1756      Data left_d(left,*this);
1757      return left_d.powD(*this);
1758    }
1759    
1760    Data
1761  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1762  {  {
1763  #if defined DOPROF  #if defined DOPROF
# Line 1374  Data::powD(const Data& right) const Line 1781  Data::powD(const Data& right) const
1781    return result;    return result;
1782  }  }
1783    
1784    
1785  //  //
1786  // NOTE: It is essential to specify the namespace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1787  Data  Data
# Line 1587  escript::operator/(const boost::python:: Line 1995  escript::operator/(const boost::python::
1995  //  return ret;  //  return ret;
1996  //}  //}
1997    
1998    /* TODO */
1999    /* global reduction */
2000  Data  Data
2001  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
2002  {  {
# Line 1601  Data::getItem(const boost::python::objec Line 2011  Data::getItem(const boost::python::objec
2011    return getSlice(slice_region);    return getSlice(slice_region);
2012  }  }
2013    
2014    /* TODO */
2015    /* global reduction */
2016  Data  Data
2017  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
2018  {  {
# Line 1610  Data::getSlice(const DataArrayView::Regi Line 2022  Data::getSlice(const DataArrayView::Regi
2022    return Data(*this,region);    return Data(*this,region);
2023  }  }
2024    
2025    /* TODO */
2026    /* global reduction */
2027  void  void
2028  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
2029                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1618  Data::setItemO(const boost::python::obje Line 2032  Data::setItemO(const boost::python::obje
2032    setItemD(key,tempData);    setItemD(key,tempData);
2033  }  }
2034    
2035    /* TODO */
2036    /* global reduction */
2037  void  void
2038  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2039                 const Data& value)                 const Data& value)
# Line 1635  Data::setItemD(const boost::python::obje Line 2051  Data::setItemD(const boost::python::obje
2051    }    }
2052  }  }
2053    
2054    /* TODO */
2055    /* global reduction */
2056  void  void
2057  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2058                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
2059  {  {
2060      if (isProtected()) {
2061            throw DataException("Error - attempt to update protected Data object.");
2062      }
2063  #if defined DOPROF  #if defined DOPROF
2064    profData->slicing++;    profData->slicing++;
2065  #endif  #endif
# Line 1676  Data::typeMatchRight(const Data& right) Line 2097  Data::typeMatchRight(const Data& right)
2097    }    }
2098  }  }
2099    
2100    /* TODO */
2101    /* global reduction */
2102  void  void
2103  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2104                       const boost::python::object& value)                       const boost::python::object& value)
2105  {  {
2106      if (isProtected()) {
2107            throw DataException("Error - attempt to update protected Data object.");
2108      }
2109    //    //
2110    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2111    tag();    tag();
# Line 1697  Data::setTaggedValue(int tagKey, Line 2123  Data::setTaggedValue(int tagKey,
2123    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2124  }  }
2125    
2126    /* TODO */
2127    /* global reduction */
2128  void  void
2129  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2130                              const DataArrayView& value)                              const DataArrayView& value)
2131  {  {
2132      if (isProtected()) {
2133            throw DataException("Error - attempt to update protected Data object.");
2134      }
2135    //    //
2136    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2137    tag();    tag();
# Line 1714  Data::setTaggedValueFromCPP(int tagKey, Line 2145  Data::setTaggedValueFromCPP(int tagKey,
2145    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
2146  }  }
2147    
2148    /* TODO */
2149    /* global reduction */
2150  int  int
2151  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2152  {  {
2153    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
2154  }  }
2155    
2156    /* TODO */
2157    /* global reduction */
2158  void  void
2159  Data::setRefValue(int ref,  Data::setRefValue(int ref,
2160                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
2161  {  {
2162      if (isProtected()) {
2163            throw DataException("Error - attempt to update protected Data object.");
2164      }
2165    //    //
2166    // Construct DataArray from boost::python::object input value    // Construct DataArray from boost::python::object input value
2167    DataArray valueDataArray(value);    DataArray valueDataArray(value);
# Line 1733  Data::setRefValue(int ref, Line 2171  Data::setRefValue(int ref,
2171    m_data->setRefValue(ref,valueDataArray);    m_data->setRefValue(ref,valueDataArray);
2172  }  }
2173    
2174    /* TODO */
2175    /* global reduction */
2176  void  void
2177  Data::getRefValue(int ref,  Data::getRefValue(int ref,
2178                    boost::python::numeric::array& value)                    boost::python::numeric::array& value)
# Line 1836  Data::archiveData(const std::string file Line 2276  Data::archiveData(const std::string file
2276    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2277    int dataLength = getLength();    int dataLength = getLength();
2278    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2279    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2280    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2281      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2282    }    }
2283    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2284    if (isTagged()) {    if (isTagged()) {
2285      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2286        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 1994  Data::extractData(const std::string file Line 2434  Data::extractData(const std::string file
2434        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2435      }      }
2436    }    }
2437    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2438    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2439      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2440    }    }
2441    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2442    if (dataType==2) {    if (dataType==2) {
2443      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2444        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 2137  ostream& escript::operator<<(ostream& o, Line 2577  ostream& escript::operator<<(ostream& o,
2577    o << data.toString();    o << data.toString();
2578    return o;    return o;
2579  }  }
2580    
2581    /* Member functions specific to the MPI implementation */
2582    
2583    void
2584    Data::print()
2585    {
2586      int i,j;
2587      
2588      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2589      for( i=0; i<getNumSamples(); i++ )
2590      {
2591        printf( "[%6d]", i );
2592        for( j=0; j<getNumDataPointsPerSample(); j++ )
2593          printf( "\t%10.7g", (getSampleData(i))[j] );
2594        printf( "\n" );
2595      }
2596    }
2597    
2598    int
2599    Data::get_MPISize() const
2600    {
2601        int error, size;
2602    #ifdef PASO_MPI
2603        error = MPI_Comm_size( get_MPIComm(), &size );
2604    #else
2605        size = 1;
2606    #endif
2607        return size;
2608    }
2609    
2610    int
2611    Data::get_MPIRank() const
2612    {
2613        int error, rank;
2614    #ifdef PASO_MPI
2615        error = MPI_Comm_rank( get_MPIComm(), &rank );
2616    #else
2617        rank = 0;
2618    #endif
2619        return rank;
2620    }
2621    
2622    MPI_Comm
2623    Data::get_MPIComm() const
2624    {
2625    #ifdef PASO_MPI
2626        return MPI_COMM_WORLD;
2627    #else
2628        return -1;
2629    #endif
2630    }
2631    

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

  ViewVC Help
Powered by ViewVC 1.1.26