/[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 1513 by gross, Tue Apr 15 08:47:57 2008 UTC revision 2271 by jfenwick, Mon Feb 16 05:08:29 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 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"
# Line 20  Line 20 
20  #ifdef USE_NETCDF  #ifdef USE_NETCDF
21  #include <netcdfcpp.h>  #include <netcdfcpp.h>
22  #endif  #endif
23    #ifdef PASO_MPI
24    #include <mpi.h>
25    #endif
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28    #include "DataMaths.h"
29    
30  using namespace std;  using namespace std;
31  using namespace boost::python;  using namespace boost::python;
32  using namespace boost;  using namespace boost;
33    using namespace escript::DataTypes;
34    
35    
36    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
37    
38    #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());}
39    
40  namespace escript {  namespace escript {
41    
42  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const WrappedArray& value,
43                             const FunctionSpace& what)                             const FunctionSpace& what)
44    : DataAbstract(what)    : parent(what,value.getShape())
45  {  {
   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]));  
   }  
46    //    //
47    // initialise the data array for this object    // initialise the data array for this object
48    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
49    //    //
50    // copy the given value to every data point    // copy the given value to every data point
51    copy(value);    copy(value);
52  }  }
53    
54  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
55    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(), other.getShape()),
56    m_data(other.m_data)    m_data(other.m_data)
57  {  {
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
58  }  }
59    
60  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
61    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
62  {  {
63    //    //
64    // initialise the data array for this object    // initialise the data array for this object
65    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
66    //    //
67    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
68    copy(other.getPointDataView());    copy(other);
69  }  }
70    
71  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& 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    // 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
79    // value from the given DataTagged object    // value from the given DataTagged object
80    int i,j;    int i,j;
81    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
82    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
83    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
84    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
85      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
86        try {        try {
87          getPointDataView().copy(getPointOffset(i,j),             DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
88                                  other.getPointDataView(),                                  other.getVectorRO(),
89                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
90        }        }
91        catch (std::exception& e) {        catch (std::exception& e) {
# Line 96  DataExpanded::DataExpanded(const DataTag Line 96  DataExpanded::DataExpanded(const DataTag
96  }  }
97    
98  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
99                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
100    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101  {  {
102    //    //
103    // get the shape of the slice    // get the shape of the slice
104    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
105    //    //
106    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
107    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
108    //    //
109    // copy the data    // copy the data
110    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
111    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
112    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
113    int i,j;    int i,j;
114    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
115    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
116      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
117        try {        try {
118          getPointDataView().copySlice(getPointOffset(i,j),          DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
119                                       other.getPointDataView(),                                       other.getVectorRO(),
120                         other.getShape(),
121                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
122                                       region_loop_range);                                       region_loop_range);
123        }        }
# Line 127  DataExpanded::DataExpanded(const DataExp Line 128  DataExpanded::DataExpanded(const DataExp
128    }    }
129  }  }
130    
131  DataExpanded::DataExpanded(const DataArrayView& value,  DataExpanded::DataExpanded(const FunctionSpace& what,
132                             const FunctionSpace& what)                             const DataTypes::ShapeType &shape,
133    : DataAbstract(what)                             const DataTypes::ValueType &data)
134  {    : parent(what,shape)
135    //  {
136    // get the shape of the given data value    EsysAssert(data.size()%getNoValues()==0,
137    DataArrayView::ShapeType tempShape=value.getShape();                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
138    //  
139    // initialise this Data object to the shape of the given data value    if (data.size()==getNoValues())
140    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    {
141    //       ValueType& vec=m_data.getData();
142    // copy the given value to every data point       //
143    copy(value);       // create the view of the data
144         initialise(what.getNumSamples(),what.getNumDPPSample());
145         // now we copy this value to all elements
146         for (int i=0;i<getLength();)
147         {
148        for (unsigned int j=0;j<getNoValues();++j,++i)
149        {
150            vec[i]=data[j];
151        }
152         }
153      }
154      else
155      {
156         //
157         // copy the data in the correct format
158         m_data.getData()=data;
159      }
160    
161    
162  }  }
163    
164  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::~DataExpanded()
                            const DataArrayView::ShapeType &shape,  
                            const DataArrayView::ValueType &data)  
   : DataAbstract(what)  
