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

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

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

revision 584 by gross, Thu Mar 9 23:03:38 2006 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
 // $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.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2    /*****************************************************************************
3    *
4    * Copyright (c) 2003-2012 by University of Queensland
5    * http://www.uq.edu.au
6    *
7    * Primary Business: Queensland, Australia
8    * Licensed under the Open Software License version 3.0
9    * http://www.opensource.org/licenses/osl-3.0.php
10    *
11    * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    
17    #include "Data.h"
18  #include "DataExpanded.h"  #include "DataExpanded.h"
19  #include "DataException.h"  #include "DataException.h"
20  #include "DataConstant.h"  #include "DataConstant.h"
21  #include "DataTagged.h"  #include "DataTagged.h"
22    #include <limits>
23    
24    #include "esysUtils/Esys_MPI.h"
25    
26    #ifdef USE_NETCDF
27    #include <netcdfcpp.h>
28    #endif
29    
30  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
31    #include "DataMaths.h"
32    
33    //#define MKLRANDOM
34    
35    #ifdef MKLRANDOM
36    #include <mkl_vsl.h>
37    #else
38    #include <boost/random/mersenne_twister.hpp>
39    #endif
40    
41  using namespace std;  using namespace std;
42  using namespace boost::python;  using namespace boost::python;
43  using namespace boost;  using namespace boost;
44    using namespace escript::DataTypes;
45    
46    
47    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
48    
49    #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; *((int*)0)=17;throw DataException(ss.str());}
50    
51  namespace escript {  namespace escript {
52    
53  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const WrappedArray& value,
54                             const FunctionSpace& what)                             const FunctionSpace& what)
55    : DataAbstract(what)    : parent(what,value.getShape())
56  {  {
   DataArrayView::ShapeType tempShape;  
   //  
   // extract the shape of the python numarray  
   for (int i=0; i<value.getrank(); i++) {  
     tempShape.push_back(extract<int>(value.getshape()[i]));  
   }  
57    //    //
58    // initialise the data array for this object    // initialise the data array for this object
59    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
60    //    //
61    // copy the given value to every data point    // copy the given value to every data point
62    copy(value);    copy(value);
63  }  }
64    
65  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
66    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(), other.getShape()),
67    m_data(other.m_data)    m_data(other.m_data)
68  {  {
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
69  }  }
70    
71  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
72    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
73  {  {
74    //    //
75    // initialise the data array for this object    // initialise the data array for this object
76    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
77    //    //
78    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
79    copy(other.getPointDataView());    copy(other);
80  }  }
81    
82  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
83    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
84  {  {
85    //    //
86    // initialise the data array for this object    // initialise the data array for this object
87    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
88    //    //
89    // for each data point in this object, extract and copy the corresponding data    // for each data point in this object, extract and copy the corresponding data
90    // value from the given DataTagged object    // value from the given DataTagged object
91    int i,j;    int i,j;
92    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
93    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
94    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
95    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
96      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
97        try {        try {
98          getPointDataView().copy(getPointOffset(i,j),             DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
99                                  other.getPointDataView(),                                  other.getVectorRO(),
100                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
101        }        }
102        catch (std::exception& e) {        catch (std::exception& e) {
# Line 93  DataExpanded::DataExpanded(const DataTag Line 107  DataExpanded::DataExpanded(const DataTag
107  }  }
108    
109  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
110                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
111    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
112  {  {
113    //    //
114    // get the shape of the slice    // get the shape of the slice
115    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
116    //    //
117    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
118    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
119    //    //
120    // copy the data    // copy the data
121    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
122    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
123    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
124    int i,j;    int i,j;
125    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
126    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
127      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
128        try {        try {
129          getPointDataView().copySlice(getPointOffset(i,j),          DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
130                                       other.getPointDataView(),                                       other.getVectorRO(),
131                         other.getShape(),
132                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
133                                       region_loop_range);                                       region_loop_range);
134        }        }
# Line 124  DataExpanded::DataExpanded(const DataExp Line 139  DataExpanded::DataExpanded(const DataExp
139    }    }
140  }  }
141    
142  DataExpanded::DataExpanded(const DataArrayView& value,  DataExpanded::DataExpanded(const FunctionSpace& what,
143                             const FunctionSpace& what)                             const DataTypes::ShapeType &shape,
144    : DataAbstract(what)                             const DataTypes::ValueType &data)
145  {    : parent(what,shape)
146    //  {
147    // get the shape of the given data value    EsysAssert(data.size()%getNoValues()==0,
148    DataArrayView::ShapeType tempShape=value.getShape();                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
149    //  
150    // initialise this Data object to the shape of the given data value    if (data.size()==getNoValues())
151    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    {
152    //       ValueType& vec=m_data.getData();
153    // copy the given value to every data point       //
154    copy(value);       // create the view of the data
155         initialise(what.getNumSamples(),what.getNumDPPSample());
156         // now we copy this value to all elements
157         for (int i=0;i<getLength();)
158         {
159        for (unsigned int j=0;j<getNoValues();++j,++i)
160        {
161            vec[i]=data[j];
162        }
163         }
164      }
165      else
166      {
167         //
168         // copy the data in the correct format
169         m_data.getData()=data;
170      }
171    
172    
173  }  }
174    
175  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
176                             const DataArrayView::ShapeType &shape,                             const DataTypes::ShapeType &shape,
177                             const DataArrayView::ValueType &data)                             const double v)
178    : DataAbstract(what)    : parent(what,shape)
179  {  {
180    //       ValueType& vec=m_data.getData();
181    // create the view of the data       //
182    initialise(shape,what.getNumSamples(),what.getNumDPPSample());       // create the view of the data
183    //       initialise(what.getNumSamples(),what.getNumDPPSample());
184    // copy the data in the correct format       // now we copy this value to all elements
185    m_data.getData()=data;       const int L=getLength();
186         int i;
187         #pragma omp parallel for schedule(static) private(i)
188         for (i=0;i<L;++i)
189         {
190            vec[i]=v;
191         }
192  }  }
193    
194    
195    
196  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
197  {  {
198  }  }
199    
200  void  DataAbstract*
201  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataExpanded::deepCopy()
202  {  {
203    if (getPointDataView().getRank()!=0) {    return new DataExpanded(*this);
     stringstream temp;  
     temp << "Error - Can only reshape Data with data points of rank 0. "  
          << "This Data has data points with rank: "  
          << getPointDataView().getRank();  
     throw DataException(temp.str());  
   }  
   //  
   // create the new DataBlocks2D data array, and a corresponding DataArrayView  
   DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));  
   DataArrayView newView(newData.getData(),shape);  
   //  
   // Copy the original data to every value for the new shape  
   int i,j;  
   int nRows=m_data.getNumRows();  
   int nCols=m_data.getNumCols();  
   #pragma omp parallel for private(i,j) schedule(static)  
   for (i=0;i<nRows;i++) {  
     for (j=0;j<nCols;j++) {  
       // NOTE: An exception may be thown from this call if  
       // DOASSERT is on which of course will play  
       // havoc with the omp threads. Run single threaded  
       // if using DOASSERT.  
       newView.copy(newData.index(i,j),m_data(i,j));  
     }  
   }  
   // swap the new data array for the original  
   m_data.Swap(newData);  
   // set the corresponding DataArrayView  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
