/[escript]/branches/trilinos_from_5897/escriptcore/src/DataExpanded.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/escriptcore/src/DataExpanded.cpp

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

revision 5962 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 5963 by caltinay, Mon Feb 22 06:59:27 2016 UTC
# Line 41  using namespace escript::DataTypes; Line 41  using namespace escript::DataTypes;
41          abort();\          abort();\
42          throw DataException(ss.str());\          throw DataException(ss.str());\
43      }\      }\
44  } while(0)  } while(0);
45    
46    
47  namespace escript {  namespace escript {
# Line 51  DataExpanded::DataExpanded(const Wrapped Line 51  DataExpanded::DataExpanded(const Wrapped
51    : parent(what,value.getShape())    : parent(what,value.getShape())
52  {  {
53      // initialise the data array for this object      // initialise the data array for this object
54      initialise(what.getNumSamples(),what.getNumDPPSample());      initialise(what.getNumSamples(),what.getNumDPPSample(), value.isComplex());
55      // copy the given value to every data point      // copy the given value to every data point
56      copy(value);      copy(value);
57  }  }
58    
59  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
60    : parent(other.getFunctionSpace(), other.getShape()),    : parent(other.getFunctionSpace(), other.getShape()),
61      m_data(other.m_data)      m_data_r(other.m_data_r), m_data_c(other.m_data_c)
62  {  {
63        m_iscompl=other.m_iscompl;
64  }  }
65    
66  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
67    : parent(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
68  {  {
69      // initialise the data array for this object      // initialise the data array for this object
70      initialise(other.getNumSamples(),other.getNumDPPSample());      initialise(other.getNumSamples(),other.getNumDPPSample(), other.isComplex());
71      // DataConstant only has one value, copy this to every data point      // DataConstant only has one value, copy this to every data point
72      copy(other);      copy(other);
73  }  }
# Line 75  DataExpanded::DataExpanded(const DataTag Line 76  DataExpanded::DataExpanded(const DataTag
76    : parent(other.getFunctionSpace(), other.getShape())    : parent(other.getFunctionSpace(), other.getShape())
77  {  {
78      // initialise the data array for this object      // initialise the data array for this object
79      initialise(other.getNumSamples(),other.getNumDPPSample());      initialise(other.getNumSamples(),other.getNumDPPSample(), other.isComplex());
80      // for each data point in this object, extract and copy the corresponding      // for each data point in this object, extract and copy the corresponding
81      // data value from the given DataTagged object      // data value from the given DataTagged object
82      DataTypes::ValueType::size_type numRows=m_data.getNumRows();      if (isComplex())
83      DataTypes::ValueType::size_type numCols=m_data.getNumCols();      {
84  #pragma omp parallel for      DataTypes::cplx_t dummy=0;
85      for (int i=0; i<numRows; i++) {      #pragma omp parallel for
86          for (int j=0; j<numCols; j++) {      for (int i=0; i<m_noSamples; i++) {
87              try {          for (int j=0; j<m_noDataPointsPerSample; j++) {
88                  DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j),          try {
89                                       getNoValues(), other.getVectorRO(),              DataTypes::copyPoint(getTypedVectorRW(dummy), getPointOffset(i,j),
90                                       other.getPointOffset(i,j));                      getNoValues(), other.getTypedVectorRO(dummy),
91              } catch (std::exception& e) {                      other.getPointOffset(i,j));
92                  cerr << e.what() << endl;          } catch (std::exception& e) {
93              }              cerr << e.what() << endl;
94          }          }
95            }
96        }      
97          
98          
99        }
100        else
101        {
102        DataTypes::real_t dummy=0;
103        #pragma omp parallel for
104        for (int i=0; i<m_noSamples; i++) {
105            for (int j=0; j<m_noDataPointsPerSample; j++) {
106            try {
107                DataTypes::copyPoint(getTypedVectorRW(dummy), getPointOffset(i,j),
108                        getNoValues(), other.getTypedVectorRO(dummy),
109                        other.getPointOffset(i,j));
110            } catch (std::exception& e) {
111                cerr << e.what() << endl;
112            }
113            }
114        }
115      }      }
116  }  }
117    
# Line 99  DataExpanded::DataExpanded(const DataExp Line 120  DataExpanded::DataExpanded(const DataExp
120    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
121  {  {
122      // initialise this Data object to the shape of the slice      // initialise this Data object to the shape of the slice
123      initialise(other.getNumSamples(),other.getNumDPPSample());      initialise(other.getNumSamples(),other.getNumDPPSample(), other.isComplex());
124      // copy the data      // copy the data
125      DataTypes::RegionLoopRangeType region_loop_range =      DataTypes::RegionLoopRangeType region_loop_range =
126                                      DataTypes::getSliceRegionLoopRange(region);                                      DataTypes::getSliceRegionLoopRange(region);
127      DataTypes::ValueType::size_type numRows=m_data.getNumRows();      if (isComplex())
128      DataTypes::ValueType::size_type numCols=m_data.getNumCols();      {
129  #pragma omp parallel for      DataTypes::cplx_t dummy=0;
130      for (int i=0; i<numRows; i++) {      #pragma omp parallel for
131          for (int j=0; j<numCols; j++) {      for (int i=0; i<m_noSamples; i++) {
132              try {          for (int j=0; j<m_noDataPointsPerSample; j++) {
133                  DataTypes::copySlice(getVectorRW(), getShape(),          try {
134                                       getPointOffset(i,j), other.getVectorRO(),              DataTypes::copySlice(getTypedVectorRW(dummy), getShape(),
135                                       other.getShape(),                      getPointOffset(i,j), other.getTypedVectorRO(dummy),
136                                       other.getPointOffset(i,j),                      other.getShape(),
137                                       region_loop_range);                      other.getPointOffset(i,j),
138              } catch (std::exception& e) {                      region_loop_range);
139                  cerr << e.what() << endl;          } catch (std::exception& e) {
140                cerr << e.what() << endl;
141            }
142            }
143        }      
144        }
145        else
146        {
147        DataTypes::real_t dummy=0;
148        #pragma omp parallel for
149        for (int i=0; i<m_noSamples; i++) {
150            for (int j=0; j<m_noDataPointsPerSample; j++) {
151            try {
152                DataTypes::copySlice(getTypedVectorRW(dummy), getShape(),
153                        getPointOffset(i,j), other.getTypedVectorRO(dummy),
154                        other.getShape(),
155                        other.getPointOffset(i,j),
156                        region_loop_range);
157            } catch (std::exception& e) {
158                cerr << e.what() << endl;
159            }
160            }
161        }
162        }
163    }
164    
165    DataExpanded::DataExpanded(const FunctionSpace& what,
166                               const DataTypes::ShapeType &shape,
167                               const DataTypes::RealVectorType &data)
168      : parent(what,shape)
169    {
170        EsysAssert(data.size()%getNoValues()==0,
171                     "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
172    
173        if (data.size() == getNoValues()) {
174            RealVectorType& vec=m_data_r;
175            // create the view of the data
176            initialise(what.getNumSamples(),what.getNumDPPSample(), false);
177            // now we copy this value to all elements
178            for (int i=0; i<getLength();) {
179                for (unsigned int j=0;j<getNoValues();++j,++i) {
180                    vec[i]=data[j];
181              }              }
182          }          }
183        } else {
184            // copy the data in the correct format
185            m_data_r = data;
186      }      }
187  }  }
188    
189  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
190                             const DataTypes::ShapeType &shape,                             const DataTypes::ShapeType &shape,
191                             const DataTypes::ValueType &data)                             const DataTypes::CplxVectorType &data)
192    : parent(what,shape)    : parent(what,shape)
193  {  {
194      EsysAssert(data.size()%getNoValues()==0,      EsysAssert(data.size()%getNoValues()==0,
195                   "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.");
196    
197      if (data.size() == getNoValues()) {      if (data.size() == getNoValues()) {
198          ValueType& vec=m_data.getData();          CplxVectorType& vec=m_data_c;
199          // create the view of the data          // create the view of the data
200          initialise(what.getNumSamples(),what.getNumDPPSample());          initialise(what.getNumSamples(),what.getNumDPPSample(), true);
201          // now we copy this value to all elements          // now we copy this value to all elements
202          for (int i=0; i<getLength();) {          for (int i=0; i<getLength();) {
203              for (unsigned int j=0;j<getNoValues();++j,++i) {              for (unsigned int j=0;j<getNoValues();++j,++i) {
# Line 141  DataExpanded::DataExpanded(const Functio Line 206  DataExpanded::DataExpanded(const Functio
206          }          }
207      } else {      } else {
208          // copy the data in the correct format          // copy the data in the correct format
209          m_data.getData() = data;          m_data_c = data;
210      }      }
211  }  }
212    
213    
214  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
215                             const DataTypes::ShapeType &shape,                             const DataTypes::ShapeType &shape,
216                             const double v)                             const DataTypes::real_t v)
217    : parent(what,shape)    : parent(what,shape)
218  {  {
219      ValueType& vec=m_data.getData();      initialise(what.getNumSamples(),what.getNumDPPSample(), false);
220      // create the view of the data      DataTypes::RealVectorType& vec=m_data_r;
     initialise(what.getNumSamples(),what.getNumDPPSample());  
221      // now we copy this value to all elements      // now we copy this value to all elements
222      const int L=getLength();      const int L=getLength();
223  #pragma omp parallel for  #pragma omp parallel for
# Line 161  DataExpanded::DataExpanded(const Functio Line 226  DataExpanded::DataExpanded(const Functio
226      }      }
227  }  }
228    
229    DataExpanded::DataExpanded(const FunctionSpace& what,
230                               const DataTypes::ShapeType &shape,
231                               const DataTypes::cplx_t v)
232      : parent(what,shape)
233    {
234        initialise(what.getNumSamples(),what.getNumDPPSample(), true);
235        DataTypes::CplxVectorType& vec=m_data_c;
236        
237        // now we copy this value to all elements
238        const int L=getLength();
239    #pragma omp parallel for
240        for (int i=0; i<L; ++i) {
241            vec[i]=v;
242        }
243    }
244    
245    
246  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
247  {  {
248  }  }
249    
250  DataAbstract* DataExpanded::deepCopy()  DataAbstract* DataExpanded::deepCopy() const
251  {  {
252      return new DataExpanded(*this);      return new DataExpanded(*this);
253  }  }
# Line 198  void DataExpanded::setSlice(const DataAb Line 280  void DataExpanded::setSlice(const DataAb
280                  shape, value->getShape()));                  shape, value->getShape()));
281    
282      // copy the data from the slice into this object      // copy the data from the slice into this object
283      DataTypes::ValueType::size_type numRows = m_data.getNumRows();      DataTypes::RealVectorType& vec=getVectorRW();
     DataTypes::ValueType::size_type numCols = m_data.getNumCols();  
     ValueType& vec=getVectorRW();  
