/[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 113 by jgs, Mon Feb 28 07:06:33 2005 UTC trunk/escript/src/DataExpanded.cpp revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 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>  /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #include "Data.h"
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 <iostream>  #include <boost/python/extract.hpp>
28    #include "DataMaths.h"
29    
30  using namespace std;  using namespace std;
31  using namespace boost::python;  using namespace boost::python;
32  using namespace boost;  using namespace boost;
33    using namespace escript::DataTypes;
34    
35    
36    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
37    
38    #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; *((int*)0)=17;throw DataException(ss.str());}
39    
40  namespace escript {  namespace escript {
41    
42  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const WrappedArray& value,
43                             const FunctionSpace& what)                             const FunctionSpace& what)
44    : DataAbstract(what)    : parent(what,value.getShape())
45  {  {
   DataArrayView::ShapeType tempShape;  
46    //    //
47    // extract the shape from the python array    // initialise the data array for this object
48    for (int i=0; i<value.getrank(); ++i) {    initialise(what.getNumSamples(),what.getNumDPPSample());
     //cout << extract<int>(value.getshape()[i]) << endl;  
     tempShape.push_back(extract<int>(value.getshape()[i]));  
   }  
   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  
49    //    //
50    // copy the given value to every data point    // copy the given value to every data point
51    copy(value);    copy(value);
52  }  }
53    
54  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
55    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(), other.getShape()),
56    m_data(other.m_data)    m_data(other.m_data)
57  {      {
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
58  }  }
59    
60  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
61    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
62  {      {
   initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());  
63    //    //
64    // DataConstant only has one value    // initialise the data array for this object
65    copy(other.getPointDataView());    initialise(other.getNumSamples(),other.getNumDPPSample());
66      //
67      // DataConstant only has one value, copy this to every data point
68      copy(other);
69  }  }
70    
71  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
72    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
73  {      {
74    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    //
75      // initialise the data array for this object
76      initialise(other.getNumSamples(),other.getNumDPPSample());
77      //
78      // for each data point in this object, extract and copy the corresponding data
79      // value from the given DataTagged object
80    int i,j;    int i,j;
81    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
82    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
83  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
84    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
85      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
86        try {        try {
87          getPointDataView().copy(getPointOffset(i,j), other.getPointDataView(), other.getPointOffset(i,j));             DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
88                                    other.getVectorRO(),
89                                    other.getPointOffset(i,j));
90        }        }
91        catch (std::exception& e) {        catch (std::exception& e) {
92          cout << e.what() << endl;          cout << e.what() << endl;
# Line 86  DataExpanded::DataExpanded(const DataTag Line 96  DataExpanded::DataExpanded(const DataTag
96  }  }
97    
98  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
99                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
100    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101  {  {
102    //    //
103    // get the shape of the slice    // get the shape of the slice
104    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
105    //    //
106    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
107    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
   //  
   DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);  
108    //    //
109    // copy the data    // copy the data
110    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
111    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
112      DataTypes::ValueType::size_type numCols=m_data.getNumCols();
113    int i,j;    int i,j;
114  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
115    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
116      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
117        try {        try {
118          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region_loop_range);          DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
119                                         other.getVectorRO(),
120                         other.getShape(),
121                                         other.getPointOffset(i,j),
122                                         region_loop_range);
123        }        }
124        catch (std::exception& e) {        catch (std::exception& e) {
125          cout << e.what() << endl;          cout << e.what() << endl;
# Line 115  DataExpanded::DataExpanded(const DataExp Line 128  DataExpanded::DataExpanded(const DataExp
128    }    }
129  }  }
130    
131  DataExpanded::DataExpanded(const DataArrayView& value,  DataExpanded::DataExpanded(const FunctionSpace& what,
132                             const FunctionSpace& what)                             const DataTypes::ShapeType &shape,
133    : DataAbstract(what)                             const DataTypes::ValueType &data)
134  {    : parent(what,shape)
135    DataArrayView::ShapeType tempShape=value.getShape();  {
136    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    EsysAssert(data.size()%getNoValues()==0,
137    copy(value);                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
138    
139      if (data.size()==getNoValues())
140      {
141         ValueType& vec=m_data.getData();
142         //
143         // create the view of the data
144         initialise(what.getNumSamples(),what.getNumDPPSample());
145         // now we copy this value to all elements
146         for (int i=0;i<getLength();)
147         {
148        for (unsigned int j=0;j<getNoValues();++j,++i)
149        {
150            vec[i]=data[j];
151        }
152         }
153      }
154      else
155      {
156         //
157         // copy the data in the correct format
158         m_data.getData()=data;
159      }
160    
161    
162  }  }
163    
164  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
165  {  {
166  }  }
167    
168  void  DataAbstract*
169  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataExpanded::deepCopy()
170  {  {
171    //    return new DataExpanded(*this);
   // 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);  
172  }  }
173    
174    
175  DataAbstract*  DataAbstract*
176  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataTypes::RegionType& region) const
177  {  {
178    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
179  }  }
180    
181  void  void
182  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
183                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
184  {  {
185    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
186    if (tempDataExp==0) {    if (tempDataExp==0) {
187      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
188    }    }
189      CHECK_FOR_EX_WRITE
190    //    //
191    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    // get shape of slice
192    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193      DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194    //    //
195    // check shape:    // check shape
196    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
197      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
198    }    }
199    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
200      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
201          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
202    }    }
203    //    //
204    // copy the data    // copy the data from the slice into this object
205    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
206    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
207    int i, j;    int i, j;
208  #pragma omp parallel for private(i,j) schedule(static)    ValueType& vec=getVectorRW();
209      const ShapeType& mshape=getShape();
210      const ValueType& tVec=tempDataExp->getVectorRO();
211      const ShapeType& tShape=tempDataExp->getShape();
212      #pragma omp parallel for private(i,j) schedule(static)
213    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
214      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
215        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region_loop_range);          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
216                                           tVec,
217                           tShape,
218                                           tempDataExp->getPointOffset(i,j),
219                                           region_loop_range);
220    
221      }      }
222    }    }
223  }  }
224    
225  void  void
226  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
227  {  {
228      EsysAssert((checkShape(getShape(), value.getShape())),
229                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
230    
231    //    //
232    // Copy a single value to every data point    // copy a single value to every data point in this object
233    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
234    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
235    int i, j;    int i,j;
236  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
237    for (i=0;i<nRows;++i) {    for (i=0;i<nRows;i++) {
238      for (j=0;j<nCols;++j) {      for (j=0;j<nCols;j++) {
239        //        DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
       // 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.  
       getPointDataView().copy(m_data.index(i,j),value);  
240      }      }
241    }    }
242  }  }
243    
244  void  void
245  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const WrappedArray& value)
246  {  {
247    //    // check the input shape matches this shape
248    // first convert the numarray into a DataArrayView format    if (!DataTypes::checkShape(getShape(),value.getShape())) {
249    DataArray temp(value);      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(  
250                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
251                          temp.getView().getShape()));                          value.getShape(),getShape()));
252    }    }
253    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the entire data structure  
   copy(temp.getView());  
