/[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 3506 by jfenwick, Wed May 11 01:59:45 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    #include <limits>
21    
22    #include "esysUtils/Esys_MPI.h"
23    
24  #ifdef USE_NETCDF  #ifdef USE_NETCDF
25  #include <netcdfcpp.h>  #include <netcdfcpp.h>
26  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
27    
28  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
29  #include "DataMaths.h"  #include "DataMaths.h"
30    
31    //#define MKLRANDOM
32    
33    #ifdef MKLRANDOM
34    #include <mkl_vsl.h>
35    
36    #endif
37    
38  using namespace std;  using namespace std;
39  using namespace boost::python;  using namespace boost::python;
40  using namespace boost;  using namespace boost;
41  using namespace escript::DataTypes;  using namespace escript::DataTypes;
42    
43    
44    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
45    
46    #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());}
47    
48  namespace escript {  namespace escript {
49    
50  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const WrappedArray& value,
51                             const FunctionSpace& what)                             const FunctionSpace& what)
52    : DataAbstract(what,DataTypes::shapeFromNumArray(value))    : parent(what,value.getShape())
53  {  {
54    //    //
55    // initialise the data array for this object    // initialise the data array for this object
# Line 47  DataExpanded::DataExpanded(const boost:: Line 60  DataExpanded::DataExpanded(const boost::
60  }  }
61    
62  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
63    : DataAbstract(other.getFunctionSpace(), other.getShape()),    : parent(other.getFunctionSpace(), other.getShape()),
64    m_data(other.m_data)    m_data(other.m_data)
65  {  {
66  }  }
67    
68  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
69    : DataAbstract(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
70  {  {
71    //    //
72    // initialise the data array for this object    // initialise the data array for this object
# Line 64  DataExpanded::DataExpanded(const DataCon Line 77  DataExpanded::DataExpanded(const DataCon
77  }  }
78    
79  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
80    : DataAbstract(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
81  {  {
82    //    //
83    // initialise the data array for this object    // initialise the data array for this object
# Line 79  DataExpanded::DataExpanded(const DataTag Line 92  DataExpanded::DataExpanded(const DataTag
92    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
93      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
94        try {        try {
95             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),             DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
96                                  other.getVector(),                                  other.getVectorRO(),
97                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
98        }        }
99        catch (std::exception& e) {        catch (std::exception& e) {
# Line 92  DataExpanded::DataExpanded(const DataTag Line 105  DataExpanded::DataExpanded(const DataTag
105    
106  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
107                             const DataTypes::RegionType& region)                             const DataTypes::RegionType& region)
108    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
109  {  {
110    //    //
111    // get the shape of the slice    // get the shape of the slice
# Line 110  DataExpanded::DataExpanded(const DataExp Line 123  DataExpanded::DataExpanded(const DataExp
123    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
124      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
125        try {        try {
126  //         getPointDataView().copySlice(getPointOffset(i,j),          DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
127  //                                      other.getPointDataView(),                                       other.getVectorRO(),
 //                                      other.getPointOffset(i,j),  
 //                                      region_loop_range);  
         DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),  
                                      other.getVector(),  
128                       other.getShape(),                       other.getShape(),
129                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
130                                       region_loop_range);                                       region_loop_range);
# Line 127  DataExpanded::DataExpanded(const DataExp Line 136  DataExpanded::DataExpanded(const DataExp
136    }    }
137  }  }
138    
 // 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);  
 // }  
   
139  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
140                             const DataTypes::ShapeType &shape,                             const DataTypes::ShapeType &shape,
141                             const DataTypes::ValueType &data)                             const DataTypes::ValueType &data)
142    : DataAbstract(what,shape)    : parent(what,shape)
143  {  {
144    EsysAssert(data.size()%getNoValues()==0,    EsysAssert(data.size()%getNoValues()==0,
145                   "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 153  DataExpanded::DataExpanded(const Functio
153       // now we copy this value to all elements       // now we copy this value to all elements
154       for (int i=0;i<getLength();)       for (int i=0;i<getLength();)
155       {       {
156      for (int j=0;j<getNoValues();++j,++i)      for (unsigned int j=0;j<getNoValues();++j,++i)
157      {      {
158          vec[i]=data[j];          vec[i]=data[j];
159      }      }
# Line 175  DataExpanded::DataExpanded(const Functio Line 169  DataExpanded::DataExpanded(const Functio
169    
170  }  }
171    
172    DataExpanded::DataExpanded(const FunctionSpace& what,
173                               const DataTypes::ShapeType &shape,
174                               const double v)
175      : parent(what,shape)
176    {
177         ValueType& vec=m_data.getData();
178         //
179         // create the view of the data
180         initialise(what.getNumSamples(),what.getNumDPPSample());
181         // now we copy this value to all elements
182         const int L=getLength();
183         int i;
184         #pragma omp parallel for schedule(static) private(i)
185         for (i=0;i<L;++i)
186         {
187            vec[i]=v;
188         }
189    }
190    
191    
192    
193  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
194  {  {
195  }  }
# Line 200  DataExpanded::setSlice(const DataAbstrac Line 215  DataExpanded::setSlice(const DataAbstrac
215    if (tempDataExp==0) {    if (tempDataExp==0) {
216      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
217    }    }
218      CHECK_FOR_EX_WRITE
219    //    //
220    // get shape of slice    // get shape of slice
221    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
# Line 218  DataExpanded::setSlice(const DataAbstrac Line 234  DataExpanded::setSlice(const DataAbstrac
234    DataTypes::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
235    DataTypes::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
236    int i, j;    int i, j;
237    ValueType& vec=getVector();    ValueType& vec=getVectorRW();
238    const ShapeType& mshape=getShape();    const ShapeType& mshape=getShape();
239    const ValueType& tVec=tempDataExp->getVector();    const ValueType& tVec=tempDataExp->getVectorRO();
240    const ShapeType& tShape=tempDataExp->getShape();    const ShapeType& tShape=tempDataExp->getShape();
241    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
242    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
243      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);*/  
244          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
245                                         tVec,                                         tVec,
246                         tShape,                         tShape,
# Line 253  DataExpanded::copy(const DataConstant& v Line 265  DataExpanded::copy(const DataConstant& v
265    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
266    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
267      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
268        // 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);  
269      }      }
270    }    }
271  }  }
272    
273    void
274  // void  DataExpanded::copy(const WrappedArray& value)
275  // 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);  
   
   //  