204  }  }
205    
206    
207  DataAbstract*  DataAbstract*
208  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataTypes::RegionType& region) const
209  {  {
210    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
211  }  }
212    
213  void  void
214  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
215                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
216  {  {
217    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
218    if (tempDataExp==0) {    if (tempDataExp==0) {
219      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
220    }    }
221      CHECK_FOR_EX_WRITE
222    //    //
223    // get shape of slice    // get shape of slice
224    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
225    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
226    //    //
227    // check shape    // check shape
228    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
229      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
230    }    }
231    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
232      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
233          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
234    }    }
235    //    //
236    // copy the data from the slice into this object    // copy the data from the slice into this object
237    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
238    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
239    int i, j;    int i, j;
240      ValueType& vec=getVectorRW();
241      const ShapeType& mshape=getShape();
242      const ValueType& tVec=tempDataExp->getVectorRO();
243      const ShapeType& tShape=tempDataExp->getShape();
244    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
245    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
246      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
247        getPointDataView().copySliceFrom(getPointOffset(i,j),          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
248                                         tempDataExp->getPointDataView(),                                         tVec,
249                           tShape,
250                                         tempDataExp->getPointOffset(i,j),                                         tempDataExp->getPointOffset(i,j),
251                                         region_loop_range);                                         region_loop_range);
252    
253      }      }
254    }    }
255  }  }
256    
257  void  void
258  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
259  {  {
260      EsysAssert((checkShape(getShape(), value.getShape())),
261                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
262    
263    //    //
264    // copy a single value to every data point in this object    // copy a single value to every data point in this object
265    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
# Line 246  DataExpanded::copy(const DataArrayView& Line 268  DataExpanded::copy(const DataArrayView&
268    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
269    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
270      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
271        // NOTE: An exception may be thown from this call if        DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
       // DOASSERT is on which of course will play  
       // havoc with the omp threads. Run single threaded  
       // if using DOASSERT.  
       getPointDataView().copy(m_data.index(i,j),value);  
272      }      }
273    }    }
274  }  }
275    
276  void  void
277  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const WrappedArray& value)
278  {  {
   //  
   // first convert the numarray into a DataArray object  
   DataArray temp(value);  
   //  
279    // check the input shape matches this shape    // check the input shape matches this shape
280    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!DataTypes::checkShape(getShape(),value.getShape())) {
281      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
282                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
283                          temp.getView().getShape()));                          value.getShape(),getShape()));
284    }    }
285    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the data  
   copy(temp.getView());  
