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

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

  ViewVC Help
Powered by ViewVC 1.1.26