284      const ShapeType& mshape=getShape();      const ShapeType& mshape=getShape();
285      const ValueType& tVec=tempDataExp->getVectorRO();      const DataTypes::RealVectorType& tVec=tempDataExp->getVectorRO();
286      const ShapeType& tShape=tempDataExp->getShape();      const ShapeType& tShape=tempDataExp->getShape();
287  #pragma omp parallel for  #pragma omp parallel for
288      for (int i=0; i<numRows; i++) {      for (int i=0; i<m_noSamples; i++) {
289          for (int j=0; j<numCols; j++) {          for (int j=0; j<m_noDataPointsPerSample; j++) {
290              DataTypes::copySliceFrom(vec, mshape, getPointOffset(i,j), tVec,              DataTypes::copySliceFrom(vec, mshape, getPointOffset(i,j), tVec,
291                                       tShape, tempDataExp->getPointOffset(i,j),                                       tShape, tempDataExp->getPointOffset(i,j),
292                                       region_loop_range);                                       region_loop_range);
# Line 218  void DataExpanded::copy(const DataConsta Line 298  void DataExpanded::copy(const DataConsta
298  {  {
299      EsysAssert((checkShape(getShape(), value.getShape())),      EsysAssert((checkShape(getShape(), value.getShape())),
300                   createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.", value.getShape(), getShape()));                   createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.", value.getShape(), getShape()));
301        if (isComplex())
302      // copy a single value to every data point in this object      {
303      int nRows=m_data.getNumRows();      if (value.isComplex())
304      int nCols=m_data.getNumCols();      {
305  #pragma omp parallel for          // copy a single value to every data point in this object
306      for (int i=0; i<nRows; i++) {          #pragma omp parallel for
307          for (int j=0; j<nCols; j++) {          for (int i=0; i<m_noSamples; i++) {
308              DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j),          for (int j=0; j<m_noDataPointsPerSample; j++) {
309                                   getNoValues(), value.getVectorRO(), 0);              DataTypes::copyPoint(getTypedVectorRW((cplx_t)(0)), getPointOffset(i,j),
310          }                      getNoValues(), value.getTypedVectorRO((cplx_t)(0)), 0);
311            }
312            }      
313        }
314        else    // value is real
315        {
316            throw DataException("Programming error - DataExpanded::copy source and target must be the same complexity.");    
317        }
318      }      }
319        else
320        {
321        if (value.isComplex())
322        {
323            throw DataException("Programming error - DataExpanded::copy source and target must be the same complexity.");        
324        }
325        else
326        {
327            real_t dummy=0;
328            // copy a single value to every data point in this object
329            #pragma omp parallel for
330            for (int i=0; i<m_noSamples; i++) {
331            for (int j=0; j<m_noDataPointsPerSample; j++) {
332                DataTypes::copyPoint(getTypedVectorRW(dummy), getPointOffset(i,j),
333                        getNoValues(), value.getTypedVectorRO(dummy), 0);
334            }
335            }
336        }
337        }
338        
339  }  }
340    
341  void DataExpanded::copy(const WrappedArray& value)  void DataExpanded::copy(const WrappedArray& value)
# Line 241  void DataExpanded::copy(const WrappedArr Line 348  void DataExpanded::copy(const WrappedArr
348      getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());      getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
349  }  }
350    
351  void DataExpanded::initialise(int noSamples, int noDataPointsPerSample)  void DataExpanded::initialise(int noSamples, int noDataPointsPerSample, bool cplx)
352  {  {
353        this->m_iscompl=cplx;
354      if (noSamples==0) //retain the default empty object      if (noSamples==0) //retain the default empty object
355          return;          return;
356    
357      // resize data array to the required size      if (cplx)
358      m_data.resize(noSamples, noDataPointsPerSample, getNoValues());      {
359  }      // resize data array to the required size
360        m_data_c.resize(noSamples*noDataPointsPerSample*getNoValues(), 0.0, noDataPointsPerSample*getNoValues());      
361  bool DataExpanded::hasNaN() const      }
362  {      else
363      bool haveNaN = false;      {
364      const ValueType& v = m_data.getData();      // resize data array to the required size
365  #pragma omp parallel for      m_data_r.resize(noSamples*noDataPointsPerSample*getNoValues(), 0.0, noDataPointsPerSample*getNoValues());
     for (ValueType::size_type i=0; i<v.size(); ++i) {  
         if (nancheck(v[i])) {  
 #pragma omp critical  
             {  
                 haveNaN=true;  
             }  
         }  
366      }      }
     return haveNaN;  
367  }  }
368    
369  void DataExpanded::replaceNaN(double value)  bool
370    DataExpanded::hasNaN() const
371  {  {
372  #pragma omp parallel for    bool haveNaN=false;
373      for (ValueType::size_type i=0; i<m_data.size(); ++i) {    if (isComplex())
374          if (nancheck(m_data[i])) {    {
375              m_data[i] = value;        #pragma omp parallel for
376          }        for (DataTypes::CplxVectorType::size_type i=0;i<m_data_r.size();++i)
377      }        {
378          if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))
379          {
380              #pragma omp critical
381              {
382              haveNaN=true;
383              }
384          }
385          }
386      }
387      else
388      {
389          #pragma omp parallel for
390          for (DataTypes::RealVectorType::size_type i=0;i<m_data_r.size();++i)
391          {
392          if (std::isnan(m_data_r[i]))
393          {
394              #pragma omp critical
395              {
396              haveNaN=true;
397              }
398          }
399          }
400      }
401      return haveNaN;
402    }
403    
404    
405    void
406    DataExpanded::replaceNaN(DataTypes::real_t value) {
407      CHECK_FOR_EX_WRITE  
408      if (isComplex())
409      {
410          #pragma omp parallel for
411          for (DataTypes::CplxVectorType::size_type i=0;i<m_data_r.size();++i)
412          {
413        if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))  
414        {
415          m_data_c[i] = value;
416        }
417          }
418      }
419      else
420      {
421          #pragma omp parallel for
422          for (DataTypes::RealVectorType::size_type i=0;i<m_data_r.size();++i)
423          {
424        if (std::isnan(m_data_r[i]))  
425        {
426          m_data_r[i] = value;
427        }
428          }    
429      }
430    }
431    
432    void
433    DataExpanded::replaceNaN(DataTypes::cplx_t value) {
434      CHECK_FOR_EX_WRITE  
435      if (isComplex())
436      {
437          #pragma omp parallel for
438          for (DataTypes::CplxVectorType::size_type i=0;i<m_data_r.size();++i)
439          {
440        if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))
441        {
442          m_data_c[i] = value;
443        }
444          }
445      }
446      else
447      {
448          complicate();
449          replaceNaN(value);
450      }
451  }  }
452    
453  string DataExpanded::toString() const  string DataExpanded::toString() const
# Line 282  string DataExpanded::toString() const Line 456  string DataExpanded::toString() const
456      FunctionSpace fs=getFunctionSpace();      FunctionSpace fs=getFunctionSpace();
457    
458      int offset=0;      int offset=0;
459      for (int i=0; i<m_data.getNumRows(); i++) {      for (int i=0; i<m_noSamples; i++) {
460          for (int j=0; j<m_data.getNumCols(); j++) {          for (int j=0; j<m_noDataPointsPerSample; j++) {
461              offset = getPointOffset(i,j);              offset = getPointOffset(i,j);
462              stringstream suffix;              stringstream suffix;
463              suffix << "( id: " << i << ", ref: "              suffix << "( id: " << i << ", ref: "
464                     << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";                     << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
465              ss << DataTypes::pointToString(getVectorRO(), getShape(), offset,              ss << (isComplex()?
466                                             suffix.str());              DataTypes::pointToString(getTypedVectorRO((cplx_t)0), getShape(), offset, suffix.str())
467              if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {             :
468                DataTypes::pointToString(getTypedVectorRO((real_t)0), getShape(), offset, suffix.str()));
469                if (!(i==(m_noSamples-1) && j==(m_noDataPointsPerSample-1))) {
470                  ss << endl;                  ss << endl;
471              }              }
472          }          }
# Line 302  string DataExpanded::toString() const Line 478  string DataExpanded::toString() const
478      return result;      return result;
479  }  }
480    
481  DataTypes::ValueType::size_type DataExpanded::getPointOffset(int sampleNo,  DataTypes::RealVectorType::size_type DataExpanded::getPointOffset(int sampleNo,
482                                                          int dataPointNo) const                                                          int dataPointNo) const
483  {  {
484      return m_data.index(sampleNo,dataPointNo);      DataTypes::RealVectorType::size_type blockSize=getNoValues();
485        EsysAssert((isComplex()?
486              ((sampleNo >= 0) && (dataPointNo >= 0) && (m_data_c.size() > 0))
487            :
488              ((sampleNo >= 0) && (dataPointNo >= 0) && (m_data_r.size() > 0))),
489               "(DataBlocks2D) Index value out of range.");
490        DataTypes::RealVectorType::size_type temp=(sampleNo*m_noDataPointsPerSample+dataPointNo)*blockSize;
491        EsysAssert((isComplex()?
492              (temp <= (m_data_c.size()-blockSize))
493            :
494              (temp <= (m_data_r.size()-blockSize))), "Index value out of range.");
495    
496        return temp;
497  }  }
498    
499  DataTypes::ValueType::size_type DataExpanded::getPointOffset(int sampleNo,  
500                                                               int dataPointNo)  void DataExpanded::complicate()
501  {  {
502      return m_data.index(sampleNo,dataPointNo);      if (!isComplex())
503        {
504            fillComplexFromReal(m_data_r, m_data_c);
505            this->m_iscompl=true;
506            m_data_r.resize(0,0,1);
507        }
508  }  }
509    
510  DataTypes::ValueType::size_type DataExpanded::getLength() const  
511    DataTypes::RealVectorType::size_type DataExpanded::getLength() const
512  {  {
513      return m_data.size();      return m_data_r.size();
514  }  }
515    
516  void DataExpanded::copyToDataPoint(int sampleNo, int dataPointNo, double value)  void DataExpanded::copyToDataPoint(int sampleNo, int dataPointNo, const DataTypes::cplx_t value)
517  {  {
518        if (!isComplex())
519        {
520        throw DataException("Programming error - attempt to set complex value on real data.");
521        }
522      CHECK_FOR_EX_WRITE;      CHECK_FOR_EX_WRITE;
523      // Get the number of samples and data-points per sample.      // Get the number of samples and data-points per sample.
524      int numSamples = getNumSamples();      int numSamples = getNumSamples();
# Line 334  void DataExpanded::copyToDataPoint(int s Line 532  void DataExpanded::copyToDataPoint(int s
532          if (dataPointNo >= numDataPointsPerSample || dataPointNo < 0)          if (dataPointNo >= numDataPointsPerSample || dataPointNo < 0)
533              throw DataException("DataExpanded::copyDataPoint: invalid dataPointNo.");              throw DataException("DataExpanded::copyDataPoint: invalid dataPointNo.");
534    
535          ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);          DataTypes::CplxVectorType::size_type offset = getPointOffset(sampleNo, dataPointNo);
536          ValueType& vec = getVectorRW();          DataTypes::CplxVectorType& vec = getTypedVectorRW(cplx_t(0));
537            if (dataPointRank==0) {
538                vec[offset] = value;
539            } else if (dataPointRank==1) {
540                for (int i=0; i<dataPointShape[0]; i++) {
541                    vec[offset+i] = value;
542                }
543            } else if (dataPointRank==2) {
544                for (int i=0; i<dataPointShape[0]; i++) {
545                    for (int j=0; j<dataPointShape[1]; j++) {
546                        vec[offset+getRelIndex(dataPointShape,i,j)] = value;
547                    }
548                }
549            } else if (dataPointRank==3) {
550                for (int i=0; i<dataPointShape[0]; i++) {
551                    for (int j=0; j<dataPointShape[1]; j++) {
552                        for (int k=0; k<dataPointShape[2]; k++) {
553                            vec[offset+getRelIndex(dataPointShape,i,j,k)] = value;
554                        }
555                    }
556                }
557            } else if (dataPointRank==4) {
558                for (int i=0; i<dataPointShape[0]; i++) {
559                    for (int j=0; j<dataPointShape[1]; j++) {
560                        for (int k=0; k<dataPointShape[2]; k++) {
561                            for (int l=0; l<dataPointShape[3]; l++) {
562                                vec[offset+getRelIndex(dataPointShape,i,j,k,l)] = value;
563                            }
564                        }
565                    }
566                }
567            }
568        }
569    }
570    
571    
572    void DataExpanded::copyToDataPoint(int sampleNo, int dataPointNo, const DataTypes::real_t value)
573    {
574        if (isComplex())
575        {
576        copyToDataPoint(sampleNo, dataPointNo, cplx_t(value));
577        return;
578        }
579        CHECK_FOR_EX_WRITE;
580        // Get the number of samples and data-points per sample.
581        int numSamples = getNumSamples();
582        int numDataPointsPerSample = getNumDPPSample();
583        int dataPointRank = getRank();
584        ShapeType dataPointShape = getShape();
585        if (numSamples*numDataPointsPerSample > 0) {
586            //TODO: global error handling
587            if (sampleNo >= numSamples || sampleNo < 0)
588                throw DataException("DataExpanded::copyDataPoint: invalid sampleNo.");
589            if (dataPointNo >= numDataPointsPerSample || dataPointNo < 0)
590                throw DataException("DataExpanded::copyDataPoint: invalid dataPointNo.");
591    
592            DataTypes::RealVectorType::size_type offset = getPointOffset(sampleNo, dataPointNo);
593            DataTypes::RealVectorType& vec = getVectorRW();
594          if (dataPointRank==0) {          if (dataPointRank==0) {
595              vec[offset] = value;              vec[offset] = value;
596          } else if (dataPointRank==1) {          } else if (dataPointRank==1) {
# Line 388  void DataExpanded::copyToDataPoint(int s Line 643  void DataExpanded::copyToDataPoint(int s
643          if (dataPointNo >= numDataPointsPerSample || dataPointNo < 0)          if (dataPointNo >= numDataPointsPerSample || dataPointNo < 0)
644              throw DataException("DataExpanded::copyDataPoint: invalid dataPointNoInSample.");              throw DataException("DataExpanded::copyDataPoint: invalid dataPointNoInSample.");
645    
646          ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);      if (isComplex())
647          ValueType& vec = getVectorRW();      {
648          vec.copyFromArrayToOffset(value, offset, 1);          DataTypes::CplxVectorType::size_type offset = getPointOffset(sampleNo, dataPointNo);
649            DataTypes::CplxVectorType& vec = getTypedVectorRW(cplx_t(0));
650            vec.copyFromArrayToOffset(value, offset, 1);
651        }
652        else
653        {
654            DataTypes::RealVectorType::size_type offset = getPointOffset(sampleNo, dataPointNo);
655            DataTypes::RealVectorType& vec = getTypedVectorRW(real_t(0));
656            vec.copyFromArrayToOffset(value, offset, 1);
657        }
658      }      }
659  }  }
660    
# Line 402  void DataExpanded::symmetric(DataAbstrac Line 666  void DataExpanded::symmetric(DataAbstrac
666      if (!temp_ev)      if (!temp_ev)
667          throw DataException("DataExpanded::symmetric: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::symmetric: casting to DataExpanded failed (probably a programming error).");
668    
669      const ValueType& vec = getVectorRO();      const DataTypes::RealVectorType& vec = getVectorRO();
670      const ShapeType& shape = getShape();      const ShapeType& shape = getShape();
671      ValueType& evVec = temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec = temp_ev->getVectorRW();
672      const ShapeType& evShape = temp_ev->getShape();      const ShapeType& evShape = temp_ev->getShape();
673  #pragma omp parallel for  #pragma omp parallel for
674      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 424  void DataExpanded::nonsymmetric(DataAbst Line 688  void DataExpanded::nonsymmetric(DataAbst
688      if (!temp_ev)      if (!temp_ev)
689          throw DataException("DataExpanded::nonsymmetric: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::nonsymmetric: casting to DataExpanded failed (probably a programming error).");
690    
691      const ValueType& vec = getVectorRO();      const DataTypes::RealVectorType& vec = getVectorRO();
692      const ShapeType& shape = getShape();      const ShapeType& shape = getShape();
693      ValueType& evVec = temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec = temp_ev->getVectorRW();
694      const ShapeType& evShape = temp_ev->getShape();      const ShapeType& evShape = temp_ev->getShape();
695  #pragma omp parallel for  #pragma omp parallel for
696      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 446  void DataExpanded::trace(DataAbstract* e Line 710  void DataExpanded::trace(DataAbstract* e
710      if (!temp_ev)      if (!temp_ev)
711          throw DataException("DataExpanded::trace: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::trace: casting to DataExpanded failed (probably a programming error).");
712    
713      const ValueType& vec=getVectorRO();      const DataTypes::RealVectorType& vec=getVectorRO();
714      const ShapeType& shape=getShape();      const ShapeType& shape=getShape();
715      ValueType& evVec=temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
716      const ShapeType& evShape=temp_ev->getShape();      const ShapeType& evShape=temp_ev->getShape();
717  #pragma omp parallel for  #pragma omp parallel for
718      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 468  void DataExpanded::transpose(DataAbstrac Line 732  void DataExpanded::transpose(DataAbstrac
732      if (!temp_ev)      if (!temp_ev)
733          throw DataException("DataExpanded::transpose: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::transpose: casting to DataExpanded failed (probably a programming error).");
734    
735      const ValueType& vec = getVectorRO();      const DataTypes::RealVectorType& vec = getVectorRO();
736      const ShapeType& shape = getShape();      const ShapeType& shape = getShape();
737      ValueType& evVec = temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec = temp_ev->getVectorRW();
738      const ShapeType& evShape = temp_ev->getShape();      const ShapeType& evShape = temp_ev->getShape();
739  #pragma omp parallel for  #pragma omp parallel for
740      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 490  void DataExpanded::swapaxes(DataAbstract Line 754  void DataExpanded::swapaxes(DataAbstract
754      if (!temp_ev)      if (!temp_ev)
755          throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (probably a programming error).");          throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (probably a programming error).");
756    
757      const ValueType& vec=getVectorRO();      const DataTypes::RealVectorType& vec=getVectorRO();
758      const ShapeType& shape=getShape();      const ShapeType& shape=getShape();
759      ValueType& evVec=temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
760      const ShapeType& evShape=temp_ev->getShape();      const ShapeType& evShape=temp_ev->getShape();
761  #pragma omp parallel for  #pragma omp parallel for
762      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 512  void DataExpanded::eigenvalues(DataAbstr Line 776  void DataExpanded::eigenvalues(DataAbstr
776      if (!temp_ev)      if (!temp_ev)
777          throw DataException("DataExpanded::eigenvalues: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::eigenvalues: casting to DataExpanded failed (probably a programming error).");
778    
779      const ValueType& vec=getVectorRO();      const DataTypes::RealVectorType& vec=getVectorRO();
780      const ShapeType& shape=getShape();      const ShapeType& shape=getShape();
781      ValueType& evVec=temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
782      const ShapeType& evShape=temp_ev->getShape();      const ShapeType& evShape=temp_ev->getShape();
783  #pragma omp parallel for  #pragma omp parallel for
784      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 539  void DataExpanded::eigenvalues_and_eigen Line 803  void DataExpanded::eigenvalues_and_eigen
803      if (!temp_V)      if (!temp_V)
804          throw DataException("DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (probably a programming error).");          throw DataException("DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (probably a programming error).");
805    
806      const ValueType& vec = getVectorRO();      const DataTypes::RealVectorType& vec = getVectorRO();
807      const ShapeType& shape = getShape();      const ShapeType& shape = getShape();
808      ValueType& evVec = temp_ev->getVectorRW();      DataTypes::RealVectorType& evVec = temp_ev->getVectorRW();
809      const ShapeType& evShape = temp_ev->getShape();      const ShapeType& evShape = temp_ev->getShape();
810      ValueType& VVec = temp_V->getVectorRW();      DataTypes::RealVectorType& VVec = temp_V->getVectorRW();
811      const ShapeType& VShape = temp_V->getShape();      const ShapeType& VShape = temp_V->getShape();
812  #pragma omp parallel for  #pragma omp parallel for
813      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 568  int DataExpanded::matrixInverse(DataAbst Line 832  int DataExpanded::matrixInverse(DataAbst
832    
833      const int numdpps=getNumDPPSample();      const int numdpps=getNumDPPSample();
834      const int numSamples = getNumSamples();      const int numSamples = getNumSamples();
835      const ValueType& vec=m_data.getData();      const DataTypes::RealVectorType& vec=m_data_r;
836      int errcode=0;      int errcode=0;
837  #pragma omp parallel  #pragma omp parallel
838      {      {
# Line 578  int DataExpanded::matrixInverse(DataAbst Line 842  int DataExpanded::matrixInverse(DataAbst
842          for (int sampleNo = 0; sampleNo < numSamples; sampleNo++)          for (int sampleNo = 0; sampleNo < numSamples; sampleNo++)
843          {          {
844              // not sure I like all those virtual calls to getPointOffset              // not sure I like all those virtual calls to getPointOffset
845              DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);              DataTypes::RealVectorType::size_type offset=getPointOffset(sampleNo,0);
846              int res=DataMaths::matrix_inverse(vec, getShape(), offset,              int res=DataMaths::matrix_inverse(vec, getShape(), offset,
847                      temp->getVectorRW(), temp->getShape(), offset, numdpps, h);                      temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
848              if (res > errorcode) {              if (res > errorcode) {
# Line 603  void DataExpanded::setToZero() Line 867  void DataExpanded::setToZero()
867      CHECK_FOR_EX_WRITE;      CHECK_FOR_EX_WRITE;
868      const int numSamples = getNumSamples();      const int numSamples = getNumSamples();
869      const int numDataPointsPerSample = getNumDPPSample();      const int numDataPointsPerSample = getNumDPPSample();
870      const DataTypes::ValueType::size_type n = getNoValues();      const DataTypes::RealVectorType::size_type n = getNoValues();
871  #pragma omp parallel for  #pragma omp parallel for
872      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
873          for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {          for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
874              double* p = &m_data[getPointOffset(sampleNo,dataPointNo)];              double* p = &m_data_r[getPointOffset(sampleNo,dataPointNo)];
875              for (int i=0; i<n; ++i)              for (int i=0; i<n; ++i)
876                  p[i] = 0.;                  p[i] = 0.;
877          }          }
# Line 624  void DataExpanded::dump(const std::strin Line 888  void DataExpanded::dump(const std::strin
888      int type=  getFunctionSpace().getTypeCode();      int type=  getFunctionSpace().getTypeCode();
889      int ndims =0;      int ndims =0;
890      long dims[ldims];      long dims[ldims];
891      const double* d_ptr=&(m_data[0]);      const double* d_ptr=&(m_data_r[0]);
892      const DataTypes::ShapeType& shape = getShape();      const DataTypes::ShapeType& shape = getShape();
893      int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();      int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
894      int mpi_num=getFunctionSpace().getDomain()->getMPISize();      int mpi_num=getFunctionSpace().getDomain()->getMPISize();
# Line 688  void DataExpanded::dump(const std::strin Line 952  void DataExpanded::dump(const std::strin
952    
953  void DataExpanded::setTaggedValue(int tagKey,  void DataExpanded::setTaggedValue(int tagKey,
954                                    const DataTypes::ShapeType& pointshape,                                    const DataTypes::ShapeType& pointshape,
955                                    const DataTypes::ValueType& value,                                    const DataTypes::RealVectorType& value,
956                                      int dataOffset)
957    {
958        CHECK_FOR_EX_WRITE;
959        if (isComplex())
960        {
961            CplxVectorType tv;
962        fillComplexFromReal(value, tv);
963        setTaggedValue(tagKey, pointshape, tv, dataOffset);
964            return;
965        }
966        const int numSamples = getNumSamples();
967        const int numDataPointsPerSample = getNumDPPSample();
968        const DataTypes::RealVectorType::size_type n = getNoValues();
969        const real_t* in = &value[0+dataOffset];
970    
971        if (value.size() != n)
972            throw DataException("DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
973    
974    #pragma omp parallel for
975        for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
976            if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey) {
977                for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
978                    real_t* p = &m_data_r[getPointOffset(sampleNo,dataPointNo)];
979                    for (int i=0; i<n; ++i)
980                        p[i] = in[i];
981                }
982            }
983        }
984    }
985    
986    
987    void DataExpanded::setTaggedValue(int tagKey,
988                                      const DataTypes::ShapeType& pointshape,
989                                      const DataTypes::CplxVectorType& value,
990                                    int dataOffset)                                    int dataOffset)
991  {  {
992      CHECK_FOR_EX_WRITE;      CHECK_FOR_EX_WRITE;
993        if (!isComplex())
994        {
995        throw DataException("Programming Error - Attempt to set a complex value on a real object.");
996        }
997      const int numSamples = getNumSamples();      const int numSamples = getNumSamples();
998      const int numDataPointsPerSample = getNumDPPSample();      const int numDataPointsPerSample = getNumDPPSample();
999      const DataTypes::ValueType::size_type n = getNoValues();      const DataTypes::CplxVectorType::size_type n = getNoValues();
1000      const double* in = &value[0+dataOffset];      const DataTypes::cplx_t* in = &value[0+dataOffset];
1001    
1002      if (value.size() != n)      if (value.size() != n)
1003          throw DataException("DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");          throw DataException("DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
# Line 704  void DataExpanded::setTaggedValue(int ta Line 1006  void DataExpanded::setTaggedValue(int ta
1006      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
1007          if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey) {          if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey) {
1008              for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {              for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
1009                  double* p = &m_data[getPointOffset(sampleNo,dataPointNo)];                  cplx_t* p = &m_data_c[getPointOffset(sampleNo,dataPointNo)];
1010                  for (int i=0; i<n; ++i)                  for (int i=0; i<n; ++i)
1011                      p[i] = in[i];                      p[i] = in[i];
1012              }              }
# Line 717  void DataExpanded::reorderByReferenceIDs Line 1019  void DataExpanded::reorderByReferenceIDs
1019  {  {
1020      CHECK_FOR_EX_WRITE;      CHECK_FOR_EX_WRITE;
1021      const int numSamples = getNumSamples();      const int numSamples = getNumSamples();
1022      const DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();      const DataTypes::RealVectorType::size_type n = getNoValues() * getNumDPPSample();
1023      FunctionSpace fs=getFunctionSpace();      FunctionSpace fs=getFunctionSpace();
1024    
1025      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
# Line 727  void DataExpanded::reorderByReferenceIDs Line 1029  void DataExpanded::reorderByReferenceIDs
1029              bool matched=false;              bool matched=false;
1030              for (int sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {              for (int sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
1031                  if (id == reference_ids[sampleNo2]) {                  if (id == reference_ids[sampleNo2]) {
1032                      double* p = &m_data[getPointOffset(sampleNo,0)];                      double* p = &m_data_r[getPointOffset(sampleNo,0)];
1033                      double* p2 = &m_data[getPointOffset(sampleNo2,0)];                      double* p2 = &m_data_r[getPointOffset(sampleNo2,0)];
1034                      for (int i=0; i<n; i++) {                      for (int i=0; i<n; i++) {
1035                          const double rtmp=p[i];                          const double rtmp=p[i];
1036                          p[i] = p2[i];                          p[i] = p2[i];
# Line 746  void DataExpanded::reorderByReferenceIDs Line 1048  void DataExpanded::reorderByReferenceIDs
1048      }      }
1049  }  }
1050    
1051  DataTypes::ValueType& DataExpanded::getVectorRW()  DataTypes::RealVectorType& DataExpanded::getVectorRW()
1052    {
1053        CHECK_FOR_EX_WRITE;
1054        return m_data_r;
1055    }
1056    
1057    const DataTypes::RealVectorType& DataExpanded::getVectorRO() const
1058    {
1059        return m_data_r;
1060    }
1061    
1062    DataTypes::CplxVectorType& DataExpanded::getVectorRWC()
1063  {  {
1064      CHECK_FOR_EX_WRITE;      CHECK_FOR_EX_WRITE;
1065      return m_data.getData();      return m_data_c;
1066  }  }
1067    
1068  const DataTypes::ValueType& DataExpanded::getVectorRO() const  const DataTypes::CplxVectorType& DataExpanded::getVectorROC() const
1069  {  {
1070      return m_data.getData();      return m_data_c;
1071  }  }
1072    
1073    DataTypes::RealVectorType& DataExpanded::getTypedVectorRW(DataTypes::real_t dummypar)
1074    {
1075        CHECK_FOR_EX_WRITE;
1076        return m_data_r;
1077    }
1078    
1079    const DataTypes::RealVectorType& DataExpanded::getTypedVectorRO(DataTypes::real_t dummypar) const
1080    {
1081        return m_data_r;
1082    }
1083    
1084    DataTypes::CplxVectorType& DataExpanded::getTypedVectorRW(DataTypes::cplx_t dummypar)
1085    {
1086        CHECK_FOR_EX_WRITE;
1087        return m_data_c;
1088    }
1089    
1090    const DataTypes::CplxVectorType& DataExpanded::getTypedVectorRO(DataTypes::cplx_t dummypar) const
1091    {
1092        return m_data_c;
1093    }
1094    
1095    
1096  //void DataExpanded::randomFill(long seed)  //void DataExpanded::randomFill(long seed)
1097  //{  //{
1098  //    CHECK_FOR_EX_WRITE;  //    CHECK_FOR_EX_WRITE;

Legend:
Removed from v.5962  
changed lines
  Added in v.5963

  ViewVC Help
Powered by ViewVC 1.1.26