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

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

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

branches/arrayview_from_1695_trunk/escript/src/DataTypes.cpp revision 1714 by jfenwick, Thu Aug 21 00:01:55 2008 UTC trunk/escript/src/DataTypes.cpp revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
 /* $Id:$ */  
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5  *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7  *  *
8  *                http://esscc.uq.edu.au  * Primary Business: Queensland, Australia
9  *        Primary Business: Queensland, Australia  * Licensed under the Open Software License version 3.0
10  *  Licensed under the Open Software License version 3.0  * http://www.opensource.org/licenses/osl-3.0.php
 *     http://www.opensource.org/licenses/osl-3.0.php  
11  *  *
12  *******************************************************/  *******************************************************/
13    
14    
15  #include <sstream>  #include <sstream>
16  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
17  #include <boost/python/tuple.hpp>  #include <boost/python/tuple.hpp>
18  #include "DataException.h"  #include "DataException.h"
19  #include "DataTypes.h"  #include "DataTypes.h"
20    
21    namespace {
22  using namespace boost::python;  using namespace boost::python;
23    using namespace escript;
24  namespace escript  using namespace escript::DataTypes;
 {  
 namespace DataTypes  
 {  
   
    int  
    noValues(const ShapeType& shape)  
    {  
       ShapeType::const_iterator i;  
       //  
       // An empty shape vector means rank 0 which contains 1 value  
       int noValues=1;  
       for (i=shape.begin();i!=shape.end();i++) {  
          noValues*=(*i);  
       }  
       return noValues;  
    }  
   
    int  
    noValues(const RegionLoopRangeType& region)  
    {  
       //  
       // An empty region vector means rank 0 which contains 1 value  
       int noValues=1;  
       unsigned int i;  
       for (i=0;i<region.size();i++) {  
          noValues*=region[i].second-region[i].first;  
       }  
       return noValues;  
    }  
   
    std::string  
    shapeToString(const DataTypes::ShapeType& shape)  
    {  
       std::stringstream temp;  
       temp << "(";  
       unsigned int i;  
       for (i=0;i<shape.size();i++) {  
          temp << shape[i];  
          if (i < shape.size()-1) {  
             temp << ",";  
          }  
       }  
       temp << ")";  
       return temp.str();  
    }  
   
25    
26  /*  /*
27    \brief    \brief
# Line 119  namespace DataTypes Line 74  namespace DataTypes
74           throw DataException("Error - lower index must less or equal upper index.");           throw DataException("Error - lower index must less or equal upper index.");
75        return std::pair<int,int>(s0,s1);        return std::pair<int,int>(s0,s1);
76     }     }
77    }
78    
79    
80    using namespace boost::python;
81    
82    namespace escript
83    {
84    namespace DataTypes
85    {
86    
87       int
88       noValues(const ShapeType& shape)
89       {
90          ShapeType::const_iterator i;
91          //
92          // An empty shape vector means rank 0 which contains 1 value
93          int noValues=1;
94          for (i=shape.begin();i!=shape.end();i++) {
95             noValues*=(*i);
96          }
97          return noValues;
98       }
99    
100       int
101       noValues(const RegionLoopRangeType& region)
102       {
103          //
104          // An empty region vector means rank 0 which contains 1 value
105          int noValues=1;
106          unsigned int i;
107          for (i=0;i<region.size();i++) {
108             noValues*=region[i].second-region[i].first;
109          }
110          return noValues;
111       }
112    
113       std::string
114       shapeToString(const DataTypes::ShapeType& shape)
115       {
116          std::stringstream temp;
117          temp << "(";
118          unsigned int i;
119          for (i=0;i<shape.size();i++) {
120             temp << shape[i];
121             if (i < shape.size()-1) {
122                temp << ",";
123             }
124          }
125          temp << ")";
126          return temp.str();
127       }
128    
129    
130    
131    
132    
# Line 173  namespace DataTypes Line 181  namespace DataTypes
181        return result;        return result;
182     }     }
183    
184       DataTypes::RegionLoopRangeType
185       getSliceRegionLoopRange(const DataTypes::RegionType& region)
186       {
187          DataTypes::RegionLoopRangeType region_loop_range(region.size());
188          unsigned int i;
189          for (i=0;i<region.size();i++) {
190             if (region[i].first==region[i].second) {
191                region_loop_range[i].first=region[i].first;
192                region_loop_range[i].second=region[i].second+1;
193             } else {
194                region_loop_range[i].first=region[i].first;
195                region_loop_range[i].second=region[i].second;
196             }
197          }
198          return region_loop_range;
199       }
200    
201    
202       std::string
203       createShapeErrorMessage(const std::string& messagePrefix,
204                                              const DataTypes::ShapeType& other,
205                          const DataTypes::ShapeType& thisShape)
206       {
207          std::stringstream temp;
208          temp << messagePrefix
209               << " This shape: " << shapeToString(thisShape)
210               << " Other shape: " << shapeToString(other);
211          return temp.str();
212       }
213    
214    
215    // Additional slice operations
216    
217       inline
218       bool
219       checkOffset(ValueType::size_type offset, int size, int noval)
220       {
221          return (size >= (offset+noval));
222       }
223    
224    
225       void
226       copySlice(ValueType& left,
227                    const ShapeType& leftShape,
228                    ValueType::size_type thisOffset,
229                                const ValueType& other,
230                    const ShapeType& otherShape,
231                                ValueType::size_type otherOffset,
232                                const RegionLoopRangeType& region)
233       {
234          //
235          // Make sure views are not empty
236    
237          EsysAssert(!left.size()==0,
238                     "Error - left data is empty.");
239          EsysAssert(!other.size()==0,
240                     "Error - other data is empty.");
241    
242          //
243          // Check the view to be sliced from is compatible with the region to be sliced,
244          // and that the region to be sliced is compatible with this view:
245          EsysAssert(checkOffset(thisOffset,left.size(),noValues(leftShape)),
246                     "Error - offset incompatible with this view.");
247          EsysAssert(otherOffset+noValues(leftShape)<=other.size(),
248                     "Error - offset incompatible with other view.");
249    
250          EsysAssert(getRank(otherShape)==region.size(),
251                     "Error - slice not same rank as view to be sliced from.");
252    
253          EsysAssert(noValues(leftShape)==noValues(getResultSliceShape(region)),
254                     "Error - slice shape not compatible shape for this view.");
255    
256          //
257          // copy the values in the specified region of the other view into this view
258    
259          // the following loops cannot be parallelised due to the numCopy counter
260          int numCopy=0;
261    
262          switch (region.size()) {
263          case 0:
264             /* this case should never be encountered,
265             as python will never pass us an empty region.
266             here for completeness only, allows slicing of a scalar */
267    //          (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
268    
269             left[thisOffset+numCopy]=other[otherOffset];
270             numCopy++;
271             break;
272          case 1:
273             for (int i=region[0].first;i<region[0].second;i++) {
274                left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i)];
275                numCopy++;
276             }
277             break;
278          case 2:
279             for (int j=region[1].first;j<region[1].second;j++) {
280                for (int i=region[0].first;i<region[0].second;i++) {
281    /*               (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];*/
282                   left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j)];
283                   numCopy++;
284                }
285             }
286             break;
287          case 3:
288             for (int k=region[2].first;k<region[2].second;k++) {
289                for (int j=region[1].first;j<region[1].second;j++) {
290                   for (int i=region[0].first;i<region[0].second;i++) {
291    //                  (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
292                      left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k)];
293                      numCopy++;
294                   }
295                }
296             }
297             break;
298          case 4:
299             for (int l=region[3].first;l<region[3].second;l++) {
300                for (int k=region[2].first;k<region[2].second;k++) {
301                   for (int j=region[1].first;j<region[1].second;j++) {
302                      for (int i=region[0].first;i<region[0].second;i++) {
303    /*                     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];*/
304                         left[thisOffset+numCopy]=other[otherOffset+getRelIndex(otherShape,i,j,k,l)];
305                         numCopy++;
306                      }
307                   }
308                }
309             }
310             break;
311          default:
312             std::stringstream mess;
313             mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
314             throw DataException(mess.str());
315          }
316       }
317    
318    
319       void
320       copySliceFrom(ValueType& left,
321                    const ShapeType& leftShape,
322                    ValueType::size_type thisOffset,
323                                    const ValueType& other,
324                    const ShapeType& otherShape,
325                                    ValueType::size_type otherOffset,
326                                    const RegionLoopRangeType& region)
327       {
328          //
329          // Make sure views are not empty
330    
331          EsysAssert(left.size()!=0,
332                     "Error - this view is empty.");
333          EsysAssert(other.size()!=0,
334                     "Error - other view is empty.");
335    
336          //
337          // Check this view is compatible with the region to be sliced,
338          // and that the region to be sliced is compatible with the other view:
339    
340          EsysAssert(checkOffset(otherOffset,other.size(),noValues(otherShape)),
341                     "Error - offset incompatible with other view.");
342          EsysAssert(thisOffset+noValues(otherShape)<=left.size(),
343                     "Error - offset incompatible with this view.");
344    
345          EsysAssert(getRank(leftShape)==region.size(),
346                     "Error - slice not same rank as this view.");
347    
348          EsysAssert(getRank(otherShape)==0 || noValues(otherShape)==noValues(getResultSliceShape(region)),
349                     "Error - slice shape not compatible shape for other view.");
350    
351          //
352          // copy the values in the other view into the specified region of this view
353    
354          // allow for case where other view is a scalar
355          if (getRank(otherShape)==0) {
356    
357             // the following loops cannot be parallelised due to the numCopy counter
358             int numCopy=0;
359    
360             switch (region.size()) {
361             case 0:
362                /* this case should never be encountered,
363                as python will never pass us an empty region.
364                here for completeness only, allows slicing of a scalar */
365                //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
366            left[thisOffset]=other[otherOffset];
367                numCopy++;
368                break;
369             case 1:
370                for (int i=region[0].first;i<region[0].second;i++) {
371                   left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset];
372                   numCopy++;
373                }
374                break;
375             case 2:
376                for (int j=region[1].first;j<region[1].second;j++) {
377                   for (int i=region[0].first;i<region[0].second;i++) {
378                      left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset];
379                      numCopy++;
380                   }
381                }
382                break;
383             case 3:
384                for (int k=region[2].first;k<region[2].second;k++) {
385                   for (int j=region[1].first;j<region[1].second;j++) {
386                      for (int i=region[0].first;i<region[0].second;i++) {
387                         left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset];
388                         numCopy++;
389                      }
390                   }
391                }
392                break;
393             case 4:
394                for (int l=region[3].first;l<region[3].second;l++) {
395                   for (int k=region[2].first;k<region[2].second;k++) {
396                      for (int j=region[1].first;j<region[1].second;j++) {
397                         for (int i=region[0].first;i<region[0].second;i++) {
398                            left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset];
399                            numCopy++;
400                         }
401                      }
402                   }
403                }
404                break;
405             default:
406                std::stringstream mess;
407                mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
408                throw DataException(mess.str());
409             }
410    
411          } else {
412    
413             // the following loops cannot be parallelised due to the numCopy counter
414             int numCopy=0;
415    
416             switch (region.size()) {
417             case 0:
418                /* this case should never be encountered,
419                as python will never pass us an empty region.
420                here for completeness only, allows slicing of a scalar */
421                //(*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
422            left[thisOffset]=other[otherOffset+numCopy];
423                numCopy++;
424                break;
425             case 1:
426                for (int i=region[0].first;i<region[0].second;i++) {
427                   left[thisOffset+getRelIndex(leftShape,i)]=other[otherOffset+numCopy];
428                   numCopy++;
429                }
430                break;
431             case 2:
432                for (int j=region[1].first;j<region[1].second;j++) {
433                   for (int i=region[0].first;i<region[0].second;i++) {
434                      left[thisOffset+getRelIndex(leftShape,i,j)]=other[otherOffset+numCopy];
435                      numCopy++;
436                   }
437                }
438                break;
439             case 3:
440                for (int k=region[2].first;k<region[2].second;k++) {
441                   for (int j=region[1].first;j<region[1].second;j++) {
442                      for (int i=region[0].first;i<region[0].second;i++) {
443                         left[thisOffset+getRelIndex(leftShape,i,j,k)]=other[otherOffset+numCopy];
444                         numCopy++;
445                      }
446                   }
447                }
448                break;
449             case 4:
450                for (int l=region[3].first;l<region[3].second;l++) {
451                   for (int k=region[2].first;k<region[2].second;k++) {
452                      for (int j=region[1].first;j<region[1].second;j++) {
453                         for (int i=region[0].first;i<region[0].second;i++) {
454                            left[thisOffset+getRelIndex(leftShape,i,j,k,l)]=other[otherOffset+numCopy];
455                            numCopy++;
456                         }
457                      }
458                   }
459                }
460                break;
461             default:
462                std::stringstream mess;
463                mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
464                throw DataException(mess.str());
465             }
466    
467          }
468    
469       }
470    
471       std::string
472       pointToString(const ValueType& data,const ShapeType& shape, int offset, const std::string& prefix)
473       {
474          using namespace std;
475          EsysAssert(data.size()>0,"Error - Data object is empty.");
476          stringstream temp;
477          string finalPrefix=prefix;
478          if (prefix.length() > 0) {
479             finalPrefix+=" ";
480          }
481          switch (getRank(shape)) {
482          case 0:
483             temp << finalPrefix << data[offset];
484             break;
485          case 1:
486             for (int i=0;i<shape[0];i++) {
487                temp << finalPrefix << "(" << i <<  ") " << data[i+offset];
488                if (i!=(shape[0]-1)) {
489                   temp << endl;
490                }
491             }
492             break;
493          case 2:
494             for (int i=0;i<shape[0];i++) {
495                for (int j=0;j<shape[1];j++) {
496                   temp << finalPrefix << "(" << i << "," << j << ") " << data[offset+getRelIndex(shape,i,j)];
497                   if (!(i==(shape[0]-1) && j==(shape[1]-1))) {
498                      temp << endl;
499                   }
500                }
501             }
502             break;
503          case 3:
504             for (int i=0;i<shape[0];i++) {
505                for (int j=0;j<shape[1];j++) {
506                   for (int k=0;k<shape[2];k++) {
507                      temp << finalPrefix << "(" << i << "," << j << "," << k << ") " << data[offset+getRelIndex(shape,i,j,k)];
508                      if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1))) {
509                         temp << endl;
510                      }
511                   }
512                }
513             }
514             break;
515          case 4:
516             for (int i=0;i<shape[0];i++) {
517                for (int j=0;j<shape[1];j++) {
518                   for (int k=0;k<shape[2];k++) {
519                      for (int l=0;l<shape[3];l++) {
520                         temp << finalPrefix << "(" << i << "," << j << "," << k << "," << l << ") " << data[offset+getRelIndex(shape,i,j,k,l)];
521                         if (!(i==(shape[0]-1) && j==(shape[1]-1) && k==(shape[2]-1) && l==(shape[3]-1))) {
522                            temp << endl;
523                         }
524                      }
525                   }
526                }
527             }
528             break;
529          default:
530             stringstream mess;
531             mess << "Error - (toString) Invalid rank: " << getRank(shape);
532             throw DataException(mess.str());
533          }
534          return temp.str();
535       }
536    
537    
538       void copyPoint(ValueType& dest, ValueType::size_type doffset, ValueType::size_type nvals, const ValueType& src, ValueType::size_type soffset)
539       {
540          EsysAssert((dest.size()>0&&src.size()>0&&checkOffset(doffset,dest.size(),nvals)),
541                     "Error - Couldn't copy due to insufficient storage.");
542    //       EsysAssert((checkShape(other.getShape())),
543    //                  createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape(),m_shape));
544          if (checkOffset(doffset,dest.size(),nvals) && checkOffset(soffset,src.size(),nvals)) {
545             memcpy(&dest[doffset],&src[soffset],sizeof(double)*nvals);
546          } else {
547             throw DataException("Error - invalid offset specified.");
548          }
549    
550    
551    
552       }
553    
554  }   // end namespace DataTypes  }   // end namespace DataTypes
 }   // end namespace escript  
555    }   // end namespace escript

Legend:
Removed from v.1714  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26