/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 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 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 200  DataExpanded::setSlice(const DataAbstrac Line 186  DataExpanded::setSlice(const DataAbstrac
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    // get shape of slice    // get shape of slice
192    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
# Line 218  DataExpanded::setSlice(const DataAbstrac Line 205  DataExpanded::setSlice(const DataAbstrac
205    DataTypes::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
206    DataTypes::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
207    int i, j;    int i, j;
208    ValueType& vec=getVector();    ValueType& vec=getVectorRW();
209    const ShapeType& mshape=getShape();    const ShapeType& mshape=getShape();
210    const ValueType& tVec=tempDataExp->getVector();    const ValueType& tVec=tempDataExp->getVectorRO();
211    const ShapeType& tShape=tempDataExp->getShape();    const ShapeType& tShape=tempDataExp->getShape();
212    #pragma omp parallel for private(i,j) schedule(static)    #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++) {
 /*      getPointDataView().copySliceFrom(getPointOffset(i,j),  
                                        tempDataExp->getPointDataView(),  
                                        tempDataExp->getPointOffset(i,j),  
                                        region_loop_range);*/  
215          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),          DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
216                                         tVec,                                         tVec,
217                         tShape,                         tShape,
# Line 253  DataExpanded::copy(const DataConstant& v Line 236  DataExpanded::copy(const DataConstant& v
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        // 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);  
240      }      }
241    }    }
242  }  }
243    
244    void
245  // void  DataExpanded::copy(const WrappedArray& value)
246  // 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);  
   
   //  
247    // check the input shape matches this shape    // check the input shape matches this shape
248    if (!DataTypes::checkShape(getShape(),tempShape)) {    if (!DataTypes::checkShape(getShape(),value.getShape())) {
249      throw DataException(DataTypes::createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
250                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
251                          tempShape,getShape()));                          value.getShape(),getShape()));
252    }    }
253    //    getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
   // now copy over the data  
   //copy(temp_dataView);  
   getVector().copyFromNumArray(value);  
