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

trunk/esys2/escript/src/Data/DataExpanded.cpp revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC trunk/escript/src/DataExpanded.cpp revision 1800 by ksteube, Thu Sep 18 05:28:20 2008 UTC
# Line 1  Line 1 
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
   
 #include "escript/Data/DataException.h"  
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataArrayView.h"  
1    
2  #include <boost/python/extract.hpp>  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16    #include "Data.h"
17    #include "DataExpanded.h"
18    #include "DataException.h"
19    #include "DataConstant.h"
20    #include "DataTagged.h"
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24    #ifdef PASO_MPI
25    #include <mpi.h>
26    #endif
27    
28  #include <iostream>  #include <boost/python/extract.hpp>
29    #include "DataMaths.h"
30    
31  using namespace std;  using namespace std;
32  using namespace boost::python;  using namespace boost::python;
33  using namespace boost;  using namespace boost;
34    using namespace escript::DataTypes;
35    
36  namespace escript {  namespace escript {
37    
38    DataExpanded::DataExpanded(const boost::python::numeric::array& value, const FunctionSpace& what) : DataAbstract(what)  DataExpanded::DataExpanded(const boost::python::numeric::array& value,
39    {                             const FunctionSpace& what)
40      //cout << "Calling DataExpanded constructor 1." << endl;    : DataAbstract(what,DataTypes::shapeFromNumArray(value))
41      DataArrayView::ShapeType tempShape;  {
42      //  /*  DataTypes::ShapeType tempShape;
43      // extract the shape from the python array    //
44      for (int i=0; i<value.getrank(); ++i) {    // extract the shape of the python numarray
45        cout << extract<int>(value.getshape()[i]) << endl;    for (int i=0; i<value.getrank(); i++) {
46        tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
47      }    }*/
48      initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    //
49      //    // initialise the data array for this object
50      // copy the given value to every data point    initialise(what.getNumSamples(),what.getNumDPPSample());
51      copy(value);    //
52    }    // copy the given value to every data point
53      copy(value);
54    DataExpanded::DataExpanded(const DataExpanded& other) : DataAbstract(other.getFunctionSpace()), m_data(other.m_data)  }
55    {      
56      //cout << "Calling DataExpanded copy constructor." << endl;  DataExpanded::DataExpanded(const DataExpanded& other)
57      //    : DataAbstract(other.getFunctionSpace(), other.getShape()),
58      // create the view for the data    m_data(other.m_data)
59      DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  {
60      setPointDataView(temp);  /*
61    }    //
62      // create the view for the data
63    DataExpanded::DataExpanded(const DataConstant& other) : DataAbstract(other.getFunctionSpace())    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
64    {        setPointDataView(temp);
65      //cout << "Creating DataExpanded from DataConstant." << endl;  
66      initialise(other.getPointDataView().getShape(),other.getNumSamples(), other.getNumDPPSample());    // keep local shape in sync
67      //    setShape(other.getPointDataView().getShape());*/
68      // DataConstant only has one value  }
69      copy(other.getPointDataView());  
70    }  DataExpanded::DataExpanded(const DataConstant& other)
71      : DataAbstract(other.getFunctionSpace(), other.getShape())
72    DataExpanded::DataExpanded(const DataTagged& other) : DataAbstract(other.getFunctionSpace())  {
73    {        //
74      //cout << "Creating DataExpanded from DataTagged." << endl;    // initialise the data array for this object
75      initialise(other.getPointDataView().getShape(),other.getNumSamples(), other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
76      int i,j;    //
77      DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    // DataConstant only has one value, copy this to every data point
78      DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    copy(other);
79  #pragma omp parallel for private(i,j) schedule(static)  }
80      for (i=0;i<numRows;++i) {  
81        for (j=0;j<numCols;++j) {  DataExpanded::DataExpanded(const DataTagged& other)
82      try {    : DataAbstract(other.getFunctionSpace(), other.getShape())
83        getPointDataView().copy(getPointOffset(i,j), other.getPointDataView(), other.getPointOffset(i,j));  {
84      }    //
85      catch (std::exception& e) {    // initialise the data array for this object
86        cout << e.what() << endl;    initialise(other.getNumSamples(),other.getNumDPPSample());
87      }    //
88      // for each data point in this object, extract and copy the corresponding data
89      // value from the given DataTagged object
90      int i,j;
91      DataTypes::ValueType::size_type numRows=m_data.getNumRows();
92      DataTypes::ValueType::size_type numCols=m_data.getNumCols();
93      #pragma omp parallel for private(i,j) schedule(static)
94      for (i=0;i<numRows;i++) {
95        for (j=0;j<numCols;j++) {
96          try {
97    //         getPointDataView().copy(getPointOffset(i,j),
98    //                                 other.getPointDataView(),
99    //                                 other.getPointOffset(i,j));
100               DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
101                                    other.getVector(),
102                                    other.getPointOffset(i,j));
103        }        }
104      }        catch (std::exception& e) {
105    }          cout << e.what() << endl;
   
   DataExpanded::DataExpanded(const DataExpanded& other, const DataArrayView::RegionType& region) : DataAbstract(other.getFunctionSpace())  
   {  
     //cout << "Creating DataExpanded by slicing from another DataExpanded." << endl;  
     //  
     // get the shape of the slice  
     DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  
     initialise(shape,other.getNumSamples(),other.getNumDPPSample());  
     //  
     // copy the data  
     DataArrayView::ValueType::size_type numRows=m_data.getNumRows();  
     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();  
     int i,j;  
 #pragma omp parallel for private(i,j) schedule(static)  
     for (i=0;i<numRows;++i) {  
       for (j=0;j<numCols;++j) {  
     try {  
       getPointDataView().copySlice(getPointOffset(i,j), other.getPointDataView(), other.getPointOffset(i,j), region);  
     }  
     catch (std::exception& e) {  
       cout << e.what() << endl;  
     }  
106        }        }
107      }      }
108    }    }
109    }
110    
111    DataExpanded::DataExpanded(const DataArrayView& value,const FunctionSpace& what): DataAbstract(what)  DataExpanded::DataExpanded(const DataExpanded& other,
112    {                             const DataTypes::RegionType& region)
113      //cout << "Calling DataExpanded constructor 2." << endl;    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
114      DataArrayView::ShapeType tempShape=value.getShape();  {
115      initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    //
116      copy(value);    // get the shape of the slice
117    }  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
118      //
119    DataExpanded::~DataExpanded()    // initialise this Data object to the shape of the slice
120    {    initialise(other.getNumSamples(),other.getNumDPPSample());
121      //cout << "Destructing DataExpanded." << endl;    //
122    }    // copy the data
123      DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
124    DataAbstract* DataExpanded::getSlice(const DataArrayView::RegionType& region) const    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
125    {    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
126      return new DataExpanded(*this,region);    int i,j;
127    }    #pragma omp parallel for private(i,j) schedule(static)
128      for (i=0;i<numRows;i++) {
129    void DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)      for (j=0;j<numCols;j++) {
130    {        try {
131      //  //         getPointDataView().copySlice(getPointOffset(i,j),
132      // reshape a rank zero data point  //                                      other.getPointDataView(),
133      if (getPointDataView().getRank()!=0) {  //                                      other.getPointOffset(i,j),
134        stringstream temp;  //                                      region_loop_range);
135        temp << "Error - Can only reshape Data with data points of rank 0. "          DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
136         << "This Data has data points with rank: "                                       other.getVector(),
137         << getPointDataView().getRank();                       other.getShape(),
138        throw DataException(temp.str());                                       other.getPointOffset(i,j),
139      }                                       region_loop_range);
     DataBlocks2D newData(getNumSamples(),getNumDPPSample(), DataArrayView::noValues(shape));  
     DataArrayView newView(newData.getData(),shape);  
     //  
     // Copy the original data to every value for the new shape  
     int i,j;  
     int nRows=m_data.getNumRows();  
     int nCols=m_data.getNumCols();  
     #pragma omp parallel for private(i,j) schedule(static)  
     for (i=0;i<nRows;++i) {  
       for (j=0;j<nCols;++j) {  
     //  
     // Copy the data into the specified offset  
     // NOTE: An exception may be thown from this call if  
     // DOASSERT is on which of course will play  
     // havoc with the omp threads. Run single threaded  
     // if using DOASSERT.  
     newView.copy(newData.index(i,j),m_data.getData()[m_data.index(i,j)]);  
140        }        }
141      }        catch (std::exception& e) {
142      m_data.Swap(newData);          cout << e.what() << endl;
     DataArrayView temp(m_data.getData(),shape);  
     setPointDataView(temp);  
   }  
   
   void DataExpanded::setSlice(const DataAbstract* value, const DataArrayView::RegionType& region)  
   {  
     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);  
     if (tempDataExp==0) {  
       throw DataException("Programming error - casting to DataExpanded.");  
     }  
     if (getPointDataView().getRank()!=region.size()) {  
       throw DataException("Error - Invalid slice region.");  
     }  
     if (!value->getPointDataView().checkShape(DataArrayView::getResultSliceShape(region))) {  
       throw DataException (value->getPointDataView().createShapeErrorMessage(  
                 "Error - Couldn't copy slice due to shape mismatch.",DataArrayView::getResultSliceShape(region)));  
     }  
     //  
     // copy the data  
     DataArrayView::ValueType::size_type numRows=m_data.getNumRows();  
     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();  
     int i,j;  
 #pragma omp parallel for private(i,j) schedule(static)  
     for (i=0;i<numRows;++i) {  
       for (j=0;j<numCols;++j) {  
     getPointDataView().copySliceFrom(getPointOffset(i,j), tempDataExp->getPointDataView(), tempDataExp->getPointOffset(i,j), region);  
143        }        }
144      }      }
145    }    }
146    }
147    
148    void DataExpanded::copy(const DataArrayView& value)  // DataExpanded::DataExpanded(const DataArrayView& value,
149    //                            const FunctionSpace& what)
150    //   : DataAbstract(what)
151    // {
152    //   //
153    //   // get the shape of the given data value
154    //   DataTypes::ShapeType tempShape=value.getShape();
155    //   //
156    //   // initialise this Data object to the shape of the given data value
157    //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
158    //   //
159    //   // copy the given value to every data point
160    //   copy(value);
161    // }
162    
163    DataExpanded::DataExpanded(const FunctionSpace& what,
164                               const DataTypes::ShapeType &shape,
165                               const DataTypes::ValueType &data)
166      : DataAbstract(what,shape)
167    {
168      EsysAssert(data.size()%getNoValues()==0,
169                     "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
170    
171      if (data.size()==getNoValues())
172      {
173         ValueType& vec=m_data.getData();
174         //
175         // create the view of the data
176         initialise(what.getNumSamples(),what.getNumDPPSample());
177         // now we copy this value to all elements
178         for (int i=0;i<getLength();)
179         {
180        for (int j=0;j<getNoValues();++j,++i)
181        {
182            vec[i]=data[j];
183        }
184         }
185      }
186      else
187    {    {
188      //       //
189      // Copy a single value to every data point       // copy the data in the correct format
190      int i,j;       m_data.getData()=data;
191      int nRows=m_data.getNumRows();    }
192      int nCols=m_data.getNumCols();  
193      #pragma omp parallel for private(i,j) schedule(static)  
194      for (i=0;i<nRows;++i) {  }
195        for (j=0;j<nCols;++j) {  
196      //  DataExpanded::~DataExpanded()
197      // Copy the data into the specified offset  {
198      // NOTE: An exception may be thown from this call if  }
199      // DOASSERT is on which of course will play  
200      // havoc with the omp threads. Run single threaded  DataAbstract*
201      // if using DOASSERT.  DataExpanded::deepCopy()
202      getPointDataView().copy(m_data.index(i,j),value);  {
203      return new DataExpanded(*this);
204    }
205    
206    
207    DataAbstract*
208    DataExpanded::getSlice(const DataTypes::RegionType& region) const
209    {
210      return new DataExpanded(*this,region);
211    }
212    
213    void
214    DataExpanded::setSlice(const DataAbstract* value,
215                           const DataTypes::RegionType& region)
216    {
217      const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
218      if (tempDataExp==0) {
219        throw DataException("Programming error - casting to DataExpanded.");
220      }
221      //
222      // get shape of slice
223      DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
224      DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
225      //
226      // check shape
227      if (getRank()!=region.size()) {
228        throw DataException("Error - Invalid slice region.");
229      }
230      if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
231        throw DataException (DataTypes::createShapeErrorMessage(
232            "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
233      }
234      //
235      // copy the data from the slice into this object
236      DataTypes::ValueType::size_type numRows=m_data.getNumRows();
237      DataTypes::ValueType::size_type numCols=m_data.getNumCols();
238      int i, j;
239      ValueType& vec=getVector();
240      const ShapeType& mshape=getShape();
241      const ValueType& tVec=tempDataExp->getVector();
242      const ShapeType& tShape=tempDataExp->getShape();
243      #pragma omp parallel for private(i,j) schedule(static)
244      for (i=0;i<numRows;i++) {
245        for (j=0;j<numCols;j++) {
246    /*      getPointDataView().copySliceFrom(getPointOffset(i,j),
247                                           tempDataExp->getPointDataView(),
248                                           tempDataExp->getPointOffset(i,j),
249                                           region_loop_range);*/
250            DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
251                                           tVec,
252                           tShape,
253                                           tempDataExp->getPointOffset(i,j),
254                                           region_loop_range);
255    
256        }
257      }
258    }
259    
260    void
261    DataExpanded::copy(const DataConstant& value)
262    {
263      EsysAssert((checkShape(getShape(), value.getShape())),
264                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
265    
266      //
267      // copy a single value to every data point in this object
268      int nRows=m_data.getNumRows();
269      int nCols=m_data.getNumCols();
270      int i,j;
271      #pragma omp parallel for private(i,j) schedule(static)
272      for (i=0;i<nRows;i++) {
273        for (j=0;j<nCols;j++) {
274          // NOTE: An exception may be thown from this call if
275          // DOASSERT is on which of course will play
276          // havoc with the omp threads. Run single threaded
277          // if using DOASSERT.
278          //getPointDataView().copy(getPointOffset(i,j),value);
279          DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
280        }
281      }
282    }
283    
284    
285    // void
286    // DataExpanded::copy(const DataArrayView& value)
287    // {
288    //   //
289    //   // copy a single value to every data point in this object
290    //   int nRows=m_data.getNumRows();
291    //   int nCols=m_data.getNumCols();
292    //   int i,j;
293    //   #pragma omp parallel for private(i,j) schedule(static)
294    //   for (i=0;i<nRows;i++) {
295    //     for (j=0;j<nCols;j++) {
296    //       // NOTE: An exception may be thown from this call if
297    //       // DOASSERT is on which of course will play
298    //       // havoc with the omp threads. Run single threaded
299    //       // if using DOASSERT.
300    //       getPointDataView().copy(getPointOffset(i,j),value);
301    //     }
302    //   }
303    // }
304    
305    void
306    DataExpanded::copy(const boost::python::numeric::array& value)
307    {
308    
309      // extract the shape of the numarray
310      DataTypes::ShapeType tempShape;
311      for (int i=0; i < value.getrank(); i++) {
312        tempShape.push_back(extract<int>(value.getshape()[i]));
313      }
314    
315      // get the space for the data vector
316    //   int len = DataTypes::noValues(tempShape);
317    //   DataVector temp_data(len, 0.0, len);
318    //   DataArrayView temp_dataView(temp_data, tempShape);
319    //   temp_dataView.copy(value);
320    
321      //
322      // check the input shape matches this shape
323      if (!DataTypes::checkShape(getShape(),tempShape)) {
324        throw DataException(DataTypes::createShapeErrorMessage(
325                            "Error - (DataExpanded) Cannot copy due to shape mismatch.",
326                            tempShape,getShape()));
327      }
328      //
329      // now copy over the data
330      //copy(temp_dataView);
331      getVector().copyFromNumArray(value);
332    }
333    
334    
335    void
336    DataExpanded::initialise(int noSamples,
337                             int noDataPointsPerSample)
338    {
339      //
340      // resize data array to the required size
341      m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
342    
343      //
344    //   // create the data view of the data array
345    //   DataArrayView temp(m_data.getData(),shape);
346    //   setPointDataView(temp);
347    
348    //   // keep shape in sync
349    //   setShape(shape);
350    }
351    
352    string
353    DataExpanded::toString() const
354    {
355      stringstream temp;
356      FunctionSpace fs=getFunctionSpace();
357      //
358      // create a temporary view as the offset will be changed
359    //  DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
360    
361      int offset=0;
362      for (int i=0;i<m_data.getNumRows();i++) {
363        for (int j=0;j<m_data.getNumCols();j++) {
364          offset=getPointOffset(i,j);
365          stringstream suffix;
366          suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
367          temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
368          if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
369            temp << endl;
370        }        }
371      }      }
372    }    }
373      return temp.str();
374    }
375    
376    void DataExpanded::copy(const boost::python::numeric::array& value)  DataTypes::ValueType::size_type
377    {  DataExpanded::getPointOffset(int sampleNo,
378      //                               int dataPointNo) const
379      // first convert the numarray into a DataArrayView format  {
380      DataArray temp(value);    return m_data.index(sampleNo,dataPointNo);
381      //  }
382      // check the input shape matches this shape, this will throw an exception  
383      if (!getPointDataView().checkShape(temp.getView().getShape())) {  // DataArrayView
384        throw DataException(getPointDataView().createShapeErrorMessage(  // DataExpanded::getDataPoint(int sampleNo,
385           "Error - (DataExpanded) Cannot copy due to shape mismatch.",  //                            int dataPointNo)
386           temp.getView().getShape()));  // {
387      }  //   DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
388      //  //   return temp;
389      // now copy over the entire data structure  // }
390      copy(temp.getView());  
391    }  DataTypes::ValueType::size_type
392    DataExpanded::getLength() const
393    void DataExpanded::initialise(const DataArrayView::ShapeType& shape, int noSamples, int noDataPointsPerSample)  {
394    {    return m_data.size();
395      //  }
396      // resize data to the required size  
397      //cout << "noDataPointsPerSample=" << noDataPointsPerSample  
398      //     << " noSamples=" << noSamples  
399      //     << " noValuesPerDataPoint=" << DataArrayView::noValues(shape) << endl;  void
400      m_data.resize(noSamples,noDataPointsPerSample, DataArrayView::noValues(shape));  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
401      //    //
402      // create a point data viewer of the data    // Get the number of samples and data-points per sample.
403      DataArrayView temp(m_data.getData(),shape);    int numSamples = getNumSamples();
404      setPointDataView(temp);    int numDataPointsPerSample = getNumDPPSample();
405      //cout << "leaving initialise." << endl;    int dataPointRank = getRank();
406    }    ShapeType dataPointShape = getShape();
407      if (numSamples*numDataPointsPerSample > 0) {
408    string DataExpanded::toString() const       //TODO: global error handling
409    {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
410      stringstream temp;            throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
411      //       }
412      // create a temporary view as the offset will be changed       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
413      DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape(), getPointDataView().getOffset());             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
414      for (int i=0;i<m_data.getNumRows();++i) {       }
415        for (int j=0;j<m_data.getNumCols();++j) {       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
416          tempView.setOffset(m_data.index(i,j));       ValueType& vec=getVector();
417      stringstream suffix;       if (dataPointRank==0) {
418      suffix << "(" << i << "," << j << ")";           vec[0]=value;
419      temp << tempView.toString(suffix.str());       } else if (dataPointRank==1) {
420      if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {          for (int i=0; i<dataPointShape[0]; i++) {
421        temp << endl;              vec[offset+i]=value;
422      }          }
423         } else if (dataPointRank==2) {
424            for (int i=0; i<dataPointShape[0]; i++) {
425               for (int j=0; j<dataPointShape[1]; j++) {
426                  vec[offset+getRelIndex(dataPointShape,i,j)]=value;
427               }
428            }
429         } else if (dataPointRank==3) {
430            for (int i=0; i<dataPointShape[0]; i++) {
431               for (int j=0; j<dataPointShape[1]; j++) {
432                  for (int k=0; k<dataPointShape[2]; k++) {
433                     vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
434                  }
435               }
436            }
437         } else if (dataPointRank==4) {
438             for (int i=0; i<dataPointShape[0]; i++) {
439               for (int j=0; j<dataPointShape[1]; j++) {
440                 for (int k=0; k<dataPointShape[2]; k++) {
441                   for (int l=0; l<dataPointShape[3]; l++) {
442                      vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
443                   }
444                 }
445               }
446             }
447         }
448      }
449    }
450    void
451    DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
452      //
453      // Get the number of samples and data-points per sample.
454      int numSamples = getNumSamples();
455      int numDataPointsPerSample = getNumDPPSample();
456      int dataPointRank = getRank();
457      const ShapeType& shape = getShape();
458      //
459      // check rank:
460      if (value.getrank()!=dataPointRank)
461           throw DataException("Rank of numarray does not match Data object rank");
462      if (numSamples*numDataPointsPerSample > 0) {
463         //TODO: global error handling
464         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
465              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
466         }
467         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
468               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
469         }
470     //    DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
471         ValueType& vec=getVector();
472         if (dataPointRank==0) {
473             vec[0]=extract<double>(value[0]);
474         } else if (dataPointRank==1) {
475            for (int i=0; i<shape[0]; i++) {
476                vec[i]=extract<double>(value[i]);
477            }
478         } else if (dataPointRank==2) {
479            for (int i=0; i<shape[0]; i++) {
480               for (int j=0; j<shape[1]; j++) {
481                  vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
482               }
483            }
484         } else if (dataPointRank==3) {
485            for (int i=0; i<shape[0]; i++) {
486               for (int j=0; j<shape[1]; j++) {
487                  for (int k=0; k<shape[2]; k++) {
488                     vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
489                  }
490               }
491            }
492         } else if (dataPointRank==4) {
493             for (int i=0; i<shape[0]; i++) {
494               for (int j=0; j<shape[1]; j++) {
495                 for (int k=0; k<shape[2]; k++) {
496                   for (int l=0; l<shape[3]; l++) {
497                      vec[getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
498                   }
499                 }
500               }
501             }
502         }
503      }
504    }
505    void
506    DataExpanded::copyAll(const boost::python::numeric::array& value) {
507      //
508      // Get the number of samples and data-points per sample.
509      int numSamples = getNumSamples();
510      int numDataPointsPerSample = getNumDPPSample();
511      int dataPointRank = getRank();
512      const ShapeType& dataPointShape = getShape();
513      //
514      // check rank:
515      if (value.getrank()!=dataPointRank+1)
516           throw DataException("Rank of numarray does not match Data object rank");
517      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
518           throw DataException("leading dimension of numarray is too small");
519      //
520      ValueType& vec=getVector();
521      int dataPoint = 0;
522      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
523        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
524    //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
525          ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
526          if (dataPointRank==0) {
527             vec[offset]=extract<double>(value[dataPoint]);
528          } else if (dataPointRank==1) {
529             for (int i=0; i<dataPointShape[0]; i++) {
530                vec[offset+i]=extract<double>(value[dataPoint][i]);
531             }
532          } else if (dataPointRank==2) {
533             for (int i=0; i<dataPointShape[0]; i++) {
534               for (int j=0; j<dataPointShape[1]; j++) {
535             vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
536    //             dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
537               }
538             }
539           } else if (dataPointRank==3) {
540             for (int i=0; i<dataPointShape[0]; i++) {
541               for (int j=0; j<dataPointShape[1]; j++) {
542                 for (int k=0; k<dataPointShape[2]; k++) {
543             vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
544                 }
545               }
546             }
547           } else if (dataPointRank==4) {
548             for (int i=0; i<dataPointShape[0]; i++) {
549               for (int j=0; j<dataPointShape[1]; j++) {
550                 for (int k=0; k<dataPointShape[2]; k++) {
551                   for (int l=0; l<dataPointShape[3]; l++) {
552                     vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
553                   }
554                 }
555               }
556             }
557        }        }
558          dataPoint++;
559      }      }
     return temp.str();  
   }  
   
   DataArrayView::ValueType::size_type DataExpanded::getPointOffset(int sampleNo, int dataPointNo) const  
   {  
     return m_data.index(sampleNo,dataPointNo);  
   }  
   
   DataArrayView DataExpanded::getDataPoint(int sampleNo, int dataPointNo)  
   {  
     DataArrayView temp(m_data.getData(), getPointDataView().getShape(), m_data.index(sampleNo,dataPointNo));  
     return temp;  
   }  
   
   DataArrayView::ValueType::size_type DataExpanded::getLength() const  
   {  
     return m_data.getData().size();  
560    }    }
561    }
562    void
563    DataExpanded::symmetric(DataAbstract* ev)
564    {
565      int sampleNo,dataPointNo;
566      int numSamples = getNumSamples();
567      int numDataPointsPerSample = getNumDPPSample();
568      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
569      if (temp_ev==0) {
570        throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
571      }
572    //  DataArrayView& thisView=getPointDataView();
573    //  DataArrayView& evView=ev->getPointDataView();
574      ValueType& vec=getVector();
575      const ShapeType& shape=getShape();
576      ValueType& evVec=temp_ev->getVector();
577      const ShapeType& evShape=temp_ev->getShape();
578      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
579      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
580        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
581             DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
582                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
583        }
584      }
585    }
586    void
587    DataExpanded::nonsymmetric(DataAbstract* ev)
588    {
589      int sampleNo,dataPointNo;
590      int numSamples = getNumSamples();
591      int numDataPointsPerSample = getNumDPPSample();
592      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
593      if (temp_ev==0) {
594        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
595      }
596    //   DataArrayView& thisView=getPointDataView();
597    //   DataArrayView& evView=ev->getPointDataView();
598      ValueType& vec=getVector();
599      const ShapeType& shape=getShape();
600      ValueType& evVec=temp_ev->getVector();
601      const ShapeType& evShape=temp_ev->getShape();
602      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
603      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
604        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
605    //          DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
606    //                                     evView,ev->getPointOffset(sampleNo,dataPointNo));
607             DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
608                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
609        }
610      }
611    }
612    void
613    DataExpanded::trace(DataAbstract* ev, int axis_offset)
614    {
615      int sampleNo,dataPointNo;
616      int numSamples = getNumSamples();
617      int numDataPointsPerSample = getNumDPPSample();
618      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
619      if (temp_ev==0) {
620        throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
621      }
622      ValueType& vec=getVector();
623      const ShapeType& shape=getShape();
624      ValueType& evVec=temp_ev->getVector();
625      const ShapeType& evShape=temp_ev->getShape();
626      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
627      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
628        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
629    /*         DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
630                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);*/
631             DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
632                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
633        }
634      }
635    }
636    
637    void
638    DataExpanded::transpose(DataAbstract* ev, int axis_offset)
639    {
640      int sampleNo,dataPointNo;
641      int numSamples = getNumSamples();
642      int numDataPointsPerSample = getNumDPPSample();
643      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
644      if (temp_ev==0) {
645        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
646      }
647      ValueType& vec=getVector();
648      const ShapeType& shape=getShape();
649      ValueType& evVec=temp_ev->getVector();
650      const ShapeType& evShape=temp_ev->getShape();
651      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
652      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
653        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
654             DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
655                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
656        }
657      }
658    }
659    
660    void
661    DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
662    {
663      int sampleNo,dataPointNo;
664      int numSamples = getNumSamples();
665      int numDataPointsPerSample = getNumDPPSample();
666      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
667      if (temp_ev==0) {
668        throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
669      }
670      ValueType& vec=getVector();
671      const ShapeType& shape=getShape();
672      ValueType& evVec=temp_ev->getVector();
673      const ShapeType& evShape=temp_ev->getShape();
674      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
675      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
676        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
677             DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
678                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
679        }
680      }
681    }
682    void
683    DataExpanded::eigenvalues(DataAbstract* ev)
684    {
685      int sampleNo,dataPointNo;
686      int numSamples = getNumSamples();
687      int numDataPointsPerSample = getNumDPPSample();
688      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
689      if (temp_ev==0) {
690        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
691      }
692      ValueType& vec=getVector();
693      const ShapeType& shape=getShape();
694      ValueType& evVec=temp_ev->getVector();
695      const ShapeType& evShape=temp_ev->getShape();
696      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
697      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
698        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
699             DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
700                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
701        }
702      }
703    }
704    void
705    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
706    {
707      int numSamples = getNumSamples();
708      int numDataPointsPerSample = getNumDPPSample();
709      int sampleNo,dataPointNo;
710      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
711      if (temp_ev==0) {
712        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
713      }
714      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
715      if (temp_V==0) {
716        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
717      }
718      ValueType& vec=getVector();
719      const ShapeType& shape=getShape();
720      ValueType& evVec=temp_ev->getVector();
721      const ShapeType& evShape=temp_ev->getShape();
722      ValueType& VVec=temp_V->getVector();
723      const ShapeType& VShape=temp_V->getShape();
724      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
725      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
726        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
727             DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
728                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
729                                        VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
730        }
731      }
732    }
733    
734    void
735    DataExpanded::setToZero(){
736    // TODO: Surely there is a more efficient way to do this????
737    // Why is there no memset here? Parallel issues?
738      int numSamples = getNumSamples();
739      int numDataPointsPerSample = getNumDPPSample();
740      DataTypes::ValueType::size_type n = getNoValues();
741      double* p;
742      int  sampleNo,dataPointNo, i;
743      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
744      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
745        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
746            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
747            for (i=0; i<n ;++i) p[i]=0.;
748        }
749      }
750    }
751    
752    /* Append MPI rank to file name if multiple MPI processes */
753    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
754      /* Make plenty of room for the mpi_rank number and terminating '\0' */
755      char *newFileName = (char *)malloc(strlen(fileName)+20);
756      strncpy(newFileName, fileName, strlen(fileName)+1);
757      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
758      return(newFileName);
759    }
760    
761    void
762    DataExpanded::dump(const std::string fileName) const
763    {
764       #ifdef USE_NETCDF
765       const int ldims=2+DataTypes::maxRank;
766       const NcDim* ncdims[ldims];
767       NcVar *var, *ids;
768       int rank = getRank();
769       int type=  getFunctionSpace().getTypeCode();
770       int ndims =0;
771       long dims[ldims];
772       const double* d_ptr=&(m_data[0]);
773       const DataTypes::ShapeType& shape = getShape();
774       int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
775       int mpi_num=getFunctionSpace().getDomain().getMPISize();
776    #ifdef PASO_MPI
777       MPI_Status status;
778    #endif
779    
780       // netCDF error handler
781       NcError err(NcError::verbose_nonfatal);
782       // Create the file.
783    #ifdef PASO_MPI
784       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, getFunctionSpace().getDomain().getMPIComm(), &status);
785    #endif
786       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
787       NcFile dataFile(newFileName, NcFile::Replace);
788       // check if writing was successful
789       if (!dataFile.is_valid())
790            throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
791       if (!dataFile.add_att("type_id",2) )
792            throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
793       if (!dataFile.add_att("rank",rank) )
794            throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
795       if (!dataFile.add_att("function_space_type",type))
796            throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
797       ndims=rank+2;
798       if ( rank >0 ) {
799           dims[0]=shape[0];
800           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
801                throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
802       }
803       if ( rank >1 ) {
804           dims[1]=shape[1];
805           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
806                throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
807       }
808       if ( rank >2 ) {
809           dims[2]=shape[2];
810           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
811                throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
812       }
813       if ( rank >3 ) {
814           dims[3]=shape[3];
815           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
816                throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
817       }
818       dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
819       if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
820                throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
821       dims[rank+1]=getFunctionSpace().getNumSamples();
822       if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
823                throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
824    
825       if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
826            throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
827       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
828       if (! (ids->put(ids_p,dims[rank+1])) )
829            throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
830    
831       if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
832            throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
833       if (! (var->put(d_ptr,dims)) )
834            throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
835    #ifdef PASO_MPI
836       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, getFunctionSpace().getDomain().getMPIComm());
837    #endif
838       #else
839       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
840       #endif
841    }
842    
843    // void
844    // DataExpanded::setTaggedValue(int tagKey,
845    //                              const DataArrayView& value)
846    // {
847    //   int numSamples = getNumSamples();
848    //   int numDataPointsPerSample = getNumDPPSample();
849    //   int sampleNo,dataPointNo, i;
850    //   DataTypes::ValueType::size_type n = getNoValues();
851    //   double* p,*in=&(value.getData()[0]);
852    //  
853    //   if (value.noValues() != n) {
854    //     throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
855    //   }
856    //
857    //   #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
858    //   for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
859    //     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
860    //         for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
861    //             p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
862    //             for (i=0; i<n ;++i) p[i]=in[i];
863    //         }
864    //     }
865    //   }
866    // }
867    
868    void  
869    DataExpanded::setTaggedValue(int tagKey,
870               const DataTypes::ShapeType& pointshape,
871                   const DataTypes::ValueType& value,
872               int dataOffset)
873    {
874      int numSamples = getNumSamples();
875      int numDataPointsPerSample = getNumDPPSample();
876      int sampleNo,dataPointNo, i;
877      DataTypes::ValueType::size_type n = getNoValues();
878      double* p;
879      const double* in=&value[0+dataOffset];
880      
881      if (value.size() != n) {
882        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
883      }
884    
885      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
886      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
887        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
888            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
889                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
890                for (i=0; i<n ;++i) p[i]=in[i];
891            }
892        }
893      }
894    }
895    
896    
897    void
898    DataExpanded::reorderByReferenceIDs(int *reference_ids)
899    {
900      int numSamples = getNumSamples();
901      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
902      int sampleNo, sampleNo2,i;
903      double* p,*p2;
904      register double rtmp;
905      FunctionSpace fs=getFunctionSpace();
906    
907      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
908         const int id_in=reference_ids[sampleNo];
909         const int id=fs.getReferenceIDOfSample(sampleNo);
910         if (id!=id_in) {
911             bool matched=false;
912             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
913                  if (id == reference_ids[sampleNo2]) {
914                     p=&(m_data[getPointOffset(sampleNo,0)]);
915                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
916                     for (i=0; i<n ;i++) {
917                             rtmp=p[i];
918                             p[i]=p2[i];
919                             p2[i]=rtmp;
920                     }
921                     reference_ids[sampleNo]=id;
922                     reference_ids[sampleNo2]=id_in;
923                     matched=true;
924                     break;
925                  }
926             }
927             if (! matched) {
928                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
929             }
930         }
931       }
932    }
933    
934    DataTypes::ValueType&
935    DataExpanded::getVector()
936    {
937        return m_data.getData();
938    }
939    
940    const DataTypes::ValueType&
941    DataExpanded::getVector() const
942    {
943        return m_data.getData();
944    }
945    
946  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.100  
changed lines
  Added in v.1800

  ViewVC Help
Powered by ViewVC 1.1.26