286  }  }
287    
288    
289  void  void
290  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
291                           int noDataPointsPerSample)                           int noDataPointsPerSample)
292  {  {
293      if (noSamples==0)     //retain the default empty object
294      {
295         return;
296      }
297    //    //
298    // resize data array to the required size    // resize data array to the required size
299    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
300    //  }
301    // create the data view of the data array  
302    DataArrayView temp(m_data.getData(),shape);  bool
303    setPointDataView(temp);  DataExpanded::hasNaN() const
304    {
305        const ValueType& v=m_data.getData();
306        for (ValueType::size_type i=0;i<v.size();++i)
307        {
308            if (nancheck(v[i]))
309            {
310                return true;
311            }
312        }
313        return false;
314  }  }
315    
316    
317  string  string
318  DataExpanded::toString() const  DataExpanded::toString() const
319  {  {
320    stringstream temp;    stringstream temp;
321    //    FunctionSpace fs=getFunctionSpace();
322    // create a temporary view as the offset will be changed  
323    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    int offset=0;
324    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
325      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
326        tempView.setOffset(m_data.index(i,j));        offset=getPointOffset(i,j);
327        stringstream suffix;        stringstream suffix;
328        suffix << "(" << i << "," << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
329        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
330        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
331          temp << endl;          temp << endl;
332        }        }
333      }      }
334    }    }
335      string result=temp.str();
336      if (result.empty())
337      {
338        return "(data contains no samples)\n";
339      }
340    return temp.str();    return temp.str();
341  }  }
342    
343  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
344  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
345                               int dataPointNo) const                               int dataPointNo) const
346  {  {
347    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
348  }  }
349    
350  DataArrayView  DataTypes::ValueType::size_type
351  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
352                             int dataPointNo)                               int dataPointNo)
353  {  {
354    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));    return m_data.index(sampleNo,dataPointNo);
   return temp;  