254  }  }
255    
256    
# Line 327  DataExpanded::initialise(int noSamples, Line 267  DataExpanded::initialise(int noSamples,
267    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());    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  {  {
# Line 339  DataExpanded::toString() const Line 294  DataExpanded::toString() const
294        offset=getPointOffset(i,j);        offset=getPointOffset(i,j);
295        stringstream suffix;        stringstream suffix;
296        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
297        temp << DataTypes::pointToString(getVector(),getShape(),offset,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        }        }
# Line 361  DataExpanded::getPointOffset(int sampleN Line 316  DataExpanded::getPointOffset(int sampleN
316  }  }
317    
318  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
319    DataExpanded::getPointOffset(int sampleNo,
320                                 int dataPointNo)
321    {
322      return m_data.index(sampleNo,dataPointNo);
323    }
324    
325    DataTypes::ValueType::size_type
326  DataExpanded::getLength() const  DataExpanded::getLength() const
327  {  {
328    return m_data.size();    return m_data.size();
# Line 370  DataExpanded::getLength() const Line 332  DataExpanded::getLength() const
332    
333  void  void
334  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
335      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();
# Line 385  DataExpanded::copyToDataPoint(const int Line 348  DataExpanded::copyToDataPoint(const int
348             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
349       }       }
350       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
351       ValueType& vec=getVector();       ValueType& vec=getVectorRW();
352       if (dataPointRank==0) {       if (dataPointRank==0) {
353           vec[0]=value;           vec[offset]=value;
354       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
355          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
356              vec[offset+i]=value;              vec[offset+i]=value;
# Line 419  DataExpanded::copyToDataPoint(const int Line 382  DataExpanded::copyToDataPoint(const int
382       }       }
383    }    }
384  }  }
385    
386  void  void
387  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) {
388      CHECK_FOR_EX_WRITE
389    //    //
390    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
391    int numSamples = getNumSamples();    int numSamples = getNumSamples();
392    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
   int dataPointRank = getRank();  
   const ShapeType& shape = getShape();  
393    //    //
394    // check rank:    // check rank:
395    if (value.getrank()!=dataPointRank)    if (value.getRank()!=getRank())
396         throw DataException("Rank of numarray does not match Data object rank");         throw DataException("Rank of value does not match Data object rank");
397    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
398       //TODO: global error handling       //TODO: global error handling
399       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 439  DataExpanded::copyToDataPoint(const int Line 402  DataExpanded::copyToDataPoint(const int
402       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
403             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
404       }       }
405       ValueType& vec=getVector();       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
406       if (dataPointRank==0) {       ValueType& vec=getVectorRW();
407           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++;  
     }  
408    }    }
409  }  }
410    
411  void  void
412  DataExpanded::symmetric(DataAbstract* ev)  DataExpanded::symmetric(DataAbstract* ev)
413  {  {
# Line 538  DataExpanded::symmetric(DataAbstract* ev Line 418  DataExpanded::symmetric(DataAbstract* ev
418    if (temp_ev==0) {    if (temp_ev==0) {
419      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).");
420    }    }
421    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
422    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
423    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
424    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
425    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
426    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 560  DataExpanded::nonsymmetric(DataAbstract* Line 440  DataExpanded::nonsymmetric(DataAbstract*
440    if (temp_ev==0) {    if (temp_ev==0) {
441      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).");
442    }    }
443    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
444    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
445    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
446    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
447    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
448    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 582  DataExpanded::trace(DataAbstract* ev, in Line 462  DataExpanded::trace(DataAbstract* ev, in
462    if (temp_ev==0) {    if (temp_ev==0) {
463      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).");
464    }    }
465    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
466    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
467    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
468    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
469    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
470    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 605  DataExpanded::transpose(DataAbstract* ev Line 485  DataExpanded::transpose(DataAbstract* ev
485    if (temp_ev==0) {    if (temp_ev==0) {
486      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).");
487    }    }
488    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
489    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
490    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
491    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
492    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
493    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 628  DataExpanded::swapaxes(DataAbstract* ev, Line 508  DataExpanded::swapaxes(DataAbstract* ev,
508    if (temp_ev==0) {    if (temp_ev==0) {
509      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).");
510    }    }
511    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
512    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
513    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
514    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
515    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
516    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 650  DataExpanded::eigenvalues(DataAbstract* Line 530  DataExpanded::eigenvalues(DataAbstract*
530    if (temp_ev==0) {    if (temp_ev==0) {
531      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).");
532    }    }
533    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
534    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
535    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
536    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
537    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 676  DataExpanded::eigenvalues_and_eigenvecto Line 556  DataExpanded::eigenvalues_and_eigenvecto
556    if (temp_V==0) {    if (temp_V==0) {
557      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).");
558    }    }
559    ValueType& vec=getVector();    const ValueType& vec=getVectorRO();
560    const ShapeType& shape=getShape();    const ShapeType& shape=getShape();
561    ValueType& evVec=temp_ev->getVector();    ValueType& evVec=temp_ev->getVectorRW();
562    const ShapeType& evShape=temp_ev->getShape();    const ShapeType& evShape=temp_ev->getShape();
563    ValueType& VVec=temp_V->getVector();    ValueType& VVec=temp_V->getVectorRW();
564    const ShapeType& VShape=temp_V->getShape();    const ShapeType& VShape=temp_V->getShape();
565    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
566    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 692  DataExpanded::eigenvalues_and_eigenvecto Line 572  DataExpanded::eigenvalues_and_eigenvecto
572    }    }
573  }  }
574    
575    
576    int
577    DataExpanded::matrixInverse(DataAbstract* out) const
578    {
579      DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
580      if (temp==0)
581      {
582        throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
583      }
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  void
623  DataExpanded::setToZero(){  DataExpanded::setToZero(){
624  // TODO: Surely there is a more efficient way to do this????  // TODO: Surely there is a more efficient way to do this????
625  // Why is there no memset here? Parallel issues?  // 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();    int numSamples = getNumSamples();
629    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
630    DataTypes::ValueType::size_type n = getNoValues();    DataTypes::ValueType::size_type n = getNoValues();
# Line 732  DataExpanded::dump(const std::string fil Line 661  DataExpanded::dump(const std::string fil
661     long dims[ldims];     long dims[ldims];
662     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
663     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
664     int mpi_iam=getFunctionSpace().getDomain().getMPIRank();     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
665     int mpi_num=getFunctionSpace().getDomain().getMPISize();     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
666  #ifdef PASO_MPI  #ifdef PASO_MPI
667     MPI_Status status;     MPI_Status status;
668  #endif  #endif
# Line 812  DataExpanded::setTaggedValue(int tagKey, Line 741  DataExpanded::setTaggedValue(int tagKey,
741                 const DataTypes::ValueType& value,                 const DataTypes::ValueType& value,
742             int dataOffset)             int dataOffset)
743  {  {
744      CHECK_FOR_EX_WRITE
745    int numSamples = getNumSamples();    int numSamples = getNumSamples();
746    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
747    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
# Line 838  DataExpanded::setTaggedValue(int tagKey, Line 768  DataExpanded::setTaggedValue(int tagKey,
768  void  void
769  DataExpanded::reorderByReferenceIDs(int *reference_ids)  DataExpanded::reorderByReferenceIDs(int *reference_ids)
770  {  {
771      CHECK_FOR_EX_WRITE
772    int numSamples = getNumSamples();    int numSamples = getNumSamples();
773    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();    DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
774    int sampleNo, sampleNo2,i;    int sampleNo, sampleNo2,i;
# Line 873  DataExpanded::reorderByReferenceIDs(int Line 804  DataExpanded::reorderByReferenceIDs(int
804  }  }
805    
806  DataTypes::ValueType&  DataTypes::ValueType&
807  DataExpanded::getVector()  DataExpanded::getVectorRW()
808  {  {
809        CHECK_FOR_EX_WRITE
810      return m_data.getData();      return m_data.getData();
811  }  }
812    
813  const DataTypes::ValueType&  const DataTypes::ValueType&
814  DataExpanded::getVector() const  DataExpanded::getVectorRO() const
815  {  {
816      return m_data.getData();      return m_data.getData();
817  }  }
818    
819    
820  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26