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

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

  ViewVC Help
Powered by ViewVC 1.1.26