165  {  {
   //  
   // create the view of the data  
   initialise(shape,what.getNumSamples(),what.getNumDPPSample());  
   //  
   // copy the data in the correct format  
   m_data.getData()=data;  
166  }  }
167    
168  DataExpanded::~DataExpanded()  DataAbstract*
169    DataExpanded::deepCopy()
170  {  {
171      return new DataExpanded(*this);
172  }  }
173    
174    
175  DataAbstract*  DataAbstract*
176  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataTypes::RegionType& region) const
177  {  {
178    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
179  }  }
180    
181  void  void
182  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
183                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
184  {  {
185    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
186    if (tempDataExp==0) {    if (tempDataExp==0) {
187      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
188    }    }
189      CHECK_FOR_EX_WRITE
190    //    //
191    // get shape of slice    // get shape of slice
192    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194    //    //
195    // check shape    // check shape
196    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
197      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
198    }    }
199    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
200      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
201          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
202    }    }
203    //    //
204    // copy the data from the slice into this object    // copy the data from the slice into this object
205    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
206    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
207    int i, j;    int i, j;
208      ValueType& vec=getVectorRW();
209      const ShapeType& mshape=getShape();
210      const ValueType& tVec=tempDataExp->getVectorRO();
211      const ShapeType& tShape=tempDataExp->getShape();
212    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
213    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
214      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
215        getPointDataView().copySliceFrom(getPointOffset(i,j),          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
216                                         tempDataExp->getPointDataView(),                                         tVec,
217                           tShape,
218                                         tempDataExp->getPointOffset(i,j),                                         tempDataExp->getPointOffset(i,j),
219                                         region_loop_range);                                         region_loop_range);
220    
221      }      }
222    }    }
223  }  }
224    
225  void  void
226  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
227  {  {
228      EsysAssert((checkShape(getShape(), value.getShape())),
229                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
230    
231    //    //
232    // copy a single value to every data point in this object    // copy a single value to every data point in this object
233    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
# Line 213  DataExpanded::copy(const DataArrayView& Line 236  DataExpanded::copy(const DataArrayView&
236    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
237    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
238      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
239        // 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);  
240      }      }
241    }    }
242  }  }
243    
244  void  void
245  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const WrappedArray& value)
246  {  {
   
   // 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);  
   
   //  
247    // check the input shape matches this shape    // check the input shape matches this shape
248    if (!getPointDataView().checkShape(temp_dataView.getShape())) {    if (!DataTypes::checkShape(getShape(),value.getShape())) {
249      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
250                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
251                          temp_dataView.getShape()));                          value.getShape(),getShape()));
252    }    }
253    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the data  
   copy(temp_dataView);  
254  }  }
255    
256    
257  void  void
258  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
259                           int noDataPointsPerSample)                           int noDataPointsPerSample)
260  {  {
261      if (noSamples==0)     //retain the default empty object
262      {
263         return;
264      }
265    //    //
266    // resize data array to the required size    // resize data array to the required size
267    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
   //  
   // create the data view of the data array  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
268  }  }
269    
270  string  string
# Line 269  DataExpanded::toString() const Line 272  DataExpanded::toString() const
272  {  {
273    stringstream temp;    stringstream temp;
274    FunctionSpace fs=getFunctionSpace();    FunctionSpace fs=getFunctionSpace();
275    //  
276    // create a temporary view as the offset will be changed    int offset=0;
   DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());  
277    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
278      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
279        tempView.setOffset(getPointOffset(i,j));        offset=getPointOffset(i,j);
280        stringstream suffix;        stringstream suffix;
281        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
282        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
283        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
284          temp << endl;          temp << endl;
285        }        }
286      }      }
287    }    }
288      string result=temp.str();
289      if (result.empty())
290      {
291        return "(data contains no samples)\n";
292      }
293    return temp.str();    return temp.str();
294  }  }
295    
296  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
297  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
298                               int dataPointNo) const                               int dataPointNo) const
299  {  {
300    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
301  }  }
302    
303  DataArrayView  DataTypes::ValueType::size_type
304  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
305                             int dataPointNo)                               int dataPointNo)
306  {  {
307    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));    return m_data.index(sampleNo,dataPointNo);
   return temp;  