355  }  }
356    
357  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
358  DataExpanded::getLength() const  DataExpanded::getLength() const
359  {  {
360    return m_data.size();    return m_data.size();
361  }  }
362    
363    
364    
365  void  void
366  DataExpanded::setRefValue(int ref,  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
367                            const DataArray& value)    CHECK_FOR_EX_WRITE
 {  
368    //    //
369    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
370    int numSamples = getNumSamples();    int numSamples = getNumSamples();
371    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
372      int dataPointRank = getRank();
373    //    ShapeType dataPointShape = getShape();
374    // Determine the sample number which corresponds to this reference number.    if (numSamples*numDataPointsPerSample > 0) {
375    int sampleNo = -1;       //TODO: global error handling
376    int tempRef = -1;       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
377    for (int n=0; n<numSamples; n++) {            throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
378      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);       }
379      if (tempRef == ref) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
380        sampleNo = n;             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
381        break;       }
382      }       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
383    }       ValueType& vec=getVectorRW();
384    if (sampleNo == -1) {       if (dataPointRank==0) {
385      throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");           vec[offset]=value;
386    }       } else if (dataPointRank==1) {
387            for (int i=0; i<dataPointShape[0]; i++) {
388    for (int n=0; n<numDPPSample; n++) {              vec[offset+i]=value;
389      //          }
390      // Get *each* data-point in the sample in turn.       } else if (dataPointRank==2) {
391      DataArrayView pointView = getDataPoint(sampleNo, n);          for (int i=0; i<dataPointShape[0]; i++) {
392      //             for (int j=0; j<dataPointShape[1]; j++) {
393      // Assign the values in the DataArray to this data-point.                vec[offset+getRelIndex(dataPointShape,i,j)]=value;
394      pointView.copy(value.getView());             }
395            }
396         } else if (dataPointRank==3) {
397            for (int i=0; i<dataPointShape[0]; i++) {
398               for (int j=0; j<dataPointShape[1]; j++) {
399                  for (int k=0; k<dataPointShape[2]; k++) {
400                     vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
401                  }
402               }
403            }
404         } else if (dataPointRank==4) {
405             for (int i=0; i<dataPointShape[0]; i++) {
406               for (int j=0; j<dataPointShape[1]; j++) {
407                 for (int k=0; k<dataPointShape[2]; k++) {
408                   for (int l=0; l<dataPointShape[3]; l++) {
409                      vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
410                   }
411                 }
412               }
413             }
414         }
415    }    }
416  }  }
417    
418  void  void
419  DataExpanded::getRefValue(int ref,  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
420                            DataArray& value)    CHECK_FOR_EX_WRITE
 {  
421    //    //
422    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
423    int numSamples = getNumSamples();    int numSamples = getNumSamples();
424    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
   
425    //    //
426    // Determine the sample number which corresponds to this reference number    // check rank:
427    int sampleNo = -1;    if (value.getRank()!=getRank())
428    int tempRef = -1;         throw DataException("Rank of value does not match Data object rank");
429    for (int n=0; n<numSamples; n++) {    if (numSamples*numDataPointsPerSample > 0) {
430      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);       //TODO: global error handling
431      if (tempRef == ref) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
432        sampleNo = n;            throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
433        break;       }
434         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
435               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
436         }
437         ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
438         ValueType& vec=getVectorRW();
439         vec.copyFromArrayToOffset(value,offset,1);
440      }
441    }
442    
443    void
444    DataExpanded::symmetric(DataAbstract* ev)
445    {
446      int sampleNo,dataPointNo;
447      int numSamples = getNumSamples();
448      int numDataPointsPerSample = getNumDPPSample();
449      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
450      if (temp_ev==0) {
451        throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
452      }
453      const ValueType& vec=getVectorRO();
454      const ShapeType& shape=getShape();
455      ValueType& evVec=temp_ev->getVectorRW();
456      const ShapeType& evShape=temp_ev->getShape();
457      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
458      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
459        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
460             DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
461                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
462      }      }
463    }    }
464    if (sampleNo == -1) {  }
465      throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");  void
466    DataExpanded::nonsymmetric(DataAbstract* ev)
467    {
468      int sampleNo,dataPointNo;
469      int numSamples = getNumSamples();
470      int numDataPointsPerSample = getNumDPPSample();
471      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
472      if (temp_ev==0) {
473        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
474      }
475      const ValueType& vec=getVectorRO();
476      const ShapeType& shape=getShape();
477      ValueType& evVec=temp_ev->getVectorRW();
478      const ShapeType& evShape=temp_ev->getShape();
479      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
480      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
481        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
482             DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
483                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
484        }
485    }    }
   
   //  
   // Get the *first* data-point associated with this sample number.  
   DataArrayView pointView = getDataPoint(sampleNo, 0);  
   
   //  
   // Load the values from this data-point into the DataArray  
   value.getView().copy(pointView);  
