/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 4599 by sshaw, Sun Dec 15 23:51:13 2013 UTC revision 4615 by caltinay, Mon Jan 13 05:05:33 2014 UTC
# Line 169  bool Rectangle::operator==(const Abstrac Line 169  bool Rectangle::operator==(const Abstrac
169  }  }
170    
171  void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,  void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
172              const vector<int>& first, const vector<int>& numValues,              const GridParameters& params) const
             const vector<int>& multiplier) const  
173  {  {
174  #ifdef USE_NETCDF  #ifdef USE_NETCDF
175      // check destination function space      // check destination function space
# Line 185  void Rectangle::readNcGrid(escript::Data Line 184  void Rectangle::readNcGrid(escript::Data
184      } else      } else
185          throw RipleyException("readNcGrid(): invalid function space for output data object");          throw RipleyException("readNcGrid(): invalid function space for output data object");
186    
187      if (first.size() != 2)      if (params.first.size() != 2)
188          throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
189    
190      if (numValues.size() != 2)      if (params.numValues.size() != 2)
191          throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
192    
193      if (multiplier.size() != 2)      if (params.multiplier.size() != 2)
194          throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");          throw RipleyException("readNcGrid(): argument 'multiplier' must have 2 entries");
195      for (size_t i=0; i<multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
196          if (multiplier[i]<1)          if (params.multiplier[i]<1)
197              throw RipleyException("readNcGrid(): all multipliers must be positive");              throw RipleyException("readNcGrid(): all multipliers must be positive");
198    
199      // check file existence and size      // check file existence and size
# Line 216  void Rectangle::readNcGrid(escript::Data Line 215  void Rectangle::readNcGrid(escript::Data
215    
216      // is this a slice of the data object (dims!=2)?      // is this a slice of the data object (dims!=2)?
217      // note the expected ordering of edges (as in numpy: y,x)      // note the expected ordering of edges (as in numpy: y,x)
218      if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))      if ( (dims==2 && (params.numValues[1] > edges[0] || params.numValues[0] > edges[1]))
219              || (dims==1 && numValues[1]>1) ) {              || (dims==1 && params.numValues[1]>1) ) {
220          throw RipleyException("readNcGrid(): not enough data in file");          throw RipleyException("readNcGrid(): not enough data in file");
221      }      }
222    
223      // check if this rank contributes anything      // check if this rank contributes anything
224      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
225              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1])              params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
226                params.first[1] >= m_offset[1]+myN1 ||
227                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1])
228          return;          return;
229    
230      // now determine how much this rank has to write      // now determine how much this rank has to write
231    
232      // first coordinates in data object to write to      // first coordinates in data object to write to
233      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
234      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
235      // indices to first value in file      // indices to first value in file
236      const int idx0 = max(0, m_offset[0]-first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
237      const int idx1 = max(0, m_offset[1]-first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
238      // number of values to read      // number of values to read
239      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
240      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
241    
242      vector<double> values(num0*num1);      vector<double> values(num0*num1);
243      if (dims==2) {      if (dims==2) {
# Line 253  void Rectangle::readNcGrid(escript::Data Line 254  void Rectangle::readNcGrid(escript::Data
254      for (index_t y=0; y<num1; y++) {      for (index_t y=0; y<num1; y++) {
255  #pragma omp parallel for  #pragma omp parallel for
256          for (index_t x=0; x<num0; x++) {          for (index_t x=0; x<num0; x++) {
257              const int baseIndex = first0+x*multiplier[0]              const int baseIndex = first0+x*params.multiplier[0]
258                                    +(first1+y*multiplier[1])*myN0;                                    +(first1+y*params.multiplier[1])*myN0;
259              const int srcIndex = y*num0+x;              const int srcIndex = y*num0+x;
260              if (!isnan(values[srcIndex])) {              if (!isnan(values[srcIndex])) {
261                  for (index_t m1=0; m1<multiplier[1]; m1++) {                  for (index_t m1=0; m1<params.multiplier[1]; m1++) {
262                      for (index_t m0=0; m0<multiplier[0]; m0++) {                      for (index_t m0=0; m0<params.multiplier[0]; m0++) {
263                          const int dataIndex = baseIndex+m0+m1*myN0;                          const int dataIndex = baseIndex+m0+m1*myN0;
264                          double* dest = out.getSampleDataRW(dataIndex);                          double* dest = out.getSampleDataRW(dataIndex);
265                          for (index_t q=0; q<dpp; q++) {                          for (index_t q=0; q<dpp; q++) {
# Line 275  void Rectangle::readNcGrid(escript::Data Line 276  void Rectangle::readNcGrid(escript::Data
276  }  }
277    
278  void Rectangle::readBinaryGrid(escript::Data& out, string filename,  void Rectangle::readBinaryGrid(escript::Data& out, string filename,
279                                 const vector<int>& first,                                 const GridParameters& params) const
                                const vector<int>& numValues,  
                                const std::vector<int>& multiplier,  
                                int byteOrder, int dataType) const  
280  {  {
281      // the mapping is not universally correct but should work on our      // the mapping is not universally correct but should work on our
282      // supported platforms      // supported platforms
283      switch (dataType) {      switch (params.dataType) {
284          case DATATYPE_INT32:          case DATATYPE_INT32:
285              readBinaryGridImpl<int>(out, filename, first, numValues,              readBinaryGridImpl<int>(out, filename, params);
                                     multiplier, byteOrder);  
286              break;              break;
287          case DATATYPE_FLOAT32:          case DATATYPE_FLOAT32:
288              readBinaryGridImpl<float>(out, filename, first, numValues,              readBinaryGridImpl<float>(out, filename, params);
                                       multiplier, byteOrder);  
289              break;              break;
290          case DATATYPE_FLOAT64:          case DATATYPE_FLOAT64:
291              readBinaryGridImpl<double>(out, filename, first, numValues,              readBinaryGridImpl<double>(out, filename, params);
                                        multiplier, byteOrder);  
292              break;              break;
293          default:          default:
294              throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");              throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
# Line 302  void Rectangle::readBinaryGrid(escript:: Line 297  void Rectangle::readBinaryGrid(escript::
297    
298  template<typename ValueType>  template<typename ValueType>
299  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
300                                     const vector<int>& first,                                     const GridParameters& params) const
                                    const vector<int>& numValues,  
                                    const std::vector<int>& multiplier,  
                                    int byteOrder) const  
301  {  {
302      // check destination function space      // check destination function space
303      int myN0, myN1;      int myN0, myN1;
# Line 327  void Rectangle::readBinaryGridImpl(escri Line 319  void Rectangle::readBinaryGridImpl(escri
319      f.seekg(0, ios::end);      f.seekg(0, ios::end);
320      const int numComp = out.getDataPointSize();      const int numComp = out.getDataPointSize();
321      const int filesize = f.tellg();      const int filesize = f.tellg();
322      const int reqsize = numValues[0]*numValues[1]*numComp*sizeof(ValueType);      const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
323      if (filesize < reqsize) {      if (filesize < reqsize) {
324          f.close();          f.close();
325          throw RipleyException("readBinaryGrid(): not enough data in file");          throw RipleyException("readBinaryGrid(): not enough data in file");
326      }      }
327    
328      // check if this rank contributes anything      // check if this rank contributes anything
329      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
330              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1] <= m_offset[1]) {              params.first[0]+params.numValues[0] <= m_offset[0] ||
331                params.first[1] >= m_offset[1]+myN1 ||
332                params.first[1]+params.numValues[1] <= m_offset[1]) {
333          f.close();          f.close();
334          return;          return;
335      }      }
# Line 343  void Rectangle::readBinaryGridImpl(escri Line 337  void Rectangle::readBinaryGridImpl(escri
337      // now determine how much this rank has to write      // now determine how much this rank has to write
338    
339      // first coordinates in data object to write to      // first coordinates in data object to write to
340      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
341      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
342      // indices to first value in file      // indices to first value in file
343      const int idx0 = max(0, m_offset[0]-first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
344      const int idx1 = max(0, m_offset[1]-first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
345      // number of values to read      // number of values to read
346      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
347      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
348    
349      out.requireWrite();      out.requireWrite();
350      vector<ValueType> values(num0*numComp);      vector<ValueType> values(num0*numComp);
351      const int dpp = out.getNumDataPointsPerSample();      const int dpp = out.getNumDataPointsPerSample();
352    
353      for (int y=0; y<num1; y++) {      for (int y=0; y<num1; y++) {
354          const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]);          const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
355          f.seekg(fileofs*sizeof(ValueType));          f.seekg(fileofs*sizeof(ValueType));
356          f.read((char*)&values[0], num0*numComp*sizeof(ValueType));          f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
357          for (int x=0; x<num0; x++) {          for (int x=0; x<num0; x++) {
358              const int baseIndex = first0+x*multiplier[0]              const int baseIndex = first0+x*params.multiplier[0]
359                                      +(first1+y*multiplier[1])*myN0;                                      +(first1+y*params.multiplier[1])*myN0;
360              for (int m1=0; m1<multiplier[1]; m1++) {              for (int m1=0; m1<params.multiplier[1]; m1++) {
361                  for (int m0=0; m0<multiplier[0]; m0++) {                  for (int m0=0; m0<params.multiplier[0]; m0++) {
362                      const int dataIndex = baseIndex+m0+m1*myN0;                      const int dataIndex = baseIndex+m0+m1*myN0;
363                      double* dest = out.getSampleDataRW(dataIndex);                      double* dest = out.getSampleDataRW(dataIndex);
364                      for (int c=0; c<numComp; c++) {                      for (int c=0; c<numComp; c++) {
365                          ValueType val = values[x*numComp+c];                          ValueType val = values[x*numComp+c];
366    
367                          if (byteOrder != BYTEORDER_NATIVE) {                          if (params.byteOrder != BYTEORDER_NATIVE) {
368                              char* cval = reinterpret_cast<char*>(&val);                              char* cval = reinterpret_cast<char*>(&val);
369                              // this will alter val!!                              // this will alter val!!
370                              byte_swap32(cval);                              byte_swap32(cval);

Legend:
Removed from v.4599  
changed lines
  Added in v.4615

  ViewVC Help
Powered by ViewVC 1.1.26