/[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

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

Legend:
Removed from v.1387  
changed lines
  Added in v.3552

  ViewVC Help
Powered by ViewVC 1.1.26