254  }  }
255    
256    
257  void  void
258  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
259                           int noDataPointsPerSample)                           int noDataPointsPerSample)
260  {  {
261    //    if (noSamples==0)     //retain the default empty object
262    // resize data to the required size    {
263    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));       return;
264    //    }
265    // create a point data viewer of the data    //
266    DataArrayView temp(m_data.getData(),shape);    // resize data array to the required size
267    setPointDataView(temp);    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
268    }
269    
270    bool
271    DataExpanded::hasNaN() const
272    {
273        const ValueType& v=m_data.getData();
274        for (ValueType::size_type i=0;i<v.size();++i)
275        {
276            if (nancheck(v[i]))
277            {
278                return true;
279            }
280        }
281        return false;
282  }  }
283    
284    
285  string  string
286  DataExpanded::toString() const  DataExpanded::toString() const
287  {  {
288    stringstream temp;    stringstream temp;
289    //    FunctionSpace fs=getFunctionSpace();
290    // create a temporary view as the offset will be changed  
291    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    int offset=0;
292    for (int i=0;i<m_data.getNumRows();++i) {    for (int i=0;i<m_data.getNumRows();i++) {
293      for (int j=0;j<m_data.getNumCols();++j) {      for (int j=0;j<m_data.getNumCols();j++) {
294        tempView.setOffset(m_data.index(i,j));        offset=getPointOffset(i,j);
295        stringstream suffix;        stringstream suffix;
296        suffix << "(" << i << "," << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
297        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
298        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
299          temp << endl;          temp << endl;
300        }        }
301      }      }
302    }    }
303      string result=temp.str();
304      if (result.empty())
305      {
306        return "(data contains no samples)\n";
307      }
308    return temp.str();    return temp.str();
309  }  }
310    
311  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
312  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
313                               int dataPointNo) const                               int dataPointNo) const
314  {  {
315    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
316  }  }
317    
318  DataArrayView  DataTypes::ValueType::size_type
319  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
320                             int dataPointNo)                               int dataPointNo)
321  {  {
322    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));    return m_data.index(sampleNo,dataPointNo);
   return temp;  