308  }  }
309    
310  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
311  DataExpanded::getLength() const  DataExpanded::getLength() const
312  {  {
313    return m_data.size();    return m_data.size();
314  }  }
315    
 int  
 DataExpanded::archiveData(ofstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues) const  
 {  
   return(m_data.archiveData(archiveFile, noValues));  
 }  
316    
 int  
 DataExpanded::extractData(ifstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues)  
 {  
   return(m_data.extractData(archiveFile, noValues));  
 }  
317    
318  void  void
319  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
320      CHECK_FOR_EX_WRITE
321    //    //
322    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
323    int numSamples = getNumSamples();    int numSamples = getNumSamples();
324    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
325    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
326    ShapeType dataPointShape = getPointDataView().getShape();    ShapeType dataPointShape = getShape();
327    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
328       //TODO: global error handling       //TODO: global error handling
329       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 337  DataExpanded::copyToDataPoint(const int Line 332  DataExpanded::copyToDataPoint(const int
332       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
333             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
334       }       }
335       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
336         ValueType& vec=getVectorRW();
337       if (dataPointRank==0) {       if (dataPointRank==0) {
338           dataPointView()=value;           vec[offset]=value;
339       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
340          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
341              dataPointView(i)=value;              vec[offset+i]=value;
342          }          }
343       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
344          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
345             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
346                dataPointView(i,j)=value;                vec[offset+getRelIndex(dataPointShape,i,j)]=value;
347             }             }
348          }          }
349       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
350          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
351             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
352                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<dataPointShape[2]; k++) {
353                   dataPointView(i,j,k)=value;                   vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
354                }                }
355             }             }
356          }          }
# Line 363  DataExpanded::copyToDataPoint(const int Line 359  DataExpanded::copyToDataPoint(const int
359             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
360               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
361                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
362                    dataPointView(i,j,k,l)=value;                    vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
363                 }                 }
364               }               }
365             }             }
# Line 371  DataExpanded::copyToDataPoint(const int Line 367  DataExpanded::copyToDataPoint(const int
367       }       }
368    }    }
369  }  }
370    
371  void  void
372  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) {
373      CHECK_FOR_EX_WRITE
374    //    //
375    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
376    int numSamples = getNumSamples();    int numSamples = getNumSamples();
377    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
   int dataPointRank = getPointDataView().getRank();  
   ShapeType dataPointShape = getPointDataView().getShape();  
378    //    //
379    // check rank:    // check rank:
380    if (value.getrank()!=dataPointRank)    if (value.getRank()!=getRank())
381         throw DataException("Rank of numarray does not match Data object rank");         throw DataException("Rank of numarray does not match Data object rank");
382    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
383       //TODO: global error handling       //TODO: global error handling
# Line 391  DataExpanded::copyToDataPoint(const int Line 387  DataExpanded::copyToDataPoint(const int
387       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
388             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
389       }       }
390       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
391       if (dataPointRank==0) {       ValueType& vec=getVectorRW();
392           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++;  
     }  
