/[escript]/trunk/escript/src/DataArrayView.h
ViewVC logotype

Diff of /trunk/escript/src/DataArrayView.h

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

trunk/esys2/escript/src/Data/DataArrayView.h revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC trunk/escript/src/DataArrayView.h revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC
# Line 1  Line 1 
1  /*  
2   ******************************************************************************  /* $Id$ */
3   *                                                                            *  
4   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  /*******************************************************
5   *                                                                            *   *
6   * This software is the property of ACcESS. No part of this code              *   *           Copyright 2003-2007 by ACceSS MNRF
7   * may be copied in any form or by any means without the expressed written    *   *       Copyright 2007 by University of Queensland
8   * consent of ACcESS.  Copying, use or modification of this software          *   *
9   * by any unauthorised person is illegal unless that person has a software    *   *                http://esscc.uq.edu.au
10   * license agreement with ACcESS.                                             *   *        Primary Business: Queensland, Australia
11   *                                                                            *   *  Licensed under the Open Software License version 3.0
12   ******************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14                                                                               *******************************************************/
15    
16  #if !defined escript_DataArrayView_20040323_H  #if !defined escript_DataArrayView_20040323_H
17  #define escript_DataArrayView_20040323_H  #define escript_DataArrayView_20040323_H
18    #include "system_dep.h"
19    
20  #include "esysUtils/EsysAssert.h"  #include "esysUtils/EsysAssert.h"
21    
22  #include <vector>  #include "DataException.h"
23    #include "DataVector.h"
24    #include "LocalOps.h"
25    
26  #include <boost/python/numeric.hpp>  #include <boost/python/numeric.hpp>
27  #include <boost/python/object.hpp>  #include <boost/python/object.hpp>
28  #include <boost/shared_ptr.hpp>  
29    #include <vector>
30    
31  namespace escript {  namespace escript {
32    
33  /**  /**
34     \brief     \brief
35     DataArrayView provides a view of data allocated externally.     DataArrayView provides a view of external data, configured according
36       to the given shape information and offset.
37    
38     Description:     Description:
39     DataArrayView provides a view of data allocated externally. The     DataArrayView provides a view of data allocated externally. The
40     external data should provide sufficient values so that the dimensions     external data should provide sufficient values so that the dimensions
41     specified for the view will be satisfied. The viewer can update     specified for the view shape will be satisfied. The viewer can update
42     values within the external data but cannot resize the external data.     values within the external data but cannot resize the external data.
43    
44       The view provided represents a single n-dimensional data-point
45       comprised of values taken from the given data array, starting at the
46       specified offset and extending for as many values as are necessary to
47       satisfy the given shape. The default offset can be changed, or different
48       offsets specified, in order to provide views of other data-points in
49       the underlying data array.
50  */  */
51    
52  class DataArrayView {  class DataArrayView {
53    
54    friend bool operator==(const DataArrayView& left, const DataArrayView& right);    ESCRIPT_DLL_API friend bool operator==(const DataArrayView& left, const DataArrayView& right);
55    friend bool operator!=(const DataArrayView& left, const DataArrayView& right);    ESCRIPT_DLL_API friend bool operator!=(const DataArrayView& left, const DataArrayView& right);
56    
57   public:   public:
58    
59    typedef std::vector<double> ValueType;    //
60    typedef std::vector<int> ShapeType;    // Some basic types which define the data values and view shapes.
61      typedef DataVector                        ValueType;
62      //typedef std::vector<double>               ValueType;
63      typedef std::vector<int>                  ShapeType;
64    typedef std::vector<std::pair<int, int> > RegionType;    typedef std::vector<std::pair<int, int> > RegionType;
65      typedef std::vector<std::pair<int, int> > RegionLoopRangeType;
66      static const int maxRank=4;
67    
68    /**    /**
69       \brief       \brief
70       Default constructor for DataArrayView.       Default constructor for DataArrayView.
71    
72       Description:       Description:
73       Default constructor for DataArrayView. Creates an       Default constructor for DataArrayView.
74       empty view with no associated data.       Creates an empty view with no associated data array.
75    
76         This is essentially useless but here for completeness.
77    */    */
78      ESCRIPT_DLL_API
79    DataArrayView();    DataArrayView();
80    
81    /**    /**
# Line 63  class DataArrayView { Line 85  class DataArrayView {
85       Description:       Description:
86       Constructor for DataArrayView.       Constructor for DataArrayView.
87    
88       \param data - Input - Container holding data to be viewed. This must       \param data - Input -
89                             be large enough for the specified shape.                  Array holding data to be viewed. This must be large enough
90       \param viewShape - Input - The shape of the view.                  to supply sufficient values for the specified shape starting at
91       \param offset - Input - Offset into the data at which view should start.                  the given offset.
92         \param viewShape - Input -
93                    The shape of the view to create.
94         \param offset - Input -
95                    Offset into the data at which view should start.
96    */    */
97      ESCRIPT_DLL_API
98    DataArrayView(ValueType& data,    DataArrayView(ValueType& data,
99                  const ShapeType& viewShape,                  const ShapeType& viewShape,
100                  int offset=0);                  int offset=0);
# Line 79  class DataArrayView { Line 106  class DataArrayView {
106       Description:       Description:
107       Copy constructor for DataArrayView.       Copy constructor for DataArrayView.
108    
109       \param other - Input - DataArrayView to copy.       \param other - Input -
110                    DataArrayView to copy.
111    
112       NOTE: The copy references the same data container.       NOTE: The copy references the same data array.
113    */    */
114      ESCRIPT_DLL_API
115    DataArrayView(const DataArrayView& other);    DataArrayView(const DataArrayView& other);
116    
117    /**    /**
118       \brief       \brief
119       Check if view is empty.       Copy from a numarray into the data array viewed by this.
120    */       This must have same shape as the given value - will throw if not.
   bool  
   isEmpty() const;  
   
   /**  
      \brief  
      Copy from a numarray into the data managed by DataArrayView.  
121    */    */
122      ESCRIPT_DLL_API
123    void    void
124    copy(const boost::python::numeric::array& value);    copy(const boost::python::numeric::array& value);
125    
126    /**    /**
127       \brief       \brief
128       Copy from another DataArrayViewinto the data managed by this       Copy data from another DataArrayView into the data array viewed by this
129       DataArrayView at the default offset.       starting at the default offset in both views.
130         The shapes of both views must be the same - will throw if not.
131         NB: views may or may not reference same underlying data array!
132    */    */
133      ESCRIPT_DLL_API
134    void    void
135    copy(const DataArrayView& other);    copy(const DataArrayView& other);
136    
137    /**    /**
138       \brief       \brief
139       Copy from another DataArrayView into this view at the       Copy data from another DataArrayView into this view starting at the
140       given offset.       given offset in this view and the default offset in the other view.
141         The shapes of both views must be the same - will throw if not.
142         NB: views may or may not reference same underlying data array!
143    */    */
144      ESCRIPT_DLL_API
145    void    void
146    copy(ValueType::size_type offset,    copy(ValueType::size_type offset,
147         const DataArrayView& other);         const DataArrayView& other);
148    
149    /**    /**
150       \brief       \brief
151       Copy from another DataArrayView into this view at the       Copy data from another DataArrayView into this view starting at the
152       given offsets.       given offsets in each view.
153         The shapes of both views must be compatible - will throw if not.
154       \param thisOffset - Input - Offset into this view's data array to copy to.       NB: views may or may not reference same underlying data array!
155       \param other - Input - View for the data to copy.  
156       \param otherOffset - Input - Offset into the other view to copy data from.       \param thisOffset - Input -
157                       Offset into this view's data array to copy to.
158         \param other - Input -
159                       View to copy data from.
160         \param otherOffset - Input -
161                       Offset into the other view's data array to copy from.
162    */    */
163      ESCRIPT_DLL_API
164    void    void
165    copy(ValueType::size_type thisOffset,    copy(ValueType::size_type thisOffset,
166         const DataArrayView& other,         const DataArrayView& other,
# Line 132  class DataArrayView { Line 168  class DataArrayView {
168    
169    /**    /**
170       \brief       \brief
171       Copy the given single value over the entire view.       Copy the given single value over the entire view starting at the given
172         offset in the view's data array.
173    
174       \param offset - Input - Offset into this view's data array to copy to.       \param offset - Input -
175       \param value - Input - Value to copy.                     Offset into this view's data array to copy to.
176         \param value - Input -
177                       Value to copy.
178    */    */
179      ESCRIPT_DLL_API
180    void    void
181    copy(ValueType::size_type offset,    copy(ValueType::size_type offset,
182         ValueType::value_type value);         ValueType::value_type value);
183    
184    /**    /**
185         \brief
186         Check if view is empty. ie: does not point to any actual data.
187      */
188      ESCRIPT_DLL_API
189      bool
190      isEmpty() const;
191    
192      /**
193      \brief      \brief
194      Return this view's offset to the start of data.      Return this view's offset into the viewed data array.
195    */    */
196      ESCRIPT_DLL_API
197    ValueType::size_type    ValueType::size_type
198    getOffset() const;    getOffset() const;
199    
200    /**    /**
201      \brief       \brief
202       Return the rank of the data.       Set the offset into the data array at which this view is to start.
203         This could be used to step through the underlying data array by incrementing
204         the offset by noValues successively. Thus this view would provide a moving
205         window on the underlying data with the given shape.
206    */    */
207    int    ESCRIPT_DLL_API
208    getRank() const;    void
209      setOffset(ValueType::size_type offset);
210    
211    /**    /**
212       \brief       \brief
213       Return the shape of the data.       Increment the offset by the number of values in the shape, if possible. Thus
214         moving the view onto the next data point of the given shape in the underlying
215         data array.
216    */    */
217    const ShapeType&    ESCRIPT_DLL_API
218    getShape() const;    void
219      incrOffset();
220    
221    /**    /**
222       \brief       \brief
223       Calculate the number of values for the given shape.       Check the (given) offset will not result in two few values being available in
224         the underlying data array for this view's shape.
225    */    */
226    static int    ESCRIPT_DLL_API
227    noValues(const ShapeType& shape);    bool
228      checkOffset() const;
229    
230      ESCRIPT_DLL_API
231      bool
232      checkOffset(ValueType::size_type offset) const;
233    
234      /**
235        \brief
236         Return the rank of the shape of this view.
237      */
238      ESCRIPT_DLL_API
239      int
240      getRank() const;
241    
242    /**    /**
243       \brief       \brief
244       Return the number of values for the current view.       Return the number of values for the shape of this view.
245    */    */
246      ESCRIPT_DLL_API
247    int    int
248    noValues() const;    noValues() const;
249    
250    /**    /**
251       \brief       \brief
252       Return true if the shapes are the same.       Calculate the number of values for the given shape or region.
253         This is purely a utility method and has no bearing on this view.
254    */    */
255    bool    ESCRIPT_DLL_API
256    checkShape(const DataArrayView::ShapeType& other) const;    static
257      int
258      noValues(const ShapeType& shape);
259    
260      ESCRIPT_DLL_API
261      static
262      int
263      noValues(const RegionLoopRangeType& region);
264    
265    /**    /**
266       \brief       \brief
267       Return a reference to the underlying data.       Return a reference to the underlying data array.
268    */    */
269      ESCRIPT_DLL_API
270    ValueType&    ValueType&
271    getData() const;    getData() const;
272    
273    /**    /**
274       \brief       \brief
275       Return the value with the given index for the view. This takes into account       Return a reference to the data value with the given
276       the offset. Effectively this returns each value of the view in sequential order.       index in this view. This takes into account the offset.
277    */    */
278      ESCRIPT_DLL_API
279    ValueType::reference    ValueType::reference
280    getData(ValueType::size_type i) const;    getData(ValueType::size_type i) const;
281    
282    /**    /**
283       \brief       \brief
284       Create a shape error message. Normally used when there is a shape       Return the shape of this view.
      mismatch.  
   
      \param messagePrefix - Input - First part of the error message.  
      \param other - Input - The other shape.  
285    */    */
286    std::string    ESCRIPT_DLL_API
287    createShapeErrorMessage(const std::string& messagePrefix,    const
288                            const DataArrayView::ShapeType& other) const;    ShapeType&
289      getShape() const;
290    
291    /**    /**
292       \brief       \brief
293       Set the offset into the array.       Return true if the given shape is the same as this view's shape.
294    */    */
295    void    ESCRIPT_DLL_API
296    setOffset(ValueType::size_type offset);    bool
297      checkShape(const ShapeType& other) const;
298    
299    /**    /**
300       \brief       \brief
301       Return the shape of a data view following the slice specified.       Create a shape error message. Normally used when there is a shape
302         mismatch between this shape and the other shape.
303    
304       \param region - Input - Slice region.       \param messagePrefix - Input -
305                           First part of the error message.
306         \param other - Input -
307                           The other shape.
308    */    */
309    static DataArrayView::ShapeType    ESCRIPT_DLL_API
310    getResultSliceShape(const RegionType& region);    std::string
311      createShapeErrorMessage(const std::string& messagePrefix,
312                              const ShapeType& other) const;
313    
314    /**    /**
315       \brief       \brief
316       Copy a slice from the given view. This view must be the right shape for the slice.       Return the viewed data as a formatted string.
317         Not recommended for very large objects!
318    
319       \param other - Input - Data to copy from.       \param suffix - Input -
320       \param region - Input - Slice region.                         Suffix appended to index display.
321    */    */
322    void    ESCRIPT_DLL_API
323    copySlice(const DataArrayView& other,    std::string
324              const RegionType& region);    toString(const std::string& suffix=std::string("")) const;
325    
326    /**    /**
327       \brief       \brief
328       Copy a slice from the given view. This view must be the right shape for the slice.       Return the given shape as a string.
329         This is purely a utility method and has no bearing on this view.
330    
331       \param thisOffset - Input - use this view offset instead of the default.       \param shape - Input.
      \param other - Input - Data to copy from.  
      \param otherOffset - Input - use this slice offset instead of the default.  
      \param region - Input - Slice region.  
332    */    */
333    void    ESCRIPT_DLL_API
334    copySlice(ValueType::size_type thisOffset,    static
335              const DataArrayView& other,    std::string
336              ValueType::size_type otherOffset,    shapeToString(const ShapeType& shape);
             const RegionType& region);  
   
  /**  
      \brief  
      Copy into a slice from the given view. This view must have the same rank as the slice region.  
337    
338       \param other - Input - Data to copy from.    /**
339       \param region - Input - Slice region.      \brief
340        Return the 1 dimensional index into the data array of the only element
341        in the view, *ignoring the offset*.
342        Assumes a rank 0 view.
343    */    */
344    void    ESCRIPT_DLL_API
345    copySliceFrom(const DataArrayView& other,    ValueType::size_type
346                  const RegionType& region);    relIndex() const;
347    
348    /**    /**
349       \brief      \brief
350       Copy into a slice from the given value. This view must have the same rank as the slice region.      Return the 1 dimensional index into the data array of the element at
351        position i in the view, *ignoring the offset*.
352       \param thisOffset - Input - use this view offset instead of the default.      Assumes a rank 1 view.
      \param other - Input - Data to copy from.  
      \param otherOffset - Input - use this slice offset instead of the default.  
      \param region - Input - Slice region.  
353    */    */
354    void    ESCRIPT_DLL_API
355    copySliceFrom(ValueType::size_type thisOffset,    ValueType::size_type
356                  const DataArrayView& other,    relIndex(ValueType::size_type i) const;
                 ValueType::size_type otherOffset,  
                 const RegionType& region);  
357    
358    /**    /**
359       \brief      \brief
360       Determine the shape of the result array for a matrix multiplication.      Return the 1 dimensional index into the data array of the element at
361        position i,j in the view, *ignoring the offset*.
362        Assumes a rank 2 view.
363    */    */
364    static ShapeType    ESCRIPT_DLL_API
365    determineResultShape(const DataArrayView& left,    ValueType::size_type
366                         const DataArrayView& right);    relIndex(ValueType::size_type i,
367               ValueType::size_type j) const;
368    
369    /**    /**
370       \brief      \brief
371       Returns the range for the slices defined by the key which is a Python      Return the 1 dimensional index into the data array of the element at
372       slice object or a tuple of Python slice object.      position i,j,k in the view, *ignoring the offset*.
373    */      Assumes a rank 3 view.
   DataArrayView::RegionType  
   getSliceRegion(const boost::python::object& key) const;  
   /*  
   DataArrayView::RegionType  
   getSliceRegion2(const boost::python::object& key) const;  
374    */    */
375      ESCRIPT_DLL_API
376    // *******************************************************************    ValueType::size_type
377    // NOTE: The following relIndex functions are a hack. The indexing    relIndex(ValueType::size_type i,
378    // mechanism should be split out from DataArrayView to get the flexability             ValueType::size_type j,
379    // needed.             ValueType::size_type k) const;
380    
381    /**    /**
382      \brief      \brief
383      Return the 1 dimensional index of the element at position i,j,k,m      Return the 1 dimensional index into the data array of the element at
384      ignoring the offset.      position i,j,k,m in the view, *ignoring the offset*.
385        Assumes a rank 4 view.
386    */    */
387      ESCRIPT_DLL_API
388    ValueType::size_type    ValueType::size_type
389    relIndex(ValueType::size_type i,    relIndex(ValueType::size_type i,
390             ValueType::size_type j,             ValueType::size_type j,
# Line 316  class DataArrayView { Line 393  class DataArrayView {
393    
394    /**    /**
395      \brief      \brief
396      Return the 1 dimensional index of the element at position i,j,k      Return the 1 dimensional index into the data array of the only element
397      ignoring the offset.      in the view.
398        Assumes a rank 0 view.
399    */    */
400      ESCRIPT_DLL_API
401    ValueType::size_type    ValueType::size_type
402    relIndex(ValueType::size_type i,    index() const;
            ValueType::size_type j,  
            ValueType::size_type k) const;  
403    
404    /**    /**
405      \brief      \brief
406      Return the 1 dimensional index of the element at position i,j      Return the 1 dimensional index into the data array of the element at
407      ignoring the offset.      position i in the view.
408        Assumes a rank 1 view.
409    */    */
410      ESCRIPT_DLL_API
411    ValueType::size_type    ValueType::size_type
412    relIndex(ValueType::size_type i,    index(ValueType::size_type i) const;
            ValueType::size_type j) const;  
   
   // ********************************************************************  
413    
414    /**    /**
415      \brief      \brief
416      Return the 1 dimensional index of the element at position i,j,k,m.      Return the 1 dimensional index into the data array of the element at
417        position i,j in the view.
418        Assumes a rank 2 view.
419    */    */
420      ESCRIPT_DLL_API
421    ValueType::size_type    ValueType::size_type
422    index(ValueType::size_type i,    index(ValueType::size_type i,
423          ValueType::size_type j,          ValueType::size_type j) const;
         ValueType::size_type k,  
         ValueType::size_type m) const;  
424    
425    /**    /**
426      \brief      \brief
427      Return the 1 dimensional index of the element at position i,j,k.      Return the 1 dimensional index into the data array of the element at
428        position i,j,k in the view.
429        Assumes a rank 3 view.
430    */    */
431      ESCRIPT_DLL_API
432    ValueType::size_type    ValueType::size_type
433    index(ValueType::size_type i,    index(ValueType::size_type i,
434          ValueType::size_type j,          ValueType::size_type j,
# Line 356  class DataArrayView { Line 436  class DataArrayView {
436    
437    /**    /**
438      \brief      \brief
439      Return the 1 dimensional index of the element at position i,j.      Return the 1 dimensional index into the data array of the element at
440        position i,j,k,m in the view.
441        Assumes a rank 4 view.
442    */    */
443      ESCRIPT_DLL_API
444    ValueType::size_type    ValueType::size_type
445    index(ValueType::size_type i,    index(ValueType::size_type i,
446          ValueType::size_type j) const;          ValueType::size_type j,
447            ValueType::size_type k,
448            ValueType::size_type m) const;
449    
450    /**    /**
451      \brief      \brief
452      Return the 1 dimensional index of the element at position i.      Return a reference for the only element in the view.
453        Assumes a rank 0 view.
454    */    */
455    ValueType::size_type    ESCRIPT_DLL_API
456    index(ValueType::size_type i) const;    ValueType::reference
457      operator()();
458    
459      ESCRIPT_DLL_API
460      ValueType::const_reference
461      operator()() const;
462    
463    /**    /**
464      \brief      \brief
465      Return a reference to the element at position i.      Return a reference to the element at position i in the view.
466        Assumes a rank 1 view.
467    */    */
468      ESCRIPT_DLL_API
469    ValueType::reference    ValueType::reference
470    operator()(ValueType::size_type i);    operator()(ValueType::size_type i);
471    
472      ESCRIPT_DLL_API
473    ValueType::const_reference    ValueType::const_reference
474    operator()(ValueType::size_type i) const;    operator()(ValueType::size_type i) const;
475    
476    /**    /**
477      \brief      \brief
478      Return a reference to the element at position i,j.      Return a reference to the element at position i,j in the view.
479        Assumes a rank 2 view.
480    */    */
481      ESCRIPT_DLL_API
482    ValueType::reference    ValueType::reference
483    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
484               ValueType::size_type j);               ValueType::size_type j);
485    
486      ESCRIPT_DLL_API
487    ValueType::const_reference    ValueType::const_reference
488    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
489               ValueType::size_type j) const;               ValueType::size_type j) const;
490    
491    /**    /**
492      \brief      \brief
493      Return a reference to the element at position i,j,k.      Return a reference to the element at position i,j,k in the view.
494        Assumes a rank 3 view.
495    */    */
496      ESCRIPT_DLL_API
497    ValueType::reference    ValueType::reference
498    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
499               ValueType::size_type j,               ValueType::size_type j,
500               ValueType::size_type k);               ValueType::size_type k);
501    
502      ESCRIPT_DLL_API
503    ValueType::const_reference    ValueType::const_reference
504    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
505               ValueType::size_type j,               ValueType::size_type j,
# Line 407  class DataArrayView { Line 507  class DataArrayView {
507    
508   /**   /**
509      \brief      \brief
510      Return a reference to the element at position i,j,k,m.      Return a reference to the element at position i,j,k,m in the view.
511        Assumes a rank 4 view.
512    */    */
513      ESCRIPT_DLL_API
514    ValueType::reference    ValueType::reference
515    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
516               ValueType::size_type j,               ValueType::size_type j,
517               ValueType::size_type k,               ValueType::size_type k,
518               ValueType::size_type m);               ValueType::size_type m);
519    
520      ESCRIPT_DLL_API
521    ValueType::const_reference    ValueType::const_reference
522    operator()(ValueType::size_type i,    operator()(ValueType::size_type i,
523               ValueType::size_type j,               ValueType::size_type j,
# Line 422  class DataArrayView { Line 525  class DataArrayView {
525               ValueType::size_type m) const;               ValueType::size_type m) const;
526    
527    /**    /**
528      \brief       \brief
529      Return a reference for the only element, assuming rank 0.       Determine the shape of the specified slice region.
530         This is purely a utility method and has no bearing on this view.
531    
532         \param region - Input -
533                           Slice region.
534    */    */
535    ValueType::reference    ESCRIPT_DLL_API
536    operator()();    static
537      ShapeType
538      getResultSliceShape(const RegionType& region);
539    
540    ValueType::const_reference    /**
541    operator()() const;       \brief
542         Determine the region specified by the given python slice object.
543    
544         \param key - Input -
545                        python slice object specifying region to be returned.
546    
547         The slice object is a tuple of n python slice specifiers, where
548         n <= the rank of this Data object. Each slice specifier specifies the
549         range of indexes to be sliced from the corresponding dimension. The
550         first specifier corresponds to the first dimension, the second to the
551         second and so on. Where n < the rank, the remaining dimensions are
552         sliced across the full range of their indicies.
553    
554         Each slice specifier is of the form "a:b", which specifies a slice
555         from index a, up to but not including index b. Where index a is ommitted
556         a is assumed to be 0. Where index b is ommitted, b is assumed to be the
557         length of this dimension. Where both are ommitted (eg: ":") the slice is
558         assumed to encompass that entire dimension.
559    
560         Where one of the slice specifiers is a single integer, eg: [1], we
561         want to generate a rank-1 dimension object, as opposed to eg: [1,2]
562         which implies we want to take a rank dimensional object with one
563         dimension of size 1.
564    
565         The return value is a vector of pairs with length equal to the rank of
566         this object. Each pair corresponds to the range of indicies from the
567         corresponding dimension to be sliced from, as specified in the input
568         slice object.
569    
570         Examples:
571    
572           For a rank 1 object of shape(5):
573    
574             getSliceRegion(:)   => < <0,5> >
575             getSliceRegion(2:3) => < <2,3> >
576             getSliceRegion(:3)  => < <0,3> >
577             getSliceRegion(2:)  => < <2,5> >
578    
579           For a rank 2 object of shape(4,5):
580    
581             getSliceRegion(2:3) => < <2,3> <0,5> >
582             getSliceRegion(2)   => < <2,3> <0,5> >
583               NB: but return object requested will have rank 1, shape(5), with
584                   values taken from index 2 of this object's first dimension.
585    
586           For a rank 3 object of shape (2,4,6):
587    
588             getSliceRegion(0:2,0:4,0:6) => < <0,2> <0,4> <0,6> >
589             getSliceRegion(:,:,:)       => < <0,2> <0,4> <0,6> >
590             getSliceRegion(0:1)         => < <0,1> <0,4> <0,6> >
591             getSliceRegion(:1,0:2)      => < <0,1> <0,2> <0,6> >
592    
593      */
594      ESCRIPT_DLL_API
595      RegionType
596      getSliceRegion(const boost::python::object& key) const;
597    
598    /**    /**
599       \brief       \brief
600       Perform the unary operation using the given offset       Copy a data slice specified by the given region from the given view
601       instead of the offset defined within the view.       into this view, using the default offsets in both views.
602    
603         \param other - Input -
604                          View to copy data from.
605         \param region - Input -
606                          Region in other view to copy data from.
607    */    */
608    template <class UnaryFunction>    ESCRIPT_DLL_API
609    void    void
610    unaryOp(ValueType::size_type leftOffset,    copySlice(const DataArrayView& other,
611            UnaryFunction operation);              const RegionLoopRangeType& region);
612    
613    /**    /**
614       \brief       \brief
615       Perform the unary operation using the view's offset.       Copy a data slice specified by the given region and offset from the
616         given view into this view at the given offset.
617    
618         \param thisOffset - Input -
619                          Copy the slice to this offset in this view.
620         \param other - Input -
621                          View to copy data from.
622         \param otherOffset - Input -
623                          Copy the slice from this offset in the given view.
624         \param region - Input -
625                          Region in other view to copy data from.
626    */    */
627    template <class UnaryFunction>    ESCRIPT_DLL_API
628    void    void
629    unaryOp(UnaryFunction operation);    copySlice(ValueType::size_type thisOffset,
630                const DataArrayView& other,
631                ValueType::size_type otherOffset,
632                const RegionLoopRangeType& region);
633    
634    /**    /**
635       \brief       \brief
636       Perform the given operation and return a result.       Copy data into a slice specified by the given region in this view from
637         the given view, using the default offsets in both views.
638    
639         \param other - Input -
640                      View to copy data from.
641         \param region - Input -
642                      Region in this view to copy data to.
643    */    */
644    template <class UnaryFunction>    ESCRIPT_DLL_API
645    double    void
646    algorithm(ValueType::size_type leftOffset,    copySliceFrom(const DataArrayView& other,
647              UnaryFunction operation) const;                  const RegionLoopRangeType& region);
648    
649      /**
650         \brief
651         Copy data into a slice specified by the given region and offset in
652         this view from the given view at the given offset.
653    
654         \param thisOffset - Input -
655                        Copy the slice to this offset in this view.
656         \param other - Input -
657                        View to copy data from.
658         \param otherOffset - Input -
659                        Copy the slice from this offset in the given view.
660         \param region - Input -
661                        Region in this view to copy data to.
662      */
663      ESCRIPT_DLL_API
664      void
665      copySliceFrom(ValueType::size_type thisOffset,
666                    const DataArrayView& other,
667                    ValueType::size_type otherOffset,
668                    const RegionLoopRangeType& region);
669    
670    /**    /**
671       \brief       \brief
672       Perform the given operation and return a result. Use the default offset.       Perform the unary operation on the data point specified by the view's
673         default offset. Applies the specified operation to each value in the data
674         point.
675    
676         Called by escript::unaryOp.
677    
678         \param operation - Input -
679                      Operation to apply. Operation must be a pointer to a function.
680    */    */
681    template <class UnaryFunction>    template <class UnaryFunction>
682    double    void
683    algorithm(UnaryFunction operation) const;    unaryOp(UnaryFunction operation);
684    
685    /**    /**
686       \brief       \brief
687       Perform the binary operation. Version which applies a double value       Perform the unary operation on the data point specified by the given
688       to all values within the view. The given offset is used instead of       offset. Applies the specified operation to each value in the data
689       the default offset specified within the view.       point. Operation must be a pointer to a function.
690    
691         Called by escript::unaryOp.
692    
693         \param offset - Input -
694                      Apply the operation to data starting at this offset in this view.
695         \param operation - Input -
696                      Operation to apply. Must be a pointer to a function.
697    */    */
698    template <class BinaryFunction>    template <class UnaryFunction>
699    void    void
700    binaryOp(ValueType::size_type leftOffset,    unaryOp(ValueType::size_type offset,
701             double right,            UnaryFunction operation);
            BinaryFunction operation);  
702    
703    /**    /**
704       \brief       \brief
705       Perform the binary operation. Version which applies a double value       Perform the binary operation on the data points specified by the default
706       to all values within the view.       offsets in this view and in view "right". Applies the specified operation
707         to corresponding values in both data points. Operation must be a pointer
708         to a function.
709    
710         Called by escript::binaryOp.
711    
712         \param right - Input -
713                      View to act as RHS in given binary operation.
714         \param operation - Input -
715                      Operation to apply. Must be a pointer to a function.
716    */    */
717    template <class BinaryFunction>    template <class BinaryFunction>
718    void    void
719    binaryOp(double right,    binaryOp(const DataArrayView& right,
720             BinaryFunction operation);             BinaryFunction operation);
721    
722    /**    /**
723       \brief       \brief
724       Perform the binary operation. The given offsets override the default       Perform the binary operation on the data points specified by the given
725       offsets specified within the views.       offsets in this view and in view "right". Applies the specified operation
726         to corresponding values in both data points. Operation must be a pointer
727         to a function.
728    
729         Called by escript::binaryOp.
730    
731         \param leftOffset - Input -
732                      Apply the operation to data starting at this offset in this view.
733         \param right - Input -
734                      View to act as RHS in given binary operation.
735         \param rightOffset - Input -
736                      Apply the operation to data starting at this offset in right.
737         \param operation - Input -
738                      Operation to apply. Must be a pointer to a function.
739    */    */
740    template <class BinaryFunction>    template <class BinaryFunction>
741    void    void
# Line 502  class DataArrayView { Line 746  class DataArrayView {
746    
747    /**    /**
748       \brief       \brief
749       Perform the binary operation.       Perform the binary operation on the data point specified by the default
750         offset in this view using the scalar value "right". Applies the specified
751         operation to values in the data point. Operation must be a pointer to
752         a function.
753    
754         Called by escript::binaryOp.
755    
756         \param right - Input -
757                      Value to act as RHS in given binary operation.
758         \param operation - Input -
759                      Operation to apply. Must be a pointer to a function.
760    */    */
761    template <class BinaryFunction>    template <class BinaryFunction>
762    void    void
763    binaryOp(const DataArrayView& right,    binaryOp(double right,
764             BinaryFunction operation);             BinaryFunction operation);
765    
766    /**    /**
767       \brief       \brief
768       Return the data as a string. Not recommended for very large objects.       Perform the binary operation on the data point specified by the given
769       \param suffix - Input - Suffix appended to index display.       offset in this view using the scalar value "right". Applies the specified
770         operation to values in the data point. Operation must be a pointer
771         to a function.
772    
773         Called by escript::binaryOp.
774    
775         \param offset - Input -
776                      Apply the operation to data starting at this offset in this view.
777         \param right - Input -
778                      Value to act as RHS in given binary operation.
779         \param operation - Input -
780                      Operation to apply. Must be a pointer to a function.
781    */    */
782    std::string    template <class BinaryFunction>
783    toString(const std::string& suffix=std::string("")) const;    void
784      binaryOp(ValueType::size_type offset,
785               double right,
786               BinaryFunction operation);
787    
788    /**    /**
789       \brief       \brief
790       Return the given shape as a string.       Perform the given data point reduction operation on the data point
791         specified by the default offset into the view. Reduces all elements of
792         the data point using the given operation, returning the result as a
793         scalar. Operation must be a pointer to a function.
794    
795       \param shape - Input.       Called by escript::algorithm.
796    
797         \param operation - Input -
798                      Operation to apply. Must be a pointer to a function.
799    */    */
800    static std::string    template <class BinaryFunction>
801    shapeToString(const DataArrayView::ShapeType& shape);    double
802      reductionOp(BinaryFunction operation,
803                  double initial_value) const;
804    
805    /**    /**
806       \brief       \brief
807       Perform matrix multiply.       Perform the given data point reduction operation on the data point
808         specified by the given offset into the view. Reduces all elements of
809         the data point using the given operation, returning the result as a
810         scalar. Operation must be a pointer to a function.
811    
812       Description:       Called by escript::algorithm.
813       Perform matrix multiply.  
814         \param offset - Input -
815                      Apply the operation to data starting at this offset in this view.
816         \param operation - Input -
817                      Operation to apply. Must be a pointer to a function.
818      */
819      template <class BinaryFunction>
820      double
821      reductionOp(ValueType::size_type offset,
822                  BinaryFunction operation,
823                  double initial_value) const;
824    
825     /**
826         \brief
827         Perform a matrix multiply of the given views.
828         This is purely a utility method and has no bearing on this view.
829    
830         NB: Only multiplies together the two given views, ie: two data-points,
831         would need to call this over all data-points to multiply the entire
832         Data objects involved.
833    
834       \param left - Input - The left hand side.       \param left - Input - The left hand side.
835       \param right - Input - The right hand side.       \param right - Input - The right hand side.
836       \param result - Output - The result of the operation.       \param result - Output - The result of the operation.
837    */    */
838    static void    ESCRIPT_DLL_API
839      static
840      void
841    matMult(const DataArrayView& left,    matMult(const DataArrayView& left,
842            const DataArrayView& right,            const DataArrayView& right,
843            DataArrayView& result);            DataArrayView& result);
844    
845      /**
846         \brief
847         Determine the shape of the result array for a matrix multiplication
848         of the given views.
849         This is purely a utility method and has no bearing on this view.
850      */
851      ESCRIPT_DLL_API
852      static
853      ShapeType
854      determineResultShape(const DataArrayView& left,
855                           const DataArrayView& right);
856    
857      /**
858         \brief
859         computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
860    
861         \param in - Input - matrix
862         \param inOffset - Input - offset into in
863         \param ev - Output - The symmetric matrix
864         \param inOffset - Input - offset into ev
865      */
866      ESCRIPT_DLL_API
867      static
868      inline
869      void
870      symmetric(DataArrayView& in,
871                ValueType::size_type inOffset,
872                DataArrayView& ev,
873                ValueType::size_type evOffset)
874      {
875       if (in.getRank() == 2) {
876         int i0, i1;
877         int s0=in.getShape()[0];
878         int s1=in.getShape()[1];
879         for (i0=0; i0<s0; i0++) {
880           for (i1=0; i1<s1; i1++) {
881             (*(ev.m_data))[evOffset+ev.index(i0,i1)] = ((*(in.m_data))[inOffset+in.index(i0,i1)] + (*(in.m_data))[inOffset+in.index(i1,i0)]) / 2.0;
882           }
883         }
884       }
885       if (in.getRank() == 4) {
886         int i0, i1, i2, i3;
887         int s0=in.getShape()[0];
888         int s1=in.getShape()[1];
889         int s2=in.getShape()[2];
890         int s3=in.getShape()[3];
891         for (i0=0; i0<s0; i0++) {
892           for (i1=0; i1<s1; i1++) {
893             for (i2=0; i2<s2; i2++) {
894               for (i3=0; i3<s3; i3++) {
895                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = ((*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)] + (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)]) / 2.0;
896               }
897             }
898           }
899         }
900       }
901      }
902    
903      /**
904         \brief
905         computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
906    
907         \param in - Input - matrix
908         \param inOffset - Input - offset into in
909         \param ev - Output - The nonsymmetric matrix
910         \param inOffset - Input - offset into ev
911      */
912      ESCRIPT_DLL_API
913      static
914      inline
915      void
916      nonsymmetric(DataArrayView& in,
917                ValueType::size_type inOffset,
918                DataArrayView& ev,
919                ValueType::size_type evOffset)
920      {
921       if (in.getRank() == 2) {
922         int i0, i1;
923         int s0=in.getShape()[0];
924         int s1=in.getShape()[1];
925         for (i0=0; i0<s0; i0++) {
926           for (i1=0; i1<s1; i1++) {
927             (*(ev.m_data))[evOffset+ev.index(i0,i1)] = ((*(in.m_data))[inOffset+in.index(i0,i1)] - (*(in.m_data))[inOffset+in.index(i1,i0)]) / 2.0;
928           }
929         }
930       }
931       if (in.getRank() == 4) {
932         int i0, i1, i2, i3;
933         int s0=in.getShape()[0];
934         int s1=in.getShape()[1];
935         int s2=in.getShape()[2];
936         int s3=in.getShape()[3];
937         for (i0=0; i0<s0; i0++) {
938           for (i1=0; i1<s1; i1++) {
939             for (i2=0; i2<s2; i2++) {
940               for (i3=0; i3<s3; i3++) {
941                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = ((*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)] - (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)]) / 2.0;
942               }
943             }
944           }
945         }
946       }
947      }
948    
949      /**
950         \brief
951         computes the trace of a matrix
952    
953         \param in - Input - matrix
954         \param inOffset - Input - offset into in
955         \param ev - Output - The nonsymmetric matrix
956         \param inOffset - Input - offset into ev
957      */
958      static
959      inline
960      void
961      trace(DataArrayView& in,
962                ValueType::size_type inOffset,
963                DataArrayView& ev,
964                ValueType::size_type evOffset,
965            int axis_offset)
966      {
967       if (in.getRank() == 2) {
968         int s0=in.getShape()[0]; // Python wrapper limits to square matrix
969         int i;
970         for (i=0; i<s0; i++) {
971           (*(ev.m_data))[evOffset+ev.index()] += (*(in.m_data))[inOffset+in.index(i,i)];
972         }
973       }
974       else if (in.getRank() == 3) {
975         if (axis_offset==0) {
976           int s0=in.getShape()[0];
977           int s2=in.getShape()[2];
978           int i0, i2;
979           for (i0=0; i0<s0; i0++) {
980             for (i2=0; i2<s2; i2++) {
981               (*(ev.m_data))[evOffset+ev.index(i2)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2)];
982             }
983           }
984         }
985         else if (axis_offset==1) {
986           int s0=in.getShape()[0];
987           int s1=in.getShape()[1];
988           int i0, i1;
989           for (i0=0; i0<s0; i0++) {
990             for (i1=0; i1<s1; i1++) {
991               (*(ev.m_data))[evOffset+ev.index(i0)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1)];
992             }
993           }
994         }
995       }
996       else if (in.getRank() == 4) {
997         if (axis_offset==0) {
998           int s0=in.getShape()[0];
999           int s2=in.getShape()[2];
1000           int s3=in.getShape()[3];
1001           int i0, i2, i3;
1002           for (i0=0; i0<s0; i0++) {
1003             for (i2=0; i2<s2; i2++) {
1004               for (i3=0; i3<s3; i3++) {
1005                 (*(ev.m_data))[evOffset+ev.index(i2,i3)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2,i3)];
1006               }
1007             }
1008           }
1009         }
1010         else if (axis_offset==1) {
1011           int s0=in.getShape()[0];
1012           int s1=in.getShape()[1];
1013           int s3=in.getShape()[3];
1014           int i0, i1, i3;
1015           for (i0=0; i0<s0; i0++) {
1016             for (i1=0; i1<s1; i1++) {
1017               for (i3=0; i3<s3; i3++) {
1018                 (*(ev.m_data))[evOffset+ev.index(i0,i3)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1,i3)];
1019               }
1020             }
1021           }
1022         }
1023         else if (axis_offset==2) {
1024           int s0=in.getShape()[0];
1025           int s1=in.getShape()[1];
1026           int s2=in.getShape()[2];
1027           int i0, i1, i2;
1028           for (i0=0; i0<s0; i0++) {
1029             for (i1=0; i1<s1; i1++) {
1030               for (i2=0; i2<s2; i2++) {
1031                 (*(ev.m_data))[evOffset+ev.index(i0,i1)] += (*(in.m_data))[inOffset+in.index(i0,i1,i2,i2)];
1032               }
1033             }
1034           }
1035         }
1036       }
1037      }
1038    
1039      /**
1040         \brief
1041         Transpose each data point of this Data object around the given axis.
1042    
1043         \param in - Input - matrix
1044         \param inOffset - Input - offset into in
1045         \param ev - Output - The nonsymmetric matrix
1046         \param inOffset - Input - offset into ev
1047      */
1048      ESCRIPT_DLL_API
1049      static
1050      inline
1051      void
1052      transpose(DataArrayView& in,
1053                ValueType::size_type inOffset,
1054                DataArrayView& ev,
1055                ValueType::size_type evOffset,
1056            int axis_offset)
1057      {
1058       if (in.getRank() == 4) {
1059         int s0=ev.getShape()[0];
1060         int s1=ev.getShape()[1];
1061         int s2=ev.getShape()[2];
1062         int s3=ev.getShape()[3];
1063         int i0, i1, i2, i3;
1064         if (axis_offset==1) {
1065           for (i0=0; i0<s0; i0++) {
1066             for (i1=0; i1<s1; i1++) {
1067               for (i2=0; i2<s2; i2++) {
1068                 for (i3=0; i3<s3; i3++) {
1069                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i3,i0,i1,i2)];
1070                 }
1071               }
1072             }
1073           }
1074         }
1075         else if (axis_offset==2) {
1076           for (i0=0; i0<s0; i0++) {
1077             for (i1=0; i1<s1; i1++) {
1078               for (i2=0; i2<s2; i2++) {
1079                 for (i3=0; i3<s3; i3++) {
1080                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)];
1081                 }
1082               }
1083             }
1084           }
1085         }
1086         else if (axis_offset==3) {
1087           for (i0=0; i0<s0; i0++) {
1088             for (i1=0; i1<s1; i1++) {
1089               for (i2=0; i2<s2; i2++) {
1090                 for (i3=0; i3<s3; i3++) {
1091                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i2,i3,i0)];
1092                 }
1093               }
1094             }
1095           }
1096         }
1097         else {
1098           for (i0=0; i0<s0; i0++) {
1099             for (i1=0; i1<s1; i1++) {
1100               for (i2=0; i2<s2; i2++) {
1101                 for (i3=0; i3<s3; i3++) {
1102                   (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)];
1103                 }
1104               }
1105             }
1106           }
1107         }
1108       }
1109       else if (in.getRank() == 3) {
1110         int s0=ev.getShape()[0];
1111         int s1=ev.getShape()[1];
1112         int s2=ev.getShape()[2];
1113         int i0, i1, i2;
1114         if (axis_offset==1) {
1115           for (i0=0; i0<s0; i0++) {
1116             for (i1=0; i1<s1; i1++) {
1117               for (i2=0; i2<s2; i2++) {
1118                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i2,i0,i1)];
1119               }
1120             }
1121           }
1122         }
1123         else if (axis_offset==2) {
1124           for (i0=0; i0<s0; i0++) {
1125             for (i1=0; i1<s1; i1++) {
1126               for (i2=0; i2<s2; i2++) {
1127                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i2,i0)];
1128               }
1129             }
1130           }
1131         }
1132         else {
1133           // Copy the matrix unchanged
1134           for (i0=0; i0<s0; i0++) {
1135             for (i1=0; i1<s1; i1++) {
1136               for (i2=0; i2<s2; i2++) {
1137                 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2)];
1138               }
1139             }
1140           }
1141         }
1142       }
1143       else if (in.getRank() == 2) {
1144         int s0=ev.getShape()[0];
1145         int s1=ev.getShape()[1];
1146         int i0, i1;
1147         if (axis_offset==1) {
1148           for (i0=0; i0<s0; i0++) {
1149             for (i1=0; i1<s1; i1++) {
1150               (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1151             }
1152           }
1153         }
1154         else {
1155           for (i0=0; i0<s0; i0++) {
1156             for (i1=0; i1<s1; i1++) {
1157               (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i0,i1)];
1158             }
1159           }
1160         }
1161       }
1162       else if (in.getRank() == 1) {
1163         int s0=ev.getShape()[0];
1164         int i0;
1165         for (i0=0; i0<s0; i0++) {
1166           (*(ev.m_data))[evOffset+ev.index(i0)] = (*(in.m_data))[inOffset+in.index(i0)];
1167         }
1168       }
1169       else if (in.getRank() == 0) {
1170         (*(ev.m_data))[evOffset+ev.index()] = (*(in.m_data))[inOffset+in.index()];
1171       }
1172       else {
1173          throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
1174       }
1175      }
1176    
1177      /**
1178         \brief
1179         swaps the components axis0 and axis1.
1180    
1181         \param in - Input - matrix
1182         \param inOffset - Input - offset into in
1183         \param ev - Output - The nonsymmetric matrix
1184         \param inOffset - Input - offset into ev
1185         \param axis0 - axis index
1186         \param axis1 - axis index
1187      */
1188      ESCRIPT_DLL_API
1189      static
1190      inline
1191      void
1192      swapaxes(DataArrayView& in,
1193               ValueType::size_type inOffset,
1194               DataArrayView& ev,
1195               ValueType::size_type evOffset,
1196               int axis0, int axis1)
1197      {
1198       if (in.getRank() == 4) {
1199         int s0=ev.getShape()[0];
1200         int s1=ev.getShape()[1];
1201         int s2=ev.getShape()[2];
1202         int s3=ev.getShape()[3];
1203         int i0, i1, i2, i3;
1204         if (axis0==0) {
1205            if (axis1==1) {
1206                for (i0=0; i0<s0; i0++) {
1207                  for (i1=0; i1<s1; i1++) {
1208                    for (i2=0; i2<s2; i2++) {
1209                      for (i3=0; i3<s3; i3++) {
1210                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i0,i2,i3)];
1211                      }
1212                    }
1213                  }
1214                }
1215            } else if (axis1==2) {
1216                for (i0=0; i0<s0; i0++) {
1217                  for (i1=0; i1<s1; i1++) {
1218                    for (i2=0; i2<s2; i2++) {
1219                      for (i3=0; i3<s3; i3++) {
1220                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i2,i1,i0,i3)];
1221                      }
1222                    }
1223                  }
1224                }
1225    
1226            } else if (axis1==3) {
1227                for (i0=0; i0<s0; i0++) {
1228                  for (i1=0; i1<s1; i1++) {
1229                    for (i2=0; i2<s2; i2++) {
1230                      for (i3=0; i3<s3; i3++) {
1231                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i3,i1,i2,i0)];
1232                      }
1233                    }
1234                  }
1235                }
1236            }
1237         } else if (axis0==1) {
1238            if (axis1==2) {
1239                for (i0=0; i0<s0; i0++) {
1240                  for (i1=0; i1<s1; i1++) {
1241                    for (i2=0; i2<s2; i2++) {
1242                      for (i3=0; i3<s3; i3++) {
1243                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i2,i1,i3)];
1244                      }
1245                    }
1246                  }
1247                }
1248            } else if (axis1==3) {
1249                for (i0=0; i0<s0; i0++) {
1250                  for (i1=0; i1<s1; i1++) {
1251                    for (i2=0; i2<s2; i2++) {
1252                      for (i3=0; i3<s3; i3++) {
1253                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i3,i2,i1)];
1254                      }
1255                    }
1256                  }
1257                }
1258            }
1259         } else if (axis0==2) {
1260            if (axis1==3) {
1261                for (i0=0; i0<s0; i0++) {
1262                  for (i1=0; i1<s1; i1++) {
1263                    for (i2=0; i2<s2; i2++) {
1264                      for (i3=0; i3<s3; i3++) {
1265                        (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i3,i2)];
1266                      }
1267                    }
1268                  }
1269                }
1270            }
1271         }
1272    
1273       } else if ( in.getRank() == 3) {
1274         int s0=ev.getShape()[0];
1275         int s1=ev.getShape()[1];
1276         int s2=ev.getShape()[2];
1277         int i0, i1, i2;
1278         if (axis0==0) {
1279            if (axis1==1) {
1280               for (i0=0; i0<s0; i0++) {
1281                 for (i1=0; i1<s1; i1++) {
1282                   for (i2=0; i2<s2; i2++) {
1283                     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i0,i2)];
1284                   }
1285                 }
1286               }
1287            } else if (axis1==2) {
1288               for (i0=0; i0<s0; i0++) {
1289                 for (i1=0; i1<s1; i1++) {
1290                   for (i2=0; i2<s2; i2++) {
1291                     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i2,i1,i0)];
1292                   }
1293                 }
1294               }
1295           }
1296         } else if (axis0==1) {
1297            if (axis1==2) {
1298               for (i0=0; i0<s0; i0++) {
1299                 for (i1=0; i1<s1; i1++) {
1300                   for (i2=0; i2<s2; i2++) {
1301                     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i2,i1)];
1302                   }
1303                 }
1304               }
1305            }
1306         }
1307       } else if ( in.getRank() == 2) {
1308         int s0=ev.getShape()[0];
1309         int s1=ev.getShape()[1];
1310         int i0, i1;
1311         if (axis0==0) {
1312            if (axis1==1) {
1313               for (i0=0; i0<s0; i0++) {
1314                 for (i1=0; i1<s1; i1++) {
1315                     (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1316                 }
1317               }
1318            }
1319        }
1320      } else {
1321          throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
1322      }
1323     }
1324    
1325      /**
1326         \brief
1327         solves a local eigenvalue problem
1328    
1329         \param in - Input - matrix
1330         \param inOffset - Input - offset into in
1331         \param ev - Output - The eigenvalues
1332         \param inOffset - Input - offset into ev
1333      */
1334      ESCRIPT_DLL_API
1335      static
1336      inline
1337      void
1338      eigenvalues(DataArrayView& in,
1339                  ValueType::size_type inOffset,
1340                  DataArrayView& ev,
1341                  ValueType::size_type evOffset)
1342      {
1343       double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1344       double ev0,ev1,ev2;
1345       int s=in.getShape()[0];
1346       if (s==1) {
1347          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1348          eigenvalues1(in00,&ev0);
1349          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1350    
1351       } else  if (s==2) {
1352          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1353          in10=(*(in.m_data))[inOffset+in.index(1,0)];
1354          in01=(*(in.m_data))[inOffset+in.index(0,1)];
1355          in11=(*(in.m_data))[inOffset+in.index(1,1)];
1356          eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
1357          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1358          (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1359    
1360       } else  if (s==3) {
1361          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1362          in10=(*(in.m_data))[inOffset+in.index(1,0)];
1363          in20=(*(in.m_data))[inOffset+in.index(2,0)];
1364          in01=(*(in.m_data))[inOffset+in.index(0,1)];
1365          in11=(*(in.m_data))[inOffset+in.index(1,1)];
1366          in21=(*(in.m_data))[inOffset+in.index(2,1)];
1367          in02=(*(in.m_data))[inOffset+in.index(0,2)];
1368          in12=(*(in.m_data))[inOffset+in.index(1,2)];
1369          in22=(*(in.m_data))[inOffset+in.index(2,2)];
1370          eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1371                     &ev0,&ev1,&ev2);
1372          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1373          (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1374          (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1375    
1376       }
1377      }
1378    
1379      /**
1380         \brief
1381         solves a local eigenvalue problem
1382    
1383         \param in - Input - matrix
1384         \param inOffset - Input - offset into in
1385         \param ev - Output - The eigenvalues
1386         \param evOffset - Input - offset into ev
1387         \param V - Output - The eigenvectors
1388         \param VOffset - Input - offset into V
1389         \param tol - Input - eigenvalues with relative difference tol are treated as equal
1390      */
1391      ESCRIPT_DLL_API
1392      static
1393      inline
1394      void
1395      eigenvalues_and_eigenvectors(DataArrayView& in,
1396                                   ValueType::size_type inOffset,
1397                                   DataArrayView& ev,
1398                                   ValueType::size_type evOffset,
1399                                   DataArrayView& V,
1400                                   ValueType::size_type VOffset,
1401                                   const double tol=1.e-13)
1402      {
1403       double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1404       double V00,V10,V20,V01,V11,V21,V02,V12,V22;
1405       double ev0,ev1,ev2;
1406       int s=in.getShape()[0];
1407       if (s==1) {
1408          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1409          eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
1410          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1411          (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1412       } else  if (s==2) {
1413          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1414          in10=(*(in.m_data))[inOffset+in.index(1,0)];
1415          in01=(*(in.m_data))[inOffset+in.index(0,1)];
1416          in11=(*(in.m_data))[inOffset+in.index(1,1)];
1417          eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
1418                       &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
1419          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1420          (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1421          (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1422          (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1423          (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1424          (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1425       } else  if (s==3) {
1426          in00=(*(in.m_data))[inOffset+in.index(0,0)];
1427          in10=(*(in.m_data))[inOffset+in.index(1,0)];
1428          in20=(*(in.m_data))[inOffset+in.index(2,0)];
1429          in01=(*(in.m_data))[inOffset+in.index(0,1)];
1430          in11=(*(in.m_data))[inOffset+in.index(1,1)];
1431          in21=(*(in.m_data))[inOffset+in.index(2,1)];
1432          in02=(*(in.m_data))[inOffset+in.index(0,2)];
1433          in12=(*(in.m_data))[inOffset+in.index(1,2)];
1434          in22=(*(in.m_data))[inOffset+in.index(2,2)];
1435          eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1436                     &ev0,&ev1,&ev2,
1437                     &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
1438          (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1439          (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1440          (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1441          (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1442          (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1443          (*(V.m_data))[inOffset+V.index(2,0)]=V20;
1444          (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1445          (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1446          (*(V.m_data))[inOffset+V.index(2,1)]=V21;
1447          (*(V.m_data))[inOffset+V.index(0,2)]=V02;
1448          (*(V.m_data))[inOffset+V.index(1,2)]=V12;
1449          (*(V.m_data))[inOffset+V.index(2,2)]=V22;
1450    
1451       }
1452     }
1453   protected:   protected:
1454    
1455   private:   private:
1456    
1457    //    //
1458      // The maximum rank allowed for the shape of any view.
1459    static const int m_maxRank=4;    static const int m_maxRank=4;
1460    
1461    //    //
1462    // The data values for the view. NOTE: This points to data external to the view.    // The data values for the view.
1463      // NOTE: This points to data external to the view.
1464      // This is just a pointer to an array of ValueType.
1465    ValueType* m_data;    ValueType* m_data;
1466    
1467    //    //
1468    // The offset into the data array used by different views.    // The offset into the data array used by different views.
1469      // This is simply an integer specifying a position in the data array
1470      // pointed to by m_data.
1471    ValueType::size_type m_offset;    ValueType::size_type m_offset;
1472    
1473    //    //
1474    // The shape of the data.    // The shape of the data.
1475      // This is simply an STL vector specifying the lengths of each dimension
1476      // of the shape as ints.
1477    ShapeType m_shape;    ShapeType m_shape;
1478    
1479    //    //
1480    // The number of values needed for the array.    // The number of values needed for the array.
1481      // This can be derived from m_shape by multiplying the size of each dimension, but
1482      // is stored here for convenience.
1483    int m_noValues;    int m_noValues;
1484    
1485  };  };
1486    
1487  inline  ESCRIPT_DLL_API bool operator==(const DataArrayView& left, const DataArrayView& right);
1488  DataArrayView::ValueType::size_type  ESCRIPT_DLL_API bool operator!=(const DataArrayView& left, const DataArrayView& right);
 DataArrayView::getOffset() const  
 {  
   return m_offset;  
 }  
1489    
1490  inline  /**
1491  DataArrayView::ValueType&    \brief
1492  DataArrayView::getData() const     Modify region to copy from in order to
1493  {     deal with the case where one range in the region contains identical indexes,
1494    EsysAssert(!isEmpty(),"Error - View is empty");     eg: <<1,1><0,3><0,3>>
1495    return *m_data;     This situation implies we want to copy from an object with rank greater than that of this
1496  }     object. eg: we want to copy the values from a two dimensional slice out of a three
1497       dimensional object into a two dimensional object.
1498       We do this by taking a slice from the other object where one dimension of
1499       the slice region is of size 1. So in the above example, we modify the above
1500       region like so: <<1,2><0,3><0,3>> and take this slice.
1501    */
1502    DataArrayView::RegionLoopRangeType
1503    getSliceRegionLoopRange(const DataArrayView::RegionType& region);
1504    
1505  inline  /**
1506  DataArrayView::ValueType::reference    \brief
1507  DataArrayView::getData(ValueType::size_type i) const    Calculate the slice range from the given python key object
1508  {    Used by DataArrayView::getSliceRegion.
1509    //    Returns the python slice object key as a pair of ints where the first
1510    // don't do any checking to allow one past the end of the vector providing    member is the start and the second member is the end. the presence of a possible
1511    // the equivalent of end()    step attribute with value different from one will throw an exception
   return (*m_data)[i+m_offset];  
 }  
1512    
1513  template <class UnaryFunction>    /param key - Input - key object specifying slice range.
1514  inline  */
1515  void  std::pair<int,int>
1516  DataArrayView::unaryOp(ValueType::size_type leftOffset,  getSliceRange(const boost::python::object& key,
1517                         UnaryFunction operation)                const int shape);
1518  {  
1519    for (ValueType::size_type i=0;i<(noValues(m_shape));++i) {  /**
1520      (*m_data)[i+leftOffset]=operation((*m_data)[i+leftOffset]);     Inline function definitions.
1521    }  */
 }  
1522    
1523  template <class UnaryFunction>  template <class UnaryFunction>
1524  inline  inline
# Line 613  DataArrayView::unaryOp(UnaryFunction ope Line 1530  DataArrayView::unaryOp(UnaryFunction ope
1530    
1531  template <class UnaryFunction>  template <class UnaryFunction>
1532  inline  inline
1533  double  void
1534  DataArrayView::algorithm(ValueType::size_type leftOffset,  DataArrayView::unaryOp(ValueType::size_type offset,
1535                           UnaryFunction operation) const                         UnaryFunction operation)
1536  {  {
1537    for (ValueType::size_type i=0;i<(noValues(m_shape));++i) {    EsysAssert((!isEmpty()&&checkOffset(offset)),
1538      operation((*m_data)[i+leftOffset]);                 "Error - Couldn't perform unaryOp due to insufficient storage.");
1539      for (ValueType::size_type i=0;i<noValues();i++) {
1540        (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1541    }    }
   return operation.getResult();  
1542  }  }
1543    
1544  template <class UnaryFunction>  template <class BinaryFunction>
1545  inline  inline
1546  double  void
1547  DataArrayView::algorithm(UnaryFunction operation) const  DataArrayView::binaryOp(const DataArrayView& right,
1548                            BinaryFunction operation)
1549  {  {
1550    return algorithm(m_offset,operation);    binaryOp(m_offset,right,right.getOffset(),operation);
1551  }  }
1552    
1553  template <class BinaryFunction>  template <class BinaryFunction>
1554  inline  inline
1555  void  void
1556  DataArrayView::binaryOp(ValueType::size_type leftOffset,  DataArrayView::binaryOp(ValueType::size_type leftOffset,
1557                          const DataArrayView& right,                          const DataArrayView& right,
1558                          ValueType::size_type rightOffset,                          ValueType::size_type rightOffset,
1559                          BinaryFunction operation)                          BinaryFunction operation)
1560  {  {
1561    EsysAssert(getShape()==right.getShape(),    EsysAssert(getShape()==right.getShape(),
1562           "Error - Right hand shape: " << shapeToString(right.getShape()) << " doesn't match the left: " << shapeToString(getShape()));           "Error - Couldn't perform binaryOp due to shape mismatch,");
1563    for (ValueType::size_type i=0;i<noValues();++i) {    EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1564      (*m_data)[i+leftOffset]=operation((*m_data)[i+leftOffset],(*right.m_data)[i+rightOffset]);               "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1565      EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1566                 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1567      for (ValueType::size_type i=0;i<noValues();i++) {
1568        (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1569    }    }
1570  }  }
1571    
1572  template <class BinaryFunction>  template <class BinaryFunction>
1573  inline  inline
1574  void  void
1575  DataArrayView::binaryOp(const DataArrayView& right,  DataArrayView::binaryOp(double right,
1576                          BinaryFunction operation)                          BinaryFunction operation)
1577  {  {
1578    //    binaryOp(m_offset,right,operation);
   // version using the offsets specified within the view  
   binaryOp(m_offset,right,right.getOffset(),operation);  
1579  }  }
1580    
1581  template <class BinaryFunction>  template <class BinaryFunction>
1582  inline  inline
1583  void  void
1584  DataArrayView::binaryOp(ValueType::size_type leftOffset,  DataArrayView::binaryOp(ValueType::size_type offset,
1585                          double right,                          double right,
1586                          BinaryFunction operation)                          BinaryFunction operation)
1587  {  {
1588    //    EsysAssert((!isEmpty()&&checkOffset(offset)),
1589    // If a scalar is to be applied to the entire array force the caller to               "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1590    // explicitly specify a single value    for (ValueType::size_type i=0;i<noValues();i++) {
1591    for (ValueType::size_type i=0;i<(noValues(m_shape));++i) {      (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
     (*m_data)[i+leftOffset]=operation((*m_data)[i+leftOffset],right);  
1592    }    }
1593  }  }
1594    
1595  template <class BinaryFunction>  template <class BinaryFunction>
1596  inline  inline
1597  void  double
1598  DataArrayView::binaryOp(double right,  DataArrayView::reductionOp(BinaryFunction operation,
1599                          BinaryFunction operation)                             double initial_value) const
1600  {  {
1601    //    return reductionOp(m_offset,operation,initial_value);
1602    // use the default offset  }
1603    binaryOp(m_offset,right,operation);  
1604    template <class BinaryFunction>
1605    inline
1606    double
1607    DataArrayView::reductionOp(ValueType::size_type offset,
1608                               BinaryFunction operation,
1609                               double initial_value) const
1610    {
1611      EsysAssert((!isEmpty()&&checkOffset(offset)),
1612                   "Error - Couldn't perform reductionOp due to insufficient storage.");
1613      double current_value=initial_value;
1614      for (ValueType::size_type i=0;i<noValues();i++) {
1615        current_value=operation(current_value,(*m_data)[offset+i]);
1616      }
1617      return current_value;
1618    }
1619    
1620    inline
1621    DataArrayView::ValueType::size_type
1622    DataArrayView::relIndex() const
1623    {
1624      EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1625      return 0;
1626    }
1627    
1628    inline
1629    DataArrayView::ValueType::size_type
1630    DataArrayView::index() const
1631    {
1632      EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1633      return (m_offset);
1634    }
1635    
1636    inline
1637    DataArrayView::ValueType::size_type
1638    DataArrayView::relIndex(ValueType::size_type i) const
1639    {
1640      EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1641      EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1642      return i;
1643  }  }
1644    
1645  inline  inline
1646  DataArrayView::ValueType::size_type  DataArrayView::ValueType::size_type
1647  DataArrayView::index(ValueType::size_type i) const  DataArrayView::index(ValueType::size_type i) const
1648  {  {
1649      //EsysAssert((i >= 0) && (i < noValues(m_shape)), "Invalid index.");    EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1650      EsysAssert((i < noValues(m_shape)), "Invalid index.");    EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1651      return (i+m_offset);    return (m_offset+i);
1652  }  }
1653    
1654  inline  inline
# Line 697  DataArrayView::ValueType::size_type Line 1656  DataArrayView::ValueType::size_type
1656  DataArrayView::relIndex(ValueType::size_type i,  DataArrayView::relIndex(ValueType::size_type i,
1657                          ValueType::size_type j) const                          ValueType::size_type j) const
1658  {  {
1659      EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1660    ValueType::size_type temp=i+j*m_shape[0];    ValueType::size_type temp=i+j*m_shape[0];
1661    //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
   EsysAssert((temp < noValues(m_shape)), "Invalid index.");  
   //  
   // no offset  
1662    return temp;    return temp;
1663  }  }
1664    
# Line 710  DataArrayView::ValueType::size_type Line 1667  DataArrayView::ValueType::size_type
1667  DataArrayView::index(ValueType::size_type i,  DataArrayView::index(ValueType::size_type i,
1668               ValueType::size_type j) const               ValueType::size_type j) const
1669  {  {
1670      EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1671    ValueType::size_type temp=i+j*m_shape[0];    ValueType::size_type temp=i+j*m_shape[0];
1672    //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1673    EsysAssert((temp < noValues(m_shape)), "Invalid index.");    return (m_offset+temp);
   return (temp+m_offset);  
1674  }  }
1675    
1676  inline  inline
# Line 722  DataArrayView::relIndex(ValueType::size_ Line 1679  DataArrayView::relIndex(ValueType::size_
1679              ValueType::size_type j,              ValueType::size_type j,
1680              ValueType::size_type k) const              ValueType::size_type k) const
1681  {  {
1682      ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];    EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1683      //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1684      EsysAssert((temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1685      //    return temp;
     // no offset  
     return temp;  
1686  }  }
1687    
1688  inline  inline
# Line 736  DataArrayView::index(ValueType::size_typ Line 1691  DataArrayView::index(ValueType::size_typ
1691               ValueType::size_type j,               ValueType::size_type j,
1692               ValueType::size_type k) const               ValueType::size_type k) const
1693  {  {
1694      ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];    EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1695      //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1696      EsysAssert((temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1697      return (temp+m_offset);    return (m_offset+temp);
1698  }  }
1699    
1700  inline  inline
# Line 749  DataArrayView::relIndex(ValueType::size_ Line 1704  DataArrayView::relIndex(ValueType::size_
1704                          ValueType::size_type k,                          ValueType::size_type k,
1705                          ValueType::size_type m) const                          ValueType::size_type m) const
1706  {  {
1707      EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1708    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0]+m*m_shape[2]*m_shape[1]*m_shape[0];    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0]+m*m_shape[2]*m_shape[1]*m_shape[0];
1709    //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
   EsysAssert((temp < noValues(m_shape)), "Invalid index.");  
   //  
   // no offset  
1710    return temp;    return temp;
1711  }  }
1712    
# Line 764  DataArrayView::index(ValueType::size_typ Line 1717  DataArrayView::index(ValueType::size_typ
1717               ValueType::size_type k,               ValueType::size_type k,
1718               ValueType::size_type m) const               ValueType::size_type m) const
1719  {  {
1720      EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1721    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0]+m*m_shape[2]*m_shape[1]*m_shape[0];    ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0]+m*m_shape[2]*m_shape[1]*m_shape[0];
1722    //EsysAssert((temp >= 0 || temp < noValues(m_shape)), "Invalid index.");    EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1723    EsysAssert((temp < noValues(m_shape)), "Invalid index.");    return (m_offset+temp);
   return (temp+m_offset);  
1724  }  }
1725    
1726  inline  inline
1727  DataArrayView::ValueType::reference  DataArrayView::ValueType::reference
1728  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()()
                           ValueType::size_type j,  
                           ValueType::size_type k,  
                           ValueType::size_type m)  
1729  {  {
1730      EsysAssert((getRank()==4),    return (*m_data)[index()];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j,k,m)];  
1731  }  }
1732    
1733  inline  inline
1734  DataArrayView::ValueType::const_reference  DataArrayView::ValueType::const_reference
1735  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()() const
                           ValueType::size_type j,  
                           ValueType::size_type k,  
                           ValueType::size_type m) const  
