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

Legend:
Removed from v.97  
changed lines
  Added in v.1796

  ViewVC Help
Powered by ViewVC 1.1.26