323  }  }
324    
325  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
326  DataExpanded::getLength() const  DataExpanded::getLength() const
327  {  {
328    return m_data.getData().size();    return m_data.size();
329  }  }
330    
331    
332    
333  void  void
334  DataExpanded::setRefValue(int ref,  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
335                            const DataArray& value)    CHECK_FOR_EX_WRITE
 {  
336    //    //
337    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
338    int numSamples = getNumSamples();    int numSamples = getNumSamples();
339    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
340      int dataPointRank = getRank();
341      ShapeType dataPointShape = getShape();
342      if (numSamples*numDataPointsPerSample > 0) {
343         //TODO: global error handling
344         if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
345              throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
346         }
347         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
348               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
349         }
350         ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
351         ValueType& vec=getVectorRW();
352         if (dataPointRank==0) {
353             vec[offset]=value;
354         } else if (dataPointRank==1) {
355            for (int i=0; i<dataPointShape[0]; i++) {
356                vec[offset+i]=value;
357            }
358         } else if (dataPointRank==2) {
359            for (int i=0; i<dataPointShape[0]; i++) {
360               for (int j=0; j<dataPointShape[1]; j++) {
361                  vec[offset+getRelIndex(dataPointShape,i,j)]=value;
362               }
363            }
364         } else if (dataPointRank==3) {
365            for (int i=0; i<dataPointShape[0]; i++) {
366               for (int j=0; j<dataPointShape[1]; j++) {
367                  for (int k=0; k<dataPointShape[2]; k++) {
368                     vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
369                  }
370               }
371            }
372         } else if (dataPointRank==4) {
373             for (int i=0; i<dataPointShape[0]; i++) {
374               for (int j=0; j<dataPointShape[1]; j++) {
375                 for (int k=0; k<dataPointShape[2]; k++) {
376                   for (int l=0; l<dataPointShape[3]; l++) {
377                      vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
378                   }
379                 }
380               }
381             }
382         }
383      }
384    }
385    
386    void
387    DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
388      CHECK_FOR_EX_WRITE
389      //
390      // Get the number of samples and data-points per sample.
391      int numSamples = getNumSamples();
392      int numDataPointsPerSample = getNumDPPSample();
393    //    //
394    // Determine the sample number which corresponds to this reference number.    // check rank:
395    int sampleNo = -1;    if (value.getRank()!=getRank())
396    int tempRef = -1;         throw DataException("Rank of value does not match Data object rank");
397    for (int n=0; n<numSamples; n++) {    if (numSamples*numDataPointsPerSample > 0) {
398      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);       //TODO: global error handling
399      if (tempRef == ref) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
400        sampleNo = n;            throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
401        break;       }
402         if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
403               throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
404         }
405         ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
406         ValueType& vec=getVectorRW();
407         vec.copyFromArrayToOffset(value,offset,1);
408      }
409    }
410    
411    void
412    DataExpanded::symmetric(DataAbstract* ev)
413    {
414      int sampleNo,dataPointNo;
415      int numSamples = getNumSamples();
416      int numDataPointsPerSample = getNumDPPSample();
417      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
418      if (temp_ev==0) {
419        throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
420      }
421      const ValueType& vec=getVectorRO();
422      const ShapeType& shape=getShape();
423      ValueType& evVec=temp_ev->getVectorRW();
424      const ShapeType& evShape=temp_ev->getShape();
425      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
426      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
427        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
428             DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
429                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
430      }      }
431    }    }
432    if (sampleNo == -1) {  }
433      throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");  void
434    DataExpanded::nonsymmetric(DataAbstract* ev)
435    {
436      int sampleNo,dataPointNo;
437      int numSamples = getNumSamples();
438      int numDataPointsPerSample = getNumDPPSample();
439      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
440      if (temp_ev==0) {
441        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
442      }
443      const ValueType& vec=getVectorRO();
444      const ShapeType& shape=getShape();
445      ValueType& evVec=temp_ev->getVectorRW();
446      const ShapeType& evShape=temp_ev->getShape();
447      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
448      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
449        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
450             DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
451                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
452        }
453    }    }
454    }
455    void
456    DataExpanded::trace(DataAbstract* ev, int axis_offset)
457    {
458      int sampleNo,dataPointNo;
459      int numSamples = getNumSamples();
460      int numDataPointsPerSample = getNumDPPSample();
461      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
462      if (temp_ev==0) {
463        throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
464      }
465      const ValueType& vec=getVectorRO();
466      const ShapeType& shape=getShape();
467      ValueType& evVec=temp_ev->getVectorRW();
468      const ShapeType& evShape=temp_ev->getShape();
469      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
470      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
471        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
472             DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
473                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
474        }
475      }
476    }
477    
478    for (int n=0; n<numDPPSample; n++) {  void
479      //  DataExpanded::transpose(DataAbstract* ev, int axis_offset)
480      // Get each data-point in the sample in turn.  {
481      DataArrayView pointView = getDataPoint(sampleNo, n);    int sampleNo,dataPointNo;
482      int numSamples = getNumSamples();
483      //    int numDataPointsPerSample = getNumDPPSample();
484      // Assign the values in the DataArray to this data-point.    DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
485      pointView.copy(value.getView());    if (temp_ev==0) {
486        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
487      }
488      const ValueType& vec=getVectorRO();
489      const ShapeType& shape=getShape();
490      ValueType& evVec=temp_ev->getVectorRW();
491      const ShapeType& evShape=temp_ev->getShape();
492      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
493      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
494        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
495             DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
496                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
497        }
498    }    }
499  }  }
500    
501  void  void
502  DataExpanded::getRefValue(int ref,  DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
                           DataArray& value)  