393    }    }
394  }  }
395    
396  void  void
397  DataExpanded::symmetric(DataAbstract* ev)  DataExpanded::symmetric(DataAbstract* ev)
398  {  {
# Line 489  DataExpanded::symmetric(DataAbstract* ev Line 403  DataExpanded::symmetric(DataAbstract* ev
403    if (temp_ev==0) {    if (temp_ev==0) {
404      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).");
405    }    }
406    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
407    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
408      ValueType& evVec=temp_ev->getVectorRW();
409      const ShapeType& evShape=temp_ev->getShape();
410    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
411    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
412      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
413           DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
414                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
415      }      }
416    }    }
417  }  }
# Line 509  DataExpanded::nonsymmetric(DataAbstract* Line 425  DataExpanded::nonsymmetric(DataAbstract*
425    if (temp_ev==0) {    if (temp_ev==0) {
426      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).");
427    }    }
428    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
429    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
430      ValueType& evVec=temp_ev->getVectorRW();
431      const ShapeType& evShape=temp_ev->getShape();
432    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
433    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
434      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
435           DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
436                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
437      }      }
438    }    }
439  }  }
# Line 529  DataExpanded::trace(DataAbstract* ev, in Line 447  DataExpanded::trace(DataAbstract* ev, in
447    if (temp_ev==0) {    if (temp_ev==0) {
448      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).");
449    }    }
450    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
451    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
452      ValueType& evVec=temp_ev->getVectorRW();
453      const ShapeType& evShape=temp_ev->getShape();
454    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
455    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
456      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
457           DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
458                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
459      }      }
460    }    }
461  }  }
# Line 550  DataExpanded::transpose(DataAbstract* ev Line 470  DataExpanded::transpose(DataAbstract* ev
470    if (temp_ev==0) {    if (temp_ev==0) {
471      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).");
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::transpose(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
481                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
482      }      }
483    }    }
484  }  }
# Line 571  DataExpanded::swapaxes(DataAbstract* ev, Line 493  DataExpanded::swapaxes(DataAbstract* ev,
493    if (temp_ev==0) {    if (temp_ev==0) {
494      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).");
495    }    }
496    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
497    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
498      ValueType& evVec=temp_ev->getVectorRW();
499      const ShapeType& evShape=temp_ev->getShape();
500    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
501    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503           DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
504                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
505      }      }
506    }    }
507  }  }
# Line 591  DataExpanded::eigenvalues(DataAbstract* Line 515  DataExpanded::eigenvalues(DataAbstract*
515    if (temp_ev==0) {    if (temp_ev==0) {
516      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).");
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::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
526                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
527      }      }
528    }    }
529  }  }
# Line 615  DataExpanded::eigenvalues_and_eigenvecto Line 541  DataExpanded::eigenvalues_and_eigenvecto
541    if (temp_V==0) {    if (temp_V==0) {
542      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).");
543    }    }
544    DataArrayView& thisView=getPointDataView();    const ValueType& vec=getVectorRO();
545    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
546    DataArrayView& VView=V->getPointDataView();    ValueType& evVec=temp_ev->getVectorRW();
547      const ShapeType& evShape=temp_ev->getShape();
548      ValueType& VVec=temp_V->getVectorRW();
549      const ShapeType& VShape=temp_V->getShape();
550    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
551    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
552      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
553           DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
554                                                       evView,ev->getPointOffset(sampleNo,dataPointNo),                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
555                                                       VView,V->getPointOffset(sampleNo,dataPointNo),                                      VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
                                                      tol);  