486  }  }
487    void
488  int  DataExpanded::trace(DataAbstract* ev, int axis_offset)
 DataExpanded::archiveData(ofstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues) const  
489  {  {
490    return(m_data.archiveData(archiveFile, noValues));    int sampleNo,dataPointNo;
491      int numSamples = getNumSamples();
492      int numDataPointsPerSample = getNumDPPSample();
493      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
494      if (temp_ev==0) {
495        throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
496      }
497      const ValueType& vec=getVectorRO();
498      const ShapeType& shape=getShape();
499      ValueType& evVec=temp_ev->getVectorRW();
500      const ShapeType& evShape=temp_ev->getShape();
501      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
502      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
503        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
504             DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
505                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
506        }
507      }
508  }  }
509    
510  int  void
511  DataExpanded::extractData(ifstream& archiveFile,  DataExpanded::transpose(DataAbstract* ev, int axis_offset)
                           const DataArrayView::ValueType::size_type noValues)  
512  {  {
513    return(m_data.extractData(archiveFile, noValues));    int sampleNo,dataPointNo;
514      int numSamples = getNumSamples();
515      int numDataPointsPerSample = getNumDPPSample();
516      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
517      if (temp_ev==0) {
518        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
519      }
520      const ValueType& vec=getVectorRO();
521      const ShapeType& shape=getShape();
522      ValueType& evVec=temp_ev->getVectorRW();
523      const ShapeType& evShape=temp_ev->getShape();
524      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
525      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
526        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
527             DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
528                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
529        }
530      }
531  }  }
532    
533  void  void
534  DataExpanded::copyAll(const boost::python::numeric::array& value) {  DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
535    //  {
536    // Get the number of samples and data-points per sample.    int sampleNo,dataPointNo;
537    int numSamples = getNumSamples();    int numSamples = getNumSamples();
538    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
539    int dataPointRank = getPointDataView().getRank();    DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
540    ShapeType dataPointShape = getPointDataView().getShape();    if (temp_ev==0) {
541    //      throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
542    // check rank:    }
543    if (value.getrank()!=dataPointRank+1)    const ValueType& vec=getVectorRO();
544         throw DataException("Rank of numarray does not match Data object rank");    const ShapeType& shape=getShape();
545    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)    ValueType& evVec=temp_ev->getVectorRW();
546         throw DataException("leading dimension of numarray is too small");    const ShapeType& evShape=temp_ev->getShape();
547    //    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
548    int dataPoint = 0;    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
549    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
550      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {           DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
551        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
       if (dataPointRank==0) {  
          dataPointView()=extract<double>(value[dataPoint]);  
       } else if (dataPointRank==1) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
             dataPointView(i)=extract<double>(value[dataPoint][i]);  
          }  
       } else if (dataPointRank==2) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
              dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);  
            }  
          }  
        } else if (dataPointRank==3) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
              for (int k=0; k<dataPointShape[2]; k++) {  
                  dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);  
              }  
            }  
          }  
        } else if (dataPointRank==4) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
              for (int k=0; k<dataPointShape[2]; k++) {  
                for (int l=0; l<dataPointShape[3]; l++) {  
                  dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);  
                }  
              }  
            }  
          }  
       }  
       dataPoint++;  