503  {  {
504    //    int sampleNo,dataPointNo;
505    // Get the number of samples and data-points per sample.    int numSamples = getNumSamples();
506      int numDataPointsPerSample = getNumDPPSample();
507      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
508      if (temp_ev==0) {
509        throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
510      }
511      const ValueType& vec=getVectorRO();
512      const ShapeType& shape=getShape();
513      ValueType& evVec=temp_ev->getVectorRW();
514      const ShapeType& evShape=temp_ev->getShape();
515      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
516      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
517        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
518             DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
519                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
520        }
521      }
522    }
523    void
524    DataExpanded::eigenvalues(DataAbstract* ev)
525    {
526      int sampleNo,dataPointNo;
527    int numSamples = getNumSamples();    int numSamples = getNumSamples();
528    int numDPPSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
529      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
530      if (temp_ev==0) {
531        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
532      }
533      const ValueType& vec=getVectorRO();
534      const ShapeType& shape=getShape();
535      ValueType& evVec=temp_ev->getVectorRW();
536      const ShapeType& evShape=temp_ev->getShape();
537      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
539        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
540             DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
541                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
542        }
543      }
544    }
545    void
546    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
547    {
548      int numSamples = getNumSamples();
549      int numDataPointsPerSample = getNumDPPSample();
550      int sampleNo,dataPointNo;
551      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
552      if (temp_ev==0) {
553        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
554      }
555      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
556      if (temp_V==0) {
557        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
558      }
559      const ValueType& vec=getVectorRO();
560      const ShapeType& shape=getShape();
561      ValueType& evVec=temp_ev->getVectorRW();
562      const ShapeType& evShape=temp_ev->getShape();
563      ValueType& VVec=temp_V->getVectorRW();
564      const ShapeType& VShape=temp_V->getShape();
565      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
566      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
567        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
568             DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
569                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
570                                        VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
571        }
572      }
573    }
574    
575    //  
576    // Determine the sample number which corresponds to this reference number  int
577    int sampleNo = -1;  DataExpanded::matrixInverse(DataAbstract* out) const
578    int tempRef = -1;  {
579    for (int n=0; n<numSamples; n++) {    DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
580      tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);    if (temp==0)
581      if (tempRef == ref) {    {
582        sampleNo = n;      throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
583        break;    }
584    
585      if (getRank()!=2)
586      {
587        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
588    
589      }
590      int  sampleNo;
591      const int numdpps=getNumDPPSample();
592      const int numSamples = getNumSamples();
593      const ValueType& vec=m_data.getData();
594      int errcode=0;
595      #pragma omp parallel private(sampleNo)
596      {
597         int errorcode=0;
598         LapackInverseHelper h(getShape()[0]);
599         #pragma omp for schedule(static)
600         for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
601         {
602                // not sure I like all those virtual calls to getPointOffset
603            DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
604            int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
605        if (res>errorcode)
606        {
607            errorcode=res;
608            #pragma omp critical
609            {
610              errcode=errorcode;    // I'm not especially concerned which error gets reported as long as one is
611            }
612        }
613         }
614      }
615      return errcode;
616      if (errcode)
617      {
618        DataMaths::matrixInverseError(errcode); // throws exceptions
619      }
620    }
621    
622    void
623    DataExpanded::setToZero(){
624    // TODO: Surely there is a more efficient way to do this????
625    // Why is there no memset here? Parallel issues?
626    // A: This ensures that memory is touched by the correct thread.
627      CHECK_FOR_EX_WRITE
628      int numSamples = getNumSamples();
629      int numDataPointsPerSample = getNumDPPSample();
630      DataTypes::ValueType::size_type n = getNoValues();
631      double* p;
632      int  sampleNo,dataPointNo, i;
633      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
634      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
635        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
636            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
637            for (i=0; i<n ;++i) p[i]=0.;
638      }      }
639    }    }
640    if (sampleNo == -1) {  }
641      throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");  
642    /* Append MPI rank to file name if multiple MPI processes */
643    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
644      /* Make plenty of room for the mpi_rank number and terminating '\0' */
645      char *newFileName = (char *)malloc(strlen(fileName)+20);
646      strncpy(newFileName, fileName, strlen(fileName)+1);
647      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
648      return(newFileName);
649    }
650    
651    void
652    DataExpanded::dump(const std::string fileName) const
653    {
654       #ifdef USE_NETCDF
655       const int ldims=2+DataTypes::maxRank;
656       const NcDim* ncdims[ldims];
657       NcVar *var, *ids;
658       int rank = getRank();
659       int type=  getFunctionSpace().getTypeCode();
660       int ndims =0;
661       long dims[ldims];
662       const double* d_ptr=&(m_data[0]);
663       const DataTypes::ShapeType& shape = getShape();
664       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
665       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
666    #ifdef PASO_MPI
667       MPI_Status status;
668    #endif
669    
670    #ifdef PASO_MPI
671       /* Serialize NetCDF I/O */
672       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
673    #endif
674    
675       // netCDF error handler
676       NcError err(NcError::verbose_nonfatal);
677       // Create the file.
678       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
679       NcFile dataFile(newFileName, NcFile::Replace);
680       // check if writing was successful
681       if (!dataFile.is_valid())
682            throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
683       if (!dataFile.add_att("type_id",2) )
684            throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
685       if (!dataFile.add_att("rank",rank) )
686            throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
687       if (!dataFile.add_att("function_space_type",type))
688            throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
689       ndims=rank+2;
690       if ( rank >0 ) {
691           dims[0]=shape[0];
692           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
693                throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
694       }
695       if ( rank >1 ) {
696           dims[1]=shape[1];
697           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
698                throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
699       }
700       if ( rank >2 ) {
701           dims[2]=shape[2];
702           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
703                throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
704       }
705       if ( rank >3 ) {
706           dims[3]=shape[3];
707           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
708                throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
709       }
710       dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
711       if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
712                throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
713       dims[rank+1]=getFunctionSpace().getNumSamples();
714       if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
715                throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
716    
717       if (getFunctionSpace().getNumSamples()>0)
718       {
719    
720         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
721            throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
722         const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
723         if (! (ids->put(ids_p,dims[rank+1])) )
724            throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
725         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
726            throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
727         if (! (var->put(d_ptr,dims)) )
728            throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
729       }
730    #ifdef PASO_MPI
731       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
732    #endif
733       #else
734       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
735       #endif
736    }
737    
738    void  
739    DataExpanded::setTaggedValue(int tagKey,
740               const DataTypes::ShapeType& pointshape,
741                   const DataTypes::ValueType& value,
742               int dataOffset)
743    {
744      CHECK_FOR_EX_WRITE
745      int numSamples = getNumSamples();
746      int numDataPointsPerSample = getNumDPPSample();
747      int sampleNo,dataPointNo, i;
748      DataTypes::ValueType::size_type n = getNoValues();
749      double* p;
750      const double* in=&value[0+dataOffset];
751      
752      if (value.size() != n) {
753        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
754      }
755    
756      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
757      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
758        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
759            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
760                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
761                for (i=0; i<n ;++i) p[i]=in[i];
762            }
763        }
764    }    }
765    }
766    
   //  
   // Get the first data-point associated with this sample number.  
   DataArrayView pointView = getDataPoint(sampleNo, 0);  