556      }      }
557    }    }
558  }  }
559    
560  void  void
561  DataExpanded::setToZero(){  DataExpanded::setToZero(){
562    // TODO: Surely there is a more efficient way to do this????
563    // Why is there no memset here? Parallel issues?
564      CHECK_FOR_EX_WRITE
565    int numSamples = getNumSamples();    int numSamples = getNumSamples();
566    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
567    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
   DataArrayView::ValueType::size_type n = thisView.noValues();  
568    double* p;    double* p;
569    int  sampleNo,dataPointNo, i;    int  sampleNo,dataPointNo, i;
570    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
# Line 646  DataExpanded::setToZero(){ Line 576  DataExpanded::setToZero(){
576    }    }
577  }  }
578    
579    /* Append MPI rank to file name if multiple MPI processes */
580    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
581      /* Make plenty of room for the mpi_rank number and terminating '\0' */
582      char *newFileName = (char *)malloc(strlen(fileName)+20);
583      strncpy(newFileName, fileName, strlen(fileName)+1);
584      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
585      return(newFileName);
586    }
587    
588  void  void
589  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
590  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
591     #ifdef USE_NETCDF     #ifdef USE_NETCDF
592     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataTypes::maxRank;
593     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
594     NcVar *var, *ids;     NcVar *var, *ids;
595     int rank = getPointDataView().getRank();     int rank = getRank();
596     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
597     int ndims =0;     int ndims =0;
598     long dims[ldims];     long dims[ldims];
599     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
600     DataArrayView::ShapeType shape = getPointDataView().getShape();     const DataTypes::ShapeType& shape = getShape();
601       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
602       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
603    #ifdef PASO_MPI
604       MPI_Status status;
605    #endif
606    
607    #ifdef PASO_MPI
608       /* Serialize NetCDF I/O */
609       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
610    #endif
611    
612     // netCDF error handler     // netCDF error handler
613     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
614     // Create the file.     // Create the file.
615     NcFile dataFile(fileName.c_str(), NcFile::Replace);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
616       NcFile dataFile(newFileName, NcFile::Replace);
617     // check if writing was successful     // check if writing was successful
618     if (!dataFile.is_valid())     if (!dataFile.is_valid())
619          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 627  DataExpanded::dump(const std::string fil
627     if ( rank >0 ) {     if ( rank >0 ) {
628         dims[0]=shape[0];         dims[0]=shape[0];
629         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
630              throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
631     }     }
632     if ( rank >1 ) {     if ( rank >1 ) {
633         dims[1]=shape[1];         dims[1]=shape[1];
634         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
635              throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
636     }     }
637     if ( rank >2 ) {     if ( rank >2 ) {
638         dims[2]=shape[2];         dims[2]=shape[2];
639         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
640              throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
641     }     }
642     if ( rank >3 ) {     if ( rank >3 ) {
643         dims[3]=shape[3];         dims[3]=shape[3];
644         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
645              throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
646     }     }
647     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
648     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 651  DataExpanded::dump(const std::string fil
651     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
652              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
653    
654     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )     if (getFunctionSpace().getNumSamples()>0)
655       {
656    
657         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
658          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
659     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
660     if (! (ids->put(ids_p,dims[rank+1])) )       if (! (ids->put(ids_p,dims[rank+1])) )
661          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
662         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
    if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )  
663          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
664     if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
665          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
666       }
667    #ifdef PASO_MPI
668       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
669    #endif
670     #else     #else
671     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.");
672     #endif     #endif
673  }  }
674    
675  void  void  
676  DataExpanded::setTaggedValue(int tagKey,  DataExpanded::setTaggedValue(int tagKey,
677                               const DataArrayView& value)             const DataTypes::ShapeType& pointshape,
678                   const DataTypes::ValueType& value,
679               int dataOffset)
680  {  {
681      CHECK_FOR_EX_WRITE
682    int numSamples = getNumSamples();    int numSamples = getNumSamples();
683    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
684    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
685    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
686    DataArrayView::ValueType::size_type n = thisView.noValues();    double* p;
687    double* p,*in=&(value.getData()[0]);    const double* in=&value[0+dataOffset];
688        
689    if (value.noValues() != n) {    if (value.size() != n) {
690      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.");
691    }    }
692    
# Line 746  DataExpanded::setTaggedValue(int tagKey, Line 701  DataExpanded::setTaggedValue(int tagKey,
701    }    }
702  }  }
703    
704    
705  void  void
706  DataExpanded::reorderByReferenceIDs(int *reference_ids)  DataExpanded::reorderByReferenceIDs(int *reference_ids)
707  {  {
708      CHECK_FOR_EX_WRITE
709    int numSamples = getNumSamples();    int numSamples = getNumSamples();
710    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
   DataArrayView::ValueType::size_type n = thisView.noValues() * getNumDPPSample();  
711    int sampleNo, sampleNo2,i;    int sampleNo, sampleNo2,i;
712    double* p,*p2;    double* p,*p2;
713    register double rtmp;    register double rtmp;
# Line 777  DataExpanded::reorderByReferenceIDs(int Line 733  DataExpanded::reorderByReferenceIDs(int
733                   break;                   break;
734                }                }
735           }           }
736           if (not matched) {           if (! matched) {
737              throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");              throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
738           }           }
739       }       }
740     }     }
741  }  }
742    
743    DataTypes::ValueType&
744    DataExpanded::getVectorRW()
745    {
746        CHECK_FOR_EX_WRITE
747        return m_data.getData();
748    }
749    
750    const DataTypes::ValueType&
751    DataExpanded::getVectorRO() const
752    {
753        return m_data.getData();
754    }
755    
756    
757  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1513  
changed lines
  Added in v.2271

  ViewVC Help
Powered by ViewVC 1.1.26