1736  {  {
1737      EsysAssert((getRank()==4),    return (*m_data)[index()];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j,k,m)];  
1738  }  }
1739    
1740  inline  inline
1741  DataArrayView::ValueType::reference  DataArrayView::ValueType::reference
1742  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()(ValueType::size_type i)
                           ValueType::size_type j,  
                           ValueType::size_type k)  
1743  {  {
1744      EsysAssert((getRank()==3),    return (*m_data)[index(i)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j,k)];  
1745  }  }
1746    
1747  inline  inline
1748  DataArrayView::ValueType::const_reference  DataArrayView::ValueType::const_reference
1749  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()(ValueType::size_type i) const
                           ValueType::size_type j,  
                           ValueType::size_type k) const  
1750  {  {
1751      EsysAssert((getRank()==3),    return (*m_data)[index(i)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j,k)];  
1752  }  }
1753    
1754  inline  inline
# Line 821  DataArrayView::ValueType::reference Line 1756  DataArrayView::ValueType::reference
1756  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()(ValueType::size_type i,
1757                            ValueType::size_type j)                            ValueType::size_type j)
1758  {  {
1759      EsysAssert((getRank()==2),    return (*m_data)[index(i,j)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j)];  
1760  }  }
1761    
1762  inline  inline
# Line 831  DataArrayView::ValueType::const_referenc Line 1764  DataArrayView::ValueType::const_referenc
1764  DataArrayView::operator()(ValueType::size_type i,  DataArrayView::operator()(ValueType::size_type i,
1765                            ValueType::size_type j) const                            ValueType::size_type j) const
1766  {  {
1767      EsysAssert((getRank()==2),    return (*m_data)[index(i,j)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i,j)];  
1768  }  }
1769    
1770  inline  inline
1771  DataArrayView::ValueType::reference  DataArrayView::ValueType::reference
1772  DataArrayView::operator()(ValueType::size_type i)  DataArrayView::operator()(ValueType::size_type i,
1773                              ValueType::size_type j,
1774                              ValueType::size_type k)
1775  {  {
1776      EsysAssert((getRank()==1),    return (*m_data)[index(i,j,k)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i)];  
1777  }  }
1778    
1779  inline  inline
1780  DataArrayView::ValueType::const_reference  DataArrayView::ValueType::const_reference
1781  DataArrayView::operator()(ValueType::size_type i) const  DataArrayView::operator()(ValueType::size_type i,
1782                              ValueType::size_type j,
1783                              ValueType::size_type k) const
1784  {  {
1785      EsysAssert((getRank()==1),    return (*m_data)[index(i,j,k)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[index(i)];  
1786  }  }
1787    
1788  inline  inline
1789  DataArrayView::ValueType::reference  DataArrayView::ValueType::reference
1790  DataArrayView::operator()()  DataArrayView::operator()(ValueType::size_type i,
1791                              ValueType::size_type j,
1792                              ValueType::size_type k,
1793                              ValueType::size_type m)
1794  {  {
1795      EsysAssert((getRank()==0),    return (*m_data)[index(i,j,k,m)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[m_offset];  
1796  }  }
1797    
1798  inline  inline
1799  DataArrayView::ValueType::const_reference  DataArrayView::ValueType::const_reference
1800  DataArrayView::operator()() const  DataArrayView::operator()(ValueType::size_type i,
1801                              ValueType::size_type j,
1802                              ValueType::size_type k,
1803                              ValueType::size_type m) const
1804  {  {
1805      EsysAssert((getRank()==0),    return (*m_data)[index(i,j,k,m)];
            "Incorrect number of indices for the rank of this object.");  
     return (*m_data)[m_offset];  
1806  }  }
1807    
 bool operator==(const DataArrayView& left, const DataArrayView& right);  
 bool operator!=(const DataArrayView& left, const DataArrayView& right);  
   
 /* returns the python slice object key as a pair of ints where the first  
     member is the start and the second member is the end. the presence of a possible  
     step attribute with value different from one will truow an exception */  
 std::pair<int,int>  
 getSliceRange(const int s,  
               const boost::python::object&  key);  
   
1808  } // end of namespace  } // end of namespace
1809    
1810  #endif  #endif

Legend:
Removed from v.100  
changed lines
  Added in v.1388

  ViewVC Help
Powered by ViewVC 1.1.26