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

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

  ViewVC Help
Powered by ViewVC 1.1.26