/[escript]/trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1827 by ksteube, Thu Oct 2 04:28:07 2008 UTC revision 3468 by jfenwick, Tue Feb 22 06:38:57 2011 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 17  Line 17 
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
19  #include "DataTagged.h"  #include "DataTagged.h"
20    
21    #include "esysUtils/Esys_MPI.h"
22    
23  #ifdef USE_NETCDF  #ifdef USE_NETCDF
24  #include <netcdfcpp.h>  #include <netcdfcpp.h>
25  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28  #include "DataMaths.h"  #include "DataMaths.h"
# Line 32  using namespace boost::python; Line 32  using namespace boost::python;
32  using namespace boost;  using namespace boost;
33  using namespace escript::DataTypes;  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,DataTypes::shapeFromNumArray(value))    : parent(what,value.getShape())
45  {  {
46    //    //
47    // initialise the data array for this object    // initialise the data array for this object
# Line 47  DataExpanded::DataExpanded(const boost:: Line 52  DataExpanded::DataExpanded(const boost::
52  }  }
53    
54  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
55    : DataAbstract(other.getFunctionSpace(), other.getShape()),    : parent(other.getFunctionSpace(), other.getShape()),
56    m_data(other.m_data)    m_data(other.m_data)
57  {  {
58  }  }
59    
60  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
61    : DataAbstract(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
62  {  {
63    //    //
64    // initialise the data array for this object    // initialise the data array for this object
# Line 64  DataExpanded::DataExpanded(const DataCon Line 69  DataExpanded::DataExpanded(const DataCon
69  }  }
70    
71  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
72    : DataAbstract(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
73  {  {
74    //    //
75    // initialise the data array for this object    // initialise the data array for this object
# Line 79  DataExpanded::DataExpanded(const DataTag Line 84  DataExpanded::DataExpanded(const DataTag
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             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),             DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
88                                  other.getVector(),                                  other.getVectorRO(),
89                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
90        }        }
91        catch (std::exception& e) {        catch (std::exception& e) {
# Line 92  DataExpanded::DataExpanded(const DataTag Line 97  DataExpanded::DataExpanded(const DataTag
97    
98  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
99                             const DataTypes::RegionType& region)                             const DataTypes::RegionType& region)
100    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101  {  {
102    //    //
103    // get the shape of the slice    // get the shape of the slice
# Line 110  DataExpanded::DataExpanded(const DataExp Line 115  DataExpanded::DataExpanded(const DataExp
115    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
116      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
117        try {        try {
118  //         getPointDataView().copySlice(getPointOffset(i,j),          DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
119  //                                      other.getPointDataView(),                                       other.getVectorRO(),
 //                                      other.getPointOffset(i,j),  
 //                                      region_loop_range);  
         DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),  
                                      other.getVector(),  
120                       other.getShape(),                       other.getShape(),
121                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
122                                       region_loop_range);                                       region_loop_range);
# Line 127  DataExpanded::DataExpanded(const DataExp Line 128  DataExpanded::DataExpanded(const DataExp
128    }    }
129  }  }
130    
 // DataExpanded::DataExpanded(const DataArrayView& value,  
 //                            const FunctionSpace& what)  
 //   : DataAbstract(what)  
 // {  
 //   //  
 //   // get the shape of the given data value  
 //   DataTypes::ShapeType tempShape=value.getShape();  
 //   //  
 //   // initialise this Data object to the shape of the given data value  
 //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  
 //   //  
 //   // copy the given value to every data point  
 //   copy(value);  
 // }  
   