276    // check the input shape matches this shape    // check the input shape matches this shape
277    if (!DataTypes::checkShape(getShape(),tempShape)) {    if (!DataTypes::checkShape(getShape(),value.getShape())) {
278      throw DataException(DataTypes::createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
279                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
280                          tempShape,getShape()));                          value.getShape(),getShape()));
281    }    }
282    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the data  
   //copy(temp_dataView);  
   getVector().copyFromNumArray(value);  
283  }  }
284    
285    
# Line 327  DataExpanded::initialise(int noSamples, Line 296  DataExpanded::initialise(int noSamples,
296    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
297  }  }
298    
299    bool
300    DataExpanded::hasNaN() const
301    {
302        const ValueType& v=m_data.getData();
303        for (ValueType::size_type i=0;i<v.size();++i)
304        {
305            if (nancheck(v[i]))
306            {
307                return true;
308            }
309        }
310        return false;
311    }
312    
313    
314  string  string
315  DataExpanded::toString() const  DataExpanded::toString() const
316  {  {
# Line 339  DataExpanded::toString() const Line 323  DataExpanded::toString() const
323        offset=getPointOffset(i,j);        offset=getPointOffset(i,j);
324        stringstream suffix;        stringstream suffix;
325        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
326        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());        temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
327        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
328          temp << endl;          temp << endl;
329        }        }
# Line 361  DataExpanded::getPointOffset(int sampleN Line 345  DataExpanded::getPointOffset(int sampleN
345  }  }
346    
347  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
348    DataExpanded::getPointOffset(int sampleNo,
349                                 int dataPointNo)
350    {
351      return m_data.index(sampleNo,dataPointNo);
352    }
353    
354    DataTypes::ValueType::size_type
355  DataExpanded::getLength() const  DataExpanded::getLength() const
356  {  {
357    return m_data.size();    return m_data.size();
# Line 370  DataExpanded::getLength() const Line 361  DataExpanded::getLength() const
361    
362  void  void
363  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
364      CHECK_FOR_EX_WRITE
365    //    //
366    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
367    int numSamples = getNumSamples();    int numSamples = getNumSamples();
# Line 385  DataExpanded::copyToDataPoint(const int Line 377  DataExpanded::copyToDataPoint(const int
377             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
378       }       }
379       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
380       ValueType& vec=getVector();       ValueType& vec=getVectorRW();
381       if (dataPointRank==0) {       if (dataPointRank==0) {
382           vec[0]=value;           vec[offset]=value;
383       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
384          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
385              vec[offset+i]=value;              vec[offset+i]=value;
# Line 419  DataExpanded::copyToDataPoint(const int Line 411  DataExpanded::copyToDataPoint(const int
411       }       }
412    }    }
413  }  }
414    
415  void  void
416  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) {
417      CHECK_FOR_EX_WRITE
418    //    //
419    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
420    int numSamples = getNumSamples();    int numSamples = getNumSamples();
421    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
   int dataPointRank = getRank();  
   const ShapeType& shape = getShape();  
422    //    //
423    // check rank:    // check rank:
424    if (value.getrank()!=dataPointRank)    if (value.getRank()!=getRank())
425         throw DataException("Rank of numarray does not match Data object rank");         throw DataException("Rank of value does not match Data object rank");
426    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
427       //TODO: global error handling       //TODO: global error handling
428       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 439  DataExpanded::copyToDataPoint(const int Line 431  DataExpanded::copyToDataPoint(const int
431       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
432             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
433       }       }
434       ValueType& vec=getVector();       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
435       if (dataPointRank==0) {       ValueType& vec=getVectorRW();
436           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++;  
     }  
437    }    }
438  }  }
439    
440  void  void
441  DataExpanded::symmetric(DataAbstract* ev)  DataExpanded::symmetric(DataAbstract* ev)
442  {  {
# Line 538  DataExpanded::symmetric(DataAbstract* ev Line 447  DataExpanded::symmetric(DataAbstract* ev
447    if (temp_ev==0) {    if (temp_ev==0) {
448      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).");
449    }    }
450    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
451    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
452    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
453    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
454    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
455    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 560  DataExpanded::nonsymmetric(DataAbstract* Line 469  DataExpanded::nonsymmetric(DataAbstract*
469    if (temp_ev==0) {    if (temp_ev==0) {
470      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).");
471    }    }
472    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
473    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
474    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
475    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
476    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
477    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 582  DataExpanded::trace(DataAbstract* ev, in Line 491  DataExpanded::trace(DataAbstract* ev, in
491    if (temp_ev==0) {    if (temp_ev==0) {
492      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).");
493    }    }
494    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
495    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
496    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
497    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
498    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
499    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 605  DataExpanded::transpose(DataAbstract* ev Line 514  DataExpanded::transpose(DataAbstract* ev
514    if (temp_ev==0) {    if (temp_ev==0) {
515      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).");
516    }    }
517    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
518    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
519    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
520    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
521    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
522    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 628  DataExpanded::swapaxes(DataAbstract* ev, Line 537  DataExpanded::swapaxes(DataAbstract* ev,
537    if (temp_ev==0) {    if (temp_ev==0) {
538      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).");
539    }    }
540    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
541    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
542    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
543    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
544    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
545    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 650  DataExpanded::eigenvalues(DataAbstract* Line 559  DataExpanded::eigenvalues(DataAbstract*
559    if (temp_ev==0) {    if (temp_ev==0) {
560      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).");
561    }    }
562    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
563    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
564    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
565    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
566    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
567    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 676  DataExpanded::eigenvalues_and_eigenvecto Line 585  DataExpanded::eigenvalues_and_eigenvecto
585    if (temp_V==0) {    if (temp_V==0) {
586      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).");
587    }    }
588    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
589    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
590    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
591    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
592    ValueType& VVec=temp_V->getVector();    ValueType& VVec=temp_V->getVectorRW();
593    const ShapeType& VShape=temp_V->getShape();    const ShapeType& VShape=temp_V->getShape();
594    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
595    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 692  DataExpanded::eigenvalues_and_eigenvecto Line 601  DataExpanded::eigenvalues_and_eigenvecto
601    }    }
602  }  }
603    
604    
605    int
606    DataExpanded::matrixInverse(DataAbstract* out) const
607    {
608      DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
609      if (temp==0)
610      {
611        throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
612      }
613    
614      if (getRank()!=2)
615      {
616        throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
617    
618      }
619      int  sampleNo;
620      const int numdpps=getNumDPPSample();
621      const int numSamples = getNumSamples();
622      const ValueType& vec=m_data.getData();
623      int errcode=0;
624      #pragma omp parallel private(sampleNo)
625      {
626         int errorcode=0;
627         LapackInverseHelper h(getShape()[0]);
628         #pragma omp for schedule(static)
629         for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
630         {
631                // not sure I like all those virtual calls to getPointOffset
632            DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
633            int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
634        if (res>errorcode)
635        {
636            errorcode=res;
637            #pragma omp critical
638            {
639              errcode=errorcode;    // I'm not especially concerned which error gets reported as long as one is
640            }
641        }
642         }
643      }
644      return errcode;
645      if (errcode)
646      {
647        DataMaths::matrixInverseError(errcode); // throws exceptions
648      }
649    }
650    
651  void  void
652  DataExpanded::setToZero(){  DataExpanded::setToZero(){
653  // TODO: Surely there is a more efficient way to do this????  // TODO: Surely there is a more efficient way to do this????
654  // Why is there no memset here? Parallel issues?  // Why is there no memset here? Parallel issues?
655    // A: This ensures that memory is touched by the correct thread.
656      CHECK_FOR_EX_WRITE
657    int numSamples = getNumSamples();    int numSamples = getNumSamples();
658    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
659    DataTypes::ValueType::size_type n = getNoValues();    DataTypes::ValueType::size_type n = getNoValues();
# Line 732  DataExpanded::dump(const std::string fil Line 690  DataExpanded::dump(const std::string fil
690     long dims[ldims];     long dims[ldims];
691     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
692     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
693     int mpi_iam=getFunctionSpace().getDomain().getMPIRank();     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
694     int mpi_num=getFunctionSpace().getDomain().getMPISize();     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
695  #ifdef PASO_MPI  #ifdef ESYS_MPI
696     MPI_Status status;     MPI_Status status;
697  #endif  #endif
698    
699  #ifdef PASO_MPI  #ifdef ESYS_MPI
700     /* Serialize NetCDF I/O */     /* Serialize NetCDF I/O */
701     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);
702  #endif  #endif
# Line 798  DataExpanded::dump(const std::string fil Line 756  DataExpanded::dump(const std::string fil
756       if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
757          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
758     }     }
759  #ifdef PASO_MPI  #ifdef ESYS_MPI
760     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);
761  #endif  #endif
762     #else     #else
# Line 812  DataExpanded::setTaggedValue(int tagKey, Line 770  DataExpanded::setTaggedValue(int tagKey,
770                 const DataTypes::ValueType& value,                 const DataTypes::ValueType& value,
771             int dataOffset)             int dataOffset)
772  {  {
773      CHECK_FOR_EX_WRITE
774    int numSamples = getNumSamples();    int numSamples = getNumSamples();
775    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
776    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
# Line 838  DataExpanded::setTaggedValue(int tagKey, Line 797  DataExpanded::setTaggedValue(int tagKey,
797  void  void
798  DataExpanded::reorderByReferenceIDs(int *reference_ids)  DataExpanded::reorderByReferenceIDs(int *reference_ids)
799  {  {
800      CHECK_FOR_EX_WRITE
801    int numSamples = getNumSamples();    int numSamples = getNumSamples();
802    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
803    int sampleNo, sampleNo2,i;    int sampleNo, sampleNo2,i;
# Line 873  DataExpanded::reorderByReferenceIDs(int Line 833  DataExpanded::reorderByReferenceIDs(int
833  }  }
834    
835  DataTypes::ValueType&  DataTypes::ValueType&
836  DataExpanded::getVector()  DataExpanded::getVectorRW()
837  {  {
838        CHECK_FOR_EX_WRITE
839      return m_data.getData();      return m_data.getData();
840  }  }
841    
842  const DataTypes::ValueType&  const DataTypes::ValueType&
843  DataExpanded::getVector() const  DataExpanded::getVectorRO() const
844  {  {
845      return m_data.getData();      return m_data.getData();
846  }  }
847    
848    
849    // Idea here is to create an array of seeds by feeding the original seed into the random generator
850    // The code at the beginning of the function to compute the seed if one is given is
851    // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
852    // I make no claim about how well these initial seeds are distributed
853    void DataExpanded::randomFill(long seed)
854    {
855        CHECK_FOR_EX_WRITE
856        static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
857        if (seed==0)        // for each one
858        {
859        if (prevseed==0)
860        {
861            time_t s=time(0);
862            seed=s;
863        }
864        else
865        {
866            seed=prevseed+419;  // these numbers are arbitrary
867            if (seed>3040101)       // I want to avoid overflow on 32bit systems
868            {
869            seed=((int)(seed)%0xABCD)+1;
870            }
871        }
872        }
873        // now we need to consider MPI since we don't want each rank to start with the same seed
874        seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
875        prevseed=seed;
876    #ifdef _OPENMP
877        int numthreads=omp_get_max_threads();
878    #else
879        int numthreads=1;
880    #endif
881    
882    #ifdef MKLRANDOM
883        double* seeds=new double[numthreads];
884        VSLStreamStatePtr sstream;
885    
886        int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed);  // use a Mersenne Twister
887        numeric_limits<double> dlim;
888        vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
889        vslDeleteStream(&sstream);
890        DataVector& dv=getVectorRW();
891        size_t dvsize=dv.size();
892        #pragma omp parallel
893        {
894        int tnum=0;
895        #ifdef _OPENMP
896        tnum=omp_get_thread_num();
897        #endif
898        VSLStreamStatePtr stream;
899        // the 12345 is a hack to give us a better chance of getting different integer seeds.
900            int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345);  // use a Mersenne Twister
901        int bigchunk=(dvsize/numthreads+1);
902        int smallchunk=dvsize-bigchunk*(numthreads-1);
903        int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
904            vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
905            vslDeleteStream(&stream);
906        }
907        delete[] seeds;
908    #else
909        srand(seed);
910        unsigned* seeds=new unsigned[numthreads];
911        for (int i=0;i<numthreads;++i)
912        {
913        seeds[i]=rand();
914        }
915        DataVector&  dv=getVectorRW();
916        long i;
917        const size_t dvsize=dv.size();
918        #pragma omp parallel private(i)
919        {
920        int tnum=0;
921        #ifdef _OPENMP
922        tnum=omp_get_thread_num();
923        #endif
924        unsigned info=seeds[tnum];
925        
926            #pragma omp for schedule(static)
927            for (i=0;i<dvsize;++i)
928            {
929            dv[i]=(double)rand_r(&info)/RAND_MAX;
930            }
931        }
932        delete[] seeds;
933    #endif
934    }
935    
936  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26