767    
768    //  void
769    // Load the values from this data-point into the DataArray  DataExpanded::reorderByReferenceIDs(int *reference_ids)
770    value.getView().copy(pointView);  {
771      CHECK_FOR_EX_WRITE
772      int numSamples = getNumSamples();
773      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
774      int sampleNo, sampleNo2,i;
775      double* p,*p2;
776      register double rtmp;
777      FunctionSpace fs=getFunctionSpace();
778    
779      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
780         const int id_in=reference_ids[sampleNo];
781         const int id=fs.getReferenceIDOfSample(sampleNo);
782         if (id!=id_in) {
783             bool matched=false;
784             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
785                  if (id == reference_ids[sampleNo2]) {
786                     p=&(m_data[getPointOffset(sampleNo,0)]);
787                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
788                     for (i=0; i<n ;i++) {
789                             rtmp=p[i];
790                             p[i]=p2[i];
791                             p2[i]=rtmp;
792                     }
793                     reference_ids[sampleNo]=id;
794                     reference_ids[sampleNo2]=id_in;
795                     matched=true;
796                     break;
797                  }
798             }
799             if (! matched) {
800                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
801             }
802         }
803       }
804    }
805    
806    DataTypes::ValueType&
807    DataExpanded::getVectorRW()
808    {
809        CHECK_FOR_EX_WRITE
810        return m_data.getData();
811  }  }
812    
813    const DataTypes::ValueType&
814    DataExpanded::getVectorRO() const
815    {
816        return m_data.getData();
817    }
818    
819    
820  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.113  
changed lines
  Added in v.2881

  ViewVC Help
Powered by ViewVC 1.1.26