131  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
132                             const DataTypes::ShapeType &shape,                             const DataTypes::ShapeType &shape,
133                             const DataTypes::ValueType &data)                             const DataTypes::ValueType &data)
134    : DataAbstract(what,shape)    : parent(what,shape)
135  {  {
136    EsysAssert(data.size()%getNoValues()==0,    EsysAssert(data.size()%getNoValues()==0,
137                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
# Line 159  DataExpanded::DataExpanded(const Functio Line 145  DataExpanded::DataExpanded(const Functio
145       // now we copy this value to all elements       // now we copy this value to all elements
146       for (int i=0;i<getLength();)       for (int i=0;i<getLength();)
147       {       {
148      for (int j=0;j<getNoValues();++j,++i)      for (unsigned int j=0;j<getNoValues();++j,++i)
149      {      {
150          vec[i]=data[j];          vec[i]=data[j];
151      }      }
# Line 175  DataExpanded::DataExpanded(const Functio Line 161  DataExpanded::DataExpanded(const Functio
161    
162  }  }
163    
164    DataExpanded::DataExpanded(const FunctionSpace& what,
165                               const DataTypes::ShapeType &shape,
166                               const double v)
167      : parent(what,shape)
168    {
169         ValueType& vec=m_data.getData();
170         //
171         // create the view of the data
172         initialise(what.getNumSamples(),what.getNumDPPSample());
173         // now we copy this value to all elements
174         const int L=getLength();
175         int i;
176         #pragma omp parallel for schedule(static) private(i)
177         for (i=0;i<L;++i)
178         {
179            vec[i]=v;
180         }
181    }
182    
183    
184    
185  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
186  {  {
187  }  }
# Line 200  DataExpanded::setSlice(const DataAbstrac Line 207  DataExpanded::setSlice(const DataAbstrac
207    if (tempDataExp==0) {    if (tempDataExp==0) {
208      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
209    }    }
210      CHECK_FOR_EX_WRITE
211    //    //
212    // get shape of slice    // get shape of slice
213    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
# Line 218  DataExpanded::setSlice(const DataAbstrac Line 226  DataExpanded::setSlice(const DataAbstrac
226    DataTypes::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
227    DataTypes::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
228    int i, j;    int i, j;
229    ValueType& vec=getVector();    ValueType& vec=getVectorRW();
230    const ShapeType& mshape=getShape();    const ShapeType& mshape=getShape();
231    const ValueType& tVec=tempDataExp->getVector();    const ValueType& tVec=tempDataExp->getVectorRO();
232    const ShapeType& tShape=tempDataExp->getShape();    const ShapeType& tShape=tempDataExp->getShape();
233    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
234    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
235      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
 /*      getPointDataView().copySliceFrom(getPointOffset(i,j),  
                                        tempDataExp->getPointDataView(),  
                                        tempDataExp->getPointOffset(i,j),  
                                        region_loop_range);*/  
236          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
237                                         tVec,                                         tVec,
238                         tShape,                         tShape,
# Line 253  DataExpanded::copy(const DataConstant& v Line 257  DataExpanded::copy(const DataConstant& v
257    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
258    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
259      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
260        // NOTE: An exception may be thown from this call if        DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
       // DOASSERT is on which of course will play  
       // havoc with the omp threads. Run single threaded  
       // if using DOASSERT.  
       //getPointDataView().copy(getPointOffset(i,j),value);  
       DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);  
261      }      }
262    }    }
263  }  }
264    
265    void
266  // void  DataExpanded::copy(const WrappedArray& value)
267  // DataExpanded::copy(const DataArrayView& value)  {
 // {  
 //   //  
 //   // copy a single value to every data point in this object  
 //   int nRows=m_data.getNumRows();  
 //   int nCols=m_data.getNumCols();  
 //   int i,j;  
 //   #pragma omp parallel for private(i,j) schedule(static)  
 //   for (i=0;i<nRows;i++) {  
 //     for (j=0;j<nCols;j++) {  
 //       // 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(getPointOffset(i,j),value);  
 //     }  
 //   }  
 // }  
   
 void  
 DataExpanded::copy(const boost::python::numeric::array& value)  
 {  
   
   // extract the shape of the numarray  
   DataTypes::ShapeType tempShape;  
   for (int i=0; i < value.getrank(); i++) {  
     tempShape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 //   DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(value);  
   
   //  
268    // check the input shape matches this shape    // check the input shape matches this shape
269    if (!DataTypes::checkShape(getShape(),tempShape)) {    if (!DataTypes::checkShape(getShape(),value.getShape())) {
270      throw DataException(DataTypes::createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
271                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
272                          tempShape,getShape()));                          value.getShape(),getShape()));
273    }    }
274    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the data  
   //copy(temp_dataView);  
   getVector().copyFromNumArray(value);  