552      }      }
553    }    }
554  }  }
# Line 474  DataExpanded::eigenvalues(DataAbstract* Line 562  DataExpanded::eigenvalues(DataAbstract*
562    if (temp_ev==0) {    if (temp_ev==0) {
563      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
564    }    }
565    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
566    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
567      ValueType& evVec=temp_ev->getVectorRW();
568      const ShapeType& evShape=temp_ev->getShape();
569    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
570    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
571      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
572           DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
573                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
574      }      }
575    }    }
576  }  }
# Line 498  DataExpanded::eigenvalues_and_eigenvecto Line 588  DataExpanded::eigenvalues_and_eigenvecto
588    if (temp_V==0) {    if (temp_V==0) {
589      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
590    }    }
591    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
592    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
593    DataArrayView& VView=V->getPointDataView();    ValueType& evVec=temp_ev->getVectorRW();
594      const ShapeType& evShape=temp_ev->getShape();
595      ValueType& VVec=temp_V->getVectorRW();
596      const ShapeType& VShape=temp_V->getShape();
597    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
598    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
599      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
600           DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
601                                                       evView,ev->getPointOffset(sampleNo,dataPointNo),                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
602                                                       VView,V->getPointOffset(sampleNo,dataPointNo),                                      VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
603                                                       tol);      }
604      }
605    }
606    
607    
608    int
609    DataExpanded::matrixInverse(DataAbstract* out) const
610    {
611      DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
612      if (temp==0)
613      {
614        throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
615      }
616    
617      if (getRank()!=2)
618      {
619        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
620    
621      }
622      int  sampleNo;
623      const int numdpps=getNumDPPSample();
624      const int numSamples = getNumSamples();
625      const ValueType& vec=m_data.getData();
626      int errcode=0;
627      #pragma omp parallel private(sampleNo)
628      {
629         int errorcode=0;
630         LapackInverseHelper h(getShape()[0]);
631         #pragma omp for schedule(static)
632         for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
633         {
634                // not sure I like all those virtual calls to getPointOffset
635            DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
636            int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
637        if (res>errorcode)
638        {
639            errorcode=res;
640            #pragma omp critical
641            {
642              errcode=errorcode;    // I'm not especially concerned which error gets reported as long as one is
643            }
644        }
645         }
646      }
647      return errcode;
648      if (errcode)
649      {
650        DataMaths::matrixInverseError(errcode); // throws exceptions
651      }
652    }
653    
654    void
655    DataExpanded::setToZero(){
656    // TODO: Surely there is a more efficient way to do this????
657    // Why is there no memset here? Parallel issues?
658    // A: This ensures that memory is touched by the correct thread.
659      CHECK_FOR_EX_WRITE
660      int numSamples = getNumSamples();
661      int numDataPointsPerSample = getNumDPPSample();
662      DataTypes::ValueType::size_type n = getNoValues();
663      double* p;
664      int  sampleNo,dataPointNo, i;
665      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
666      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
667        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
668            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
669            for (i=0; i<n ;++i) p[i]=0.;
670      }      }
671    }    }
672  }  }
673    
674    /* Append MPI rank to file name if multiple MPI processes */
675    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
676      /* Make plenty of room for the mpi_rank number and terminating '\0' */
677      char *newFileName = (char *)malloc(strlen(fileName)+20);
678      strncpy(newFileName, fileName, strlen(fileName)+1);
679      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
680      return(newFileName);
681    }
682    
683    void
684    DataExpanded::dump(const std::string fileName) const
685    {
686       #ifdef USE_NETCDF
687       const int ldims=2+DataTypes::maxRank;
688       const NcDim* ncdims[ldims];
689       NcVar *var, *ids;
690       int rank = getRank();
691       int type=  getFunctionSpace().getTypeCode();
692       int ndims =0;
693       long dims[ldims];
694       const double* d_ptr=&(m_data[0]);
695       const DataTypes::ShapeType& shape = getShape();
696       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
697       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
698    #ifdef ESYS_MPI
699       MPI_Status status;
700    #endif
701    
702    #ifdef ESYS_MPI
703       /* Serialize NetCDF I/O */
704       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
705    #endif
706    
707       // netCDF error handler
708       NcError err(NcError::verbose_nonfatal);
709       // Create the file.
710       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
711       NcFile dataFile(newFileName, NcFile::Replace);
712       // check if writing was successful
713       if (!dataFile.is_valid())
714            throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
715       if (!dataFile.add_att("type_id",2) )
716            throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
717       if (!dataFile.add_att("rank",rank) )
718            throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
719       if (!dataFile.add_att("function_space_type",type))
720            throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
721       ndims=rank+2;
722       if ( rank >0 ) {
723           dims[0]=shape[0];
724           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
725                throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
726       }
727       if ( rank >1 ) {
728           dims[1]=shape[1];
729           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
730                throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
731       }
732       if ( rank >2 ) {
733           dims[2]=shape[2];
734           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
735                throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
736       }
737       if ( rank >3 ) {
738           dims[3]=shape[3];
739           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
740                throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
741       }
742       dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
743       if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
744                throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
745       dims[rank+1]=getFunctionSpace().getNumSamples();
746       if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
747                throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
748    
749       if (getFunctionSpace().getNumSamples()>0)
750       {
751    
752         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
753            throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
754         const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
755         if (! (ids->put(ids_p,dims[rank+1])) )
756            throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
757         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
758            throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
759         if (! (var->put(d_ptr,dims)) )
760            throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
761       }
762    #ifdef ESYS_MPI
763       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
764    #endif
765       #else
766       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
767       #endif
768    }
769    
770    void  
771    DataExpanded::setTaggedValue(int tagKey,
772               const DataTypes::ShapeType& pointshape,
773                   const DataTypes::ValueType& value,
774               int dataOffset)
775    {
776      CHECK_FOR_EX_WRITE
777      int numSamples = getNumSamples();
778      int numDataPointsPerSample = getNumDPPSample();
779      int sampleNo,dataPointNo, i;
780      DataTypes::ValueType::size_type n = getNoValues();
781      double* p;
782      const double* in=&value[0+dataOffset];
783      
784      if (value.size() != n) {
785        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
786      }
787    
788      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
789      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
790        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
791            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
792                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
793                for (i=0; i<n ;++i) p[i]=in[i];
794            }
795        }
796      }
797    }
798    
799    
800    void
801    DataExpanded::reorderByReferenceIDs(int *reference_ids)
802    {
803      CHECK_FOR_EX_WRITE
804      int numSamples = getNumSamples();
805      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
806      int sampleNo, sampleNo2,i;
807      double* p,*p2;
808      register double rtmp;
809      FunctionSpace fs=getFunctionSpace();
810    
811      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
812         const int id_in=reference_ids[sampleNo];
813         const int id=fs.getReferenceIDOfSample(sampleNo);
814         if (id!=id_in) {
815             bool matched=false;
816             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
817                  if (id == reference_ids[sampleNo2]) {
818                     p=&(m_data[getPointOffset(sampleNo,0)]);
819                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
820                     for (i=0; i<n ;i++) {
821                             rtmp=p[i];
822                             p[i]=p2[i];
823                             p2[i]=rtmp;
824                     }
825                     reference_ids[sampleNo]=id;
826                     reference_ids[sampleNo2]=id_in;
827                     matched=true;
828                     break;
829                  }
830             }
831             if (! matched) {
832                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
833             }
834         }
835       }
836    }
837    
838    DataTypes::ValueType&
839    DataExpanded::getVectorRW()
840    {
841        CHECK_FOR_EX_WRITE
842        return m_data.getData();
843    }
844    
845    const DataTypes::ValueType&
846    DataExpanded::getVectorRO() const
847    {
848        return m_data.getData();
849    }
850    
851    
852    #ifndef MKLRANDOM
853    namespace {
854        
855    boost::mt19937 base;        // used to seed all the other generators  
856    vector<boost::mt19937> gens;
857    
858    
859    void seedGens(long seed)
860    {
861    #ifdef _OPENMP
862        int numthreads=omp_get_max_threads();
863    #else
864        int numthreads=1;
865    #endif
866        if (gens.size()==0)     // we haven't instantiated the generators yet  
867        {
868            gens.resize(numthreads);    // yes this means all the generators will be owned by one thread in a NUMA sense      
869        }                   // I don't think it is worth a more complicated solution at this point
870        if (seed!=0)
871        {
872           base.seed((uint32_t)seed);   // without this cast, icc gets confused    
873           for (int i=0;i<numthreads;++i)
874           {
875            uint32_t b=base();
876                gens[i].seed(b);    // initialise each generator with successive random values      
877           }      
878        }
879    }
880      
881      
882    }
883    #endif
884    
885    // Idea here is to create an array of seeds by feeding the original seed into the random generator
886    // The code at the beginning of the function to compute the seed if one is given is
887    // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
888    // I make no claim about how well these initial seeds are distributed
889    void DataExpanded::randomFill(long seed)
890    {
891        CHECK_FOR_EX_WRITE
892        static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
893        if (seed==0)        // for each one
894        {
895        if (prevseed==0)
896        {
897            time_t s=time(0);
898            seed=s;
899        }
900        else
901        {
902            seed=prevseed+419;  // these numbers are arbitrary
903            if (seed>3040101)       // I want to avoid overflow on 32bit systems
904            {
905            seed=((int)(seed)%0xABCD)+1;
906            }
907        }
908        }
909        // now we need to consider MPI since we don't want each rank to start with the same seed
910        seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
911        prevseed=seed;
912    
913    #ifdef MKLRANDOM
914    
915    #ifdef _OPENMP
916        int numthreads=omp_get_max_threads();
917    #else
918        int numthreads=1;
919    #endif
920        double* seeds=new double[numthreads];
921        VSLStreamStatePtr sstream;
922    
923        int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed);  // use a Mersenne Twister
924        numeric_limits<double> dlim;
925        vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
926        vslDeleteStream(&sstream);
927        DataVector& dv=getVectorRW();
928        size_t dvsize=dv.size();
929        #pragma omp parallel
930        {
931        int tnum=0;
932        #ifdef _OPENMP
933        tnum=omp_get_thread_num();
934        #endif
935        VSLStreamStatePtr stream;
936        // the 12345 is a hack to give us a better chance of getting different integer seeds.
937            int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345);  // use a Mersenne Twister
938        int bigchunk=(dvsize/numthreads+1);
939        int smallchunk=dvsize-bigchunk*(numthreads-1);
940        int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
941            vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
942            vslDeleteStream(&stream);
943        }
944        delete[] seeds;
945    #else
946        boost::mt19937::result_type RMAX=base.max();
947        seedGens(seed);
948        DataVector&  dv=getVectorRW();
949        long i;
950        const size_t dvsize=dv.size();
951        
952        #pragma omp parallel private(i)
953        {
954        int tnum=0;
955        #ifdef _OPENMP
956        tnum=omp_get_thread_num();
957        #endif
958        boost::mt19937& generator=gens[tnum];
959        
960            #pragma omp for schedule(static)
961            for (i=0;i<dvsize;++i)
962            {
963          
964          
965        
966          
967          
968    #ifdef _WIN32
969            dv[i]=((double)generator())/RMAX;
970    #else
971            dv[i]=((double)generator())/RMAX;
972    #endif
973            }
974        }
975    #endif
976    }
977    
978  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.584  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26