275  }  }
276    
277    
# Line 327  DataExpanded::initialise(int noSamples, Line 288  DataExpanded::initialise(int noSamples,
288    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
289  }  }
290    
291    bool
292    DataExpanded::hasNaN() const
293    {
294        const ValueType& v=m_data.getData();
295        for (ValueType::size_type i=0;i<v.size();++i)
296        {
297            if (nancheck(v[i]))
298            {
299                return true;
300            }
301        }
302        return false;
303    }
304    
305    
306  string  string
307  DataExpanded::toString() const  DataExpanded::toString() const
308  {  {
# Line 339  DataExpanded::toString() const Line 315  DataExpanded::toString() const
315        offset=getPointOffset(i,j);        offset=getPointOffset(i,j);
316        stringstream suffix;        stringstream suffix;
317        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
318        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());        temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
319        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
320          temp << endl;          temp << endl;
321        }        }
# Line 361  DataExpanded::getPointOffset(int sampleN Line 337  DataExpanded::getPointOffset(int sampleN
337  }  }
338    
339  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
340    DataExpanded::getPointOffset(int sampleNo,
341                                 int dataPointNo)
342    {
343      return m_data.index(sampleNo,dataPointNo);
344    }
345    
346    DataTypes::ValueType::size_type
347  DataExpanded::getLength() const  DataExpanded::getLength() const
348  {  {
349    return m_data.size();    return m_data.size();
# Line 370  DataExpanded::getLength() const Line 353  DataExpanded::getLength() const
353    
354  void  void
355  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
356      CHECK_FOR_EX_WRITE
357    //    //
358    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
359    int numSamples = getNumSamples();    int numSamples = getNumSamples();
# Line 385  DataExpanded::copyToDataPoint(const int Line 369  DataExpanded::copyToDataPoint(const int
369             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
370       }       }
371       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
372       ValueType& vec=getVector();       ValueType& vec=getVectorRW();
373       if (dataPointRank==0) {       if (dataPointRank==0) {
374           vec[0]=value;           vec[offset]=value;
375       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
376          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
377              vec[offset+i]=value;              vec[offset+i]=value;
# Line 419  DataExpanded::copyToDataPoint(const int Line 403  DataExpanded::copyToDataPoint(const int
403       }       }
404    }    }
405  }  }
406    
407  void  void
408  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
409      CHECK_FOR_EX_WRITE
410    //    //
411    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
412    int numSamples = getNumSamples();    int numSamples = getNumSamples();
413    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
   int dataPointRank = getRank();  
   const ShapeType& shape = getShape();  
414    //    //
415    // check rank:    // check rank:
416    if (value.getrank()!=dataPointRank)    if (value.getRank()!=getRank())
417         throw DataException("Rank of numarray does not match Data object rank");         throw DataException("Rank of value does not match Data object rank");
418    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
419       //TODO: global error handling       //TODO: global error handling
420       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 439  DataExpanded::copyToDataPoint(const int Line 423  DataExpanded::copyToDataPoint(const int
423       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
424             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
425       }       }
426       ValueType& vec=getVector();       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
427       if (dataPointRank==0) {       ValueType& vec=getVectorRW();
428           vec[0]=extract<double>(value[0]);       vec.copyFromArrayToOffset(value,offset,1);
      } else if (dataPointRank==1) {  
         for (int i=0; i<shape[0]; i++) {  
             vec[i]=extract<double>(value[i]);  
         }  
      } else if (dataPointRank==2) {  
         for (int i=0; i<shape[0]; i++) {  
            for (int j=0; j<shape[1]; j++) {  
               vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);  
            }  
         }  
      } else if (dataPointRank==3) {  
         for (int i=0; i<shape[0]; i++) {  
            for (int j=0; j<shape[1]; j++) {  
               for (int k=0; k<shape[2]; k++) {  
                  vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);  
               }  
            }  
         }  
      } else if (dataPointRank==4) {  
          for (int i=0; i<shape[0]; i++) {  
            for (int j=0; j<shape[1]; j++) {  
              for (int k=0; k<shape[2]; k++) {  
                for (int l=0; l<shape[3]; l++) {  
                   vec[getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);  
                }  
              }  
            }  
          }  
      }  
   }  
 }  
 void  
 DataExpanded::copyAll(const boost::python::numeric::array& value) {  
   //  
   // Get the number of samples and data-points per sample.  
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDPPSample();  
   int dataPointRank = getRank();  
   const ShapeType& dataPointShape = getShape();  
   //  
   // check rank:  
   if (value.getrank()!=dataPointRank+1)  
        throw DataException("Rank of numarray does not match Data object rank");  
   if (value.getshape()[0]!=numSamples*numDataPointsPerSample)  
        throw DataException("leading dimension of numarray is too small");  
   //  
   ValueType& vec=getVector();  
   int dataPoint = 0;  
   for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {  
     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {  
       ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);  
       if (dataPointRank==0) {  
          vec[offset]=extract<double>(value[dataPoint]);  
       } else if (dataPointRank==1) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
             vec[offset+i]=extract<double>(value[dataPoint][i]);  
          }  
       } else if (dataPointRank==2) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
          vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);  
            }  
          }  
        } else if (dataPointRank==3) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
              for (int k=0; k<dataPointShape[2]; k++) {  
          vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);  
              }  
            }  
          }  
        } else if (dataPointRank==4) {  
          for (int i=0; i<dataPointShape[0]; i++) {  
            for (int j=0; j<dataPointShape[1]; j++) {  
              for (int k=0; k<dataPointShape[2]; k++) {  
                for (int l=0; l<dataPointShape[3]; l++) {  
                  vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);  
                }  
              }  
            }  
          }  
       }  
       dataPoint++;  
     }  
429    }    }
430  }  }
431    
432  void  void
433  DataExpanded::symmetric(DataAbstract* ev)  DataExpanded::symmetric(DataAbstract* ev)
434  {  {
# Line 538  DataExpanded::symmetric(DataAbstract* ev Line 439  DataExpanded::symmetric(DataAbstract* ev
439    if (temp_ev==0) {    if (temp_ev==0) {
440      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
441    }    }
442    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
443    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
444    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
445    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
446    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
447    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 560  DataExpanded::nonsymmetric(DataAbstract* Line 461  DataExpanded::nonsymmetric(DataAbstract*
461    if (temp_ev==0) {    if (temp_ev==0) {
462      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
463    }    }
464    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
465    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
466    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
467    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
468    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
469    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 582  DataExpanded::trace(DataAbstract* ev, in Line 483  DataExpanded::trace(DataAbstract* ev, in
483    if (temp_ev==0) {    if (temp_ev==0) {
484      throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
485    }    }
486    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
487    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
488    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
489    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
490    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
491    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 605  DataExpanded::transpose(DataAbstract* ev Line 506  DataExpanded::transpose(DataAbstract* ev
506    if (temp_ev==0) {    if (temp_ev==0) {
507      throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
508    }    }
509    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
510    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
511    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
512    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
513    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
514    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 628  DataExpanded::swapaxes(DataAbstract* ev, Line 529  DataExpanded::swapaxes(DataAbstract* ev,
529    if (temp_ev==0) {    if (temp_ev==0) {
530      throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
531    }    }
532    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
533    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
534    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
535    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
536    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
537    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 650  DataExpanded::eigenvalues(DataAbstract* Line 551  DataExpanded::eigenvalues(DataAbstract*
551    if (temp_ev==0) {    if (temp_ev==0) {
552      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
553    }    }
554    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
555    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
556    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
557    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
558    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
559    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 676  DataExpanded::eigenvalues_and_eigenvecto Line 577  DataExpanded::eigenvalues_and_eigenvecto
577    if (temp_V==0) {    if (temp_V==0) {
578      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
579    }    }
580    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
581    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
582    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
583    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
584    ValueType& VVec=temp_V->getVector();    ValueType& VVec=temp_V->getVectorRW();
585    const ShapeType& VShape=temp_V->getShape();    const ShapeType& VShape=temp_V->getShape();
586    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
587    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 692  DataExpanded::eigenvalues_and_eigenvecto Line 593  DataExpanded::eigenvalues_and_eigenvecto
593    }    }
594  }  }
595    
596    
597    int
598    DataExpanded::matrixInverse(DataAbstract* out) const
599    {
600      DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
601      if (temp==0)
602      {
603        throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
604      }
605    
606      if (getRank()!=2)
607      {
608        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
609    
610      }
611      int  sampleNo;
612      const int numdpps=getNumDPPSample();
613      const int numSamples = getNumSamples();
614      const ValueType& vec=m_data.getData();
615      int errcode=0;
616      #pragma omp parallel private(sampleNo)
617      {
618         int errorcode=0;
619         LapackInverseHelper h(getShape()[0]);
620         #pragma omp for schedule(static)
621         for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
622         {
623                // not sure I like all those virtual calls to getPointOffset
624            DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
625            int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
626        if (res>errorcode)
627        {
628            errorcode=res;
629            #pragma omp critical
630            {
631              errcode=errorcode;    // I'm not especially concerned which error gets reported as long as one is
632            }
633        }
634         }
635      }
636      return errcode;
637      if (errcode)
638      {
639        DataMaths::matrixInverseError(errcode); // throws exceptions
640      }
641    }
642    
643  void  void
644  DataExpanded::setToZero(){  DataExpanded::setToZero(){
645  // TODO: Surely there is a more efficient way to do this????  // TODO: Surely there is a more efficient way to do this????
646  // Why is there no memset here? Parallel issues?  // Why is there no memset here? Parallel issues?
647    // A: This ensures that memory is touched by the correct thread.
648      CHECK_FOR_EX_WRITE
649    int numSamples = getNumSamples();    int numSamples = getNumSamples();
650    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
651    DataTypes::ValueType::size_type n = getNoValues();    DataTypes::ValueType::size_type n = getNoValues();
# Line 732  DataExpanded::dump(const std::string fil Line 682  DataExpanded::dump(const std::string fil
682     long dims[ldims];     long dims[ldims];
683     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
684     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
685     int mpi_iam=getFunctionSpace().getDomain().getMPIRank();     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
686     int mpi_num=getFunctionSpace().getDomain().getMPISize();     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
687  #ifdef PASO_MPI  #ifdef ESYS_MPI
688     MPI_Status status;     MPI_Status status;
689  #endif  #endif
690    
691  #ifdef PASO_MPI  #ifdef ESYS_MPI
692     /* Serialize NetCDF I/O */     /* Serialize NetCDF I/O */
693     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
694  #endif  #endif
# Line 798  DataExpanded::dump(const std::string fil Line 748  DataExpanded::dump(const std::string fil
748       if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
749          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
750     }     }
751  #ifdef PASO_MPI  #ifdef ESYS_MPI
752     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
753  #endif  #endif
754     #else     #else
# Line 812  DataExpanded::setTaggedValue(int tagKey, Line 762  DataExpanded::setTaggedValue(int tagKey,
762                 const DataTypes::ValueType& value,                 const DataTypes::ValueType& value,
763             int dataOffset)             int dataOffset)
764  {  {
765      CHECK_FOR_EX_WRITE
766    int numSamples = getNumSamples();    int numSamples = getNumSamples();
767    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
768    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
# Line 838  DataExpanded::setTaggedValue(int tagKey, Line 789  DataExpanded::setTaggedValue(int tagKey,
789  void  void
790  DataExpanded::reorderByReferenceIDs(int *reference_ids)  DataExpanded::reorderByReferenceIDs(int *reference_ids)
791  {  {
792      CHECK_FOR_EX_WRITE
793    int numSamples = getNumSamples();    int numSamples = getNumSamples();
794    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
795    int sampleNo, sampleNo2,i;    int sampleNo, sampleNo2,i;
# Line 873  DataExpanded::reorderByReferenceIDs(int Line 825  DataExpanded::reorderByReferenceIDs(int
825  }  }
826    
827  DataTypes::ValueType&  DataTypes::ValueType&
828  DataExpanded::getVector()  DataExpanded::getVectorRW()
829  {  {
830        CHECK_FOR_EX_WRITE
831      return m_data.getData();      return m_data.getData();
832  }  }
833    
834  const DataTypes::ValueType&  const DataTypes::ValueType&
835  DataExpanded::getVector() const  DataExpanded::getVectorRO() const
836  {  {
837      return m_data.getData();      return m_data.getData();
838  }  }
839    
840    void DataExpanded::randomFill(double seed)
841    {
842        CHECK_FOR_EX_WRITE
843        if (seed==0)
844        {
845        time_t s=time(0);
846        seed=s;
847        }
848        srand(seed);
849        DataVector&  dv=getVectorRW();
850        for (long i=0;i<dv.size();++i)
851        {
852        dv[i]=(double)rand()/RAND_MAX;
853        }
854    }
855    
856  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1827  
changed lines
  Added in v.3468

  ViewVC Help
Powered by ViewVC 1.1.26