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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 775 - (hide annotations)
Mon Jul 10 04:00:08 2006 UTC (15 years ago) by ksteube
File MIME type: text/plain
File size: 49720 byte(s)
Modified the following python methods in escript/py_src/util.py to
call faster C++ methods:
	escript_trace
	escript_transpose
	escript_symmetric
	escript_nonsymmetric

1 jgs 102 // $Id$
2 jgs 108
3 jgs 82 /*
4 elspeth 615 ************************************************************
5     * Copyright 2006 by ACcESS MNRF *
6     * *
7     * http://www.access.edu.au *
8     * Primary Business: Queensland, Australia *
9     * Licensed under the Open Software License version 3.0 *
10     * http://www.opensource.org/licenses/osl-3.0.php *
11     * *
12     ************************************************************
13 jgs 82 */
14 jgs 474
15 jgs 82 #if !defined escript_DataArrayView_20040323_H
16     #define escript_DataArrayView_20040323_H
17 woo409 757 #include "system_dep.h"
18 jgs 82
19 robwdcock 682 #include "esysUtils/EsysAssert.h"
20 jgs 82
21 ksteube 775 #include "DataException.h"
22 jgs 474 #include "DataVector.h"
23 gross 576 #include "LocalOps.h"
24 jgs 474
25 jgs 82 #include <boost/python/numeric.hpp>
26     #include <boost/python/object.hpp>
27    
28 jgs 108 #include <vector>
29 jgs 102
30 jgs 82 namespace escript {
31    
32     /**
33     \brief
34 jgs 108 DataArrayView provides a view of external data, configured according
35     to the given shape information and offset.
36 jgs 82
37     Description:
38     DataArrayView provides a view of data allocated externally. The
39     external data should provide sufficient values so that the dimensions
40 jgs 108 specified for the view shape will be satisfied. The viewer can update
41 jgs 82 values within the external data but cannot resize the external data.
42 jgs 108
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 jgs 117 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 jgs 82 */
50    
51     class DataArrayView {
52    
53 woo409 757 ESCRIPT_DLL_API friend bool operator==(const DataArrayView& left, const DataArrayView& right);
54     ESCRIPT_DLL_API friend bool operator!=(const DataArrayView& left, const DataArrayView& right);
55 jgs 82
56     public:
57    
58 jgs 117 //
59     // Some basic types which define the data values and view shapes.
60     typedef DataVector ValueType;
61     //typedef std::vector<double> ValueType;
62 jgs 108 typedef std::vector<int> ShapeType;
63 jgs 82 typedef std::vector<std::pair<int, int> > RegionType;
64 jgs 108 typedef std::vector<std::pair<int, int> > RegionLoopRangeType;
65 jgs 82
66     /**
67     \brief
68     Default constructor for DataArrayView.
69    
70     Description:
71 jgs 108 Default constructor for DataArrayView.
72     Creates an empty view with no associated data array.
73    
74     This is essentially useless but here for completeness.
75 jgs 82 */
76 woo409 757 ESCRIPT_DLL_API
77 jgs 82 DataArrayView();
78    
79     /**
80     \brief
81     Constructor for DataArrayView.
82    
83     Description:
84     Constructor for DataArrayView.
85    
86 jgs 108 \param data - Input -
87     Array holding data to be viewed. This must be large enough
88     to supply sufficient values for the specified shape starting at
89     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 jgs 82 */
95 woo409 757 ESCRIPT_DLL_API
96 jgs 82 DataArrayView(ValueType& data,
97     const ShapeType& viewShape,
98     int offset=0);
99    
100     /**
101     \brief
102     Copy constructor for DataArrayView.
103    
104     Description:
105     Copy constructor for DataArrayView.
106    
107 jgs 108 \param other - Input -
108     DataArrayView to copy.
109 jgs 82
110 jgs 108 NOTE: The copy references the same data array.
111 jgs 82 */
112 woo409 757 ESCRIPT_DLL_API
113 jgs 82 DataArrayView(const DataArrayView& other);
114    
115     /**
116     \brief
117 jgs 108 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.
119 jgs 82 */
120 woo409 757 ESCRIPT_DLL_API
121 jgs 82 void
122     copy(const boost::python::numeric::array& value);
123    
124     /**
125     \brief
126 jgs 108 Copy data from another DataArrayView into the data array viewed by this
127     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 jgs 82 */
131 woo409 757 ESCRIPT_DLL_API
132 jgs 82 void
133     copy(const DataArrayView& other);
134    
135     /**
136     \brief
137 jgs 108 Copy data from another DataArrayView into this view starting at the
138     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 jgs 82 */
142 woo409 757 ESCRIPT_DLL_API
143 jgs 82 void
144     copy(ValueType::size_type offset,
145     const DataArrayView& other);
146    
147     /**
148     \brief
149 jgs 108 Copy data from another DataArrayView into this view starting at the
150     given offsets in each view.
151     The shapes of both views must be compatible - will throw if not.
152     NB: views may or may not reference same underlying data array!
153 jgs 82
154 jgs 108 \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 jgs 82 */
161 woo409 757 ESCRIPT_DLL_API
162 jgs 82 void
163     copy(ValueType::size_type thisOffset,
164     const DataArrayView& other,
165     ValueType::size_type otherOffset);
166    
167     /**
168     \brief
169 jgs 108 Copy the given single value over the entire view starting at the given
170     offset in the view's data array.
171 jgs 82
172 jgs 108 \param offset - Input -
173     Offset into this view's data array to copy to.
174     \param value - Input -
175     Value to copy.
176 jgs 82 */
177 woo409 757 ESCRIPT_DLL_API
178 jgs 82 void
179     copy(ValueType::size_type offset,
180     ValueType::value_type value);
181    
182     /**
183 jgs 108 \brief
184     Check if view is empty. ie: does not point to any actual data.
185     */
186 woo409 757 ESCRIPT_DLL_API
187 jgs 108 bool
188     isEmpty() const;
189    
190     /**
191 jgs 82 \brief
192 jgs 108 Return this view's offset into the viewed data array.
193 jgs 82 */
194 woo409 757 ESCRIPT_DLL_API
195 jgs 82 ValueType::size_type
196     getOffset() const;
197    
198     /**
199 jgs 108 \brief
200     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 jgs 82 */
205 woo409 757 ESCRIPT_DLL_API
206 jgs 108 void
207     setOffset(ValueType::size_type offset);
208 jgs 82
209     /**
210     \brief
211 jgs 108 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 jgs 82 */
215 woo409 757 ESCRIPT_DLL_API
216 jgs 108 void
217     incrOffset();
218 jgs 82
219     /**
220     \brief
221 jgs 108 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 jgs 82 */
224 woo409 757 ESCRIPT_DLL_API
225 jgs 108 bool
226     checkOffset() const;
227 jgs 117
228 woo409 757 ESCRIPT_DLL_API
229 jgs 108 bool
230     checkOffset(ValueType::size_type offset) const;
231 jgs 82
232     /**
233 jgs 108 \brief
234     Return the rank of the shape of this view.
235     */
236 woo409 757 ESCRIPT_DLL_API
237 jgs 108 int
238     getRank() const;
239    
240     /**
241 jgs 82 \brief
242 jgs 108 Return the number of values for the shape of this view.
243 jgs 82 */
244 woo409 757 ESCRIPT_DLL_API
245 jgs 82 int
246     noValues() const;
247    
248     /**
249     \brief
250 jgs 108 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 jgs 82 */
253 woo409 757 ESCRIPT_DLL_API
254 jgs 117 static
255     int
256 jgs 108 noValues(const ShapeType& shape);
257 jgs 117
258 woo409 757 ESCRIPT_DLL_API
259 jgs 117 static
260     int
261 jgs 108 noValues(const RegionLoopRangeType& region);
262 jgs 82
263     /**
264     \brief
265 jgs 108 Return a reference to the underlying data array.
266 jgs 82 */
267 woo409 757 ESCRIPT_DLL_API
268 jgs 82 ValueType&
269     getData() const;
270    
271     /**
272     \brief
273 jgs 108 Return a reference to the data value with the given
274     index in this view. This takes into account the offset.
275 jgs 82 */
276 woo409 757 ESCRIPT_DLL_API
277 jgs 82 ValueType::reference
278     getData(ValueType::size_type i) const;
279    
280     /**
281     \brief
282 jgs 108 Return the shape of this view.
283 jgs 82 */
284 woo409 757 ESCRIPT_DLL_API
285 jgs 117 const
286     ShapeType&
287 jgs 108 getShape() const;
288 jgs 82
289     /**
290     \brief
291 jgs 108 Return true if the given shape is the same as this view's shape.
292 jgs 82 */
293 woo409 757 ESCRIPT_DLL_API
294 jgs 108 bool
295 jgs 117 checkShape(const ShapeType& other) const;
296 jgs 82
297     /**
298     \brief
299 jgs 108 Create a shape error message. Normally used when there is a shape
300 jgs 117 mismatch between this shape and the other shape.
301 jgs 82
302 jgs 108 \param messagePrefix - Input -
303     First part of the error message.
304     \param other - Input -
305     The other shape.
306 jgs 82 */
307 woo409 757 ESCRIPT_DLL_API
308 jgs 108 std::string
309     createShapeErrorMessage(const std::string& messagePrefix,
310 jgs 117 const ShapeType& other) const;
311 jgs 82
312     /**
313     \brief
314 jgs 108 Return the viewed data as a formatted string.
315     Not recommended for very large objects!
316 jgs 82
317 jgs 108 \param suffix - Input -
318     Suffix appended to index display.
319 jgs 82 */
320 woo409 757 ESCRIPT_DLL_API
321 jgs 108 std::string
322     toString(const std::string& suffix=std::string("")) const;
323 jgs 82
324     /**
325     \brief
326 jgs 108 Return the given shape as a string.
327     This is purely a utility method and has no bearing on this view.
328 jgs 82
329 jgs 108 \param shape - Input.
330 jgs 82 */
331 woo409 757 ESCRIPT_DLL_API
332 jgs 117 static
333     std::string
334     shapeToString(const ShapeType& shape);
335 jgs 82
336 jgs 108 /**
337     \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 jgs 82 */
342 woo409 757 ESCRIPT_DLL_API
343 jgs 108 ValueType::size_type
344     relIndex() const;
345 jgs 82
346     /**
347 jgs 108 \brief
348     Return the 1 dimensional index into the data array of the element at
349     position i in the view, *ignoring the offset*.
350     Assumes a rank 1 view.
351 jgs 82 */
352 woo409 757 ESCRIPT_DLL_API
353 jgs 108 ValueType::size_type
354     relIndex(ValueType::size_type i) const;
355 jgs 82
356     /**
357 jgs 108 \brief
358     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 jgs 82 */
362 woo409 757 ESCRIPT_DLL_API
363 jgs 108 ValueType::size_type
364     relIndex(ValueType::size_type i,
365     ValueType::size_type j) const;
366 jgs 82
367     /**
368 jgs 108 \brief
369     Return the 1 dimensional index into the data array of the element at
370     position i,j,k in the view, *ignoring the offset*.
371     Assumes a rank 3 view.
372 jgs 82 */
373 woo409 757 ESCRIPT_DLL_API
374 jgs 108 ValueType::size_type
375     relIndex(ValueType::size_type i,
376     ValueType::size_type j,
377     ValueType::size_type k) const;
378 jgs 82
379     /**
380     \brief
381 jgs 108 Return the 1 dimensional index into the data array of the element at
382     position i,j,k,m in the view, *ignoring the offset*.
383     Assumes a rank 4 view.
384 jgs 82 */
385 woo409 757 ESCRIPT_DLL_API
386 jgs 82 ValueType::size_type
387     relIndex(ValueType::size_type i,
388     ValueType::size_type j,
389     ValueType::size_type k,
390     ValueType::size_type m) const;
391    
392     /**
393     \brief
394 jgs 108 Return the 1 dimensional index into the data array of the only element
395     in the view.
396     Assumes a rank 0 view.
397 jgs 82 */
398 woo409 757 ESCRIPT_DLL_API
399 jgs 82 ValueType::size_type
400 jgs 108 index() const;
401 jgs 82
402     /**
403     \brief
404 jgs 108 Return the 1 dimensional index into the data array of the element at
405     position i in the view.
406     Assumes a rank 1 view.
407 jgs 82 */
408 woo409 757 ESCRIPT_DLL_API
409 jgs 82 ValueType::size_type
410 jgs 108 index(ValueType::size_type i) const;
411 jgs 82
412     /**
413     \brief
414 jgs 108 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 jgs 82 */
418 woo409 757 ESCRIPT_DLL_API
419 jgs 82 ValueType::size_type
420     index(ValueType::size_type i,
421 jgs 108 ValueType::size_type j) const;
422 jgs 82
423     /**
424     \brief
425 jgs 108 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 jgs 82 */
429 woo409 757 ESCRIPT_DLL_API
430 jgs 82 ValueType::size_type
431     index(ValueType::size_type i,
432     ValueType::size_type j,
433     ValueType::size_type k) const;
434    
435     /**
436     \brief
437 jgs 108 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 jgs 82 */
441 woo409 757 ESCRIPT_DLL_API
442 jgs 82 ValueType::size_type
443     index(ValueType::size_type i,
444 jgs 108 ValueType::size_type j,
445     ValueType::size_type k,
446     ValueType::size_type m) const;
447 jgs 82
448     /**
449     \brief
450 jgs 108 Return a reference for the only element in the view.
451     Assumes a rank 0 view.
452 jgs 82 */
453 woo409 757 ESCRIPT_DLL_API
454 jgs 108 ValueType::reference
455     operator()();
456 jgs 82
457 woo409 757 ESCRIPT_DLL_API
458 jgs 108 ValueType::const_reference
459     operator()() const;
460    
461 jgs 82 /**
462     \brief
463 jgs 108 Return a reference to the element at position i in the view.
464     Assumes a rank 1 view.
465 jgs 82 */
466 woo409 757 ESCRIPT_DLL_API
467 jgs 82 ValueType::reference
468     operator()(ValueType::size_type i);
469    
470 woo409 757 ESCRIPT_DLL_API
471 jgs 82 ValueType::const_reference
472     operator()(ValueType::size_type i) const;
473    
474     /**
475     \brief
476 jgs 108 Return a reference to the element at position i,j in the view.
477     Assumes a rank 2 view.
478 jgs 82 */
479 woo409 757 ESCRIPT_DLL_API
480 jgs 82 ValueType::reference
481     operator()(ValueType::size_type i,
482     ValueType::size_type j);
483    
484 woo409 757 ESCRIPT_DLL_API
485 jgs 82 ValueType::const_reference
486     operator()(ValueType::size_type i,
487     ValueType::size_type j) const;
488    
489     /**
490     \brief
491 jgs 108 Return a reference to the element at position i,j,k in the view.
492     Assumes a rank 3 view.
493 jgs 82 */
494 woo409 757 ESCRIPT_DLL_API
495 jgs 82 ValueType::reference
496     operator()(ValueType::size_type i,
497     ValueType::size_type j,
498     ValueType::size_type k);
499    
500 woo409 757 ESCRIPT_DLL_API
501 jgs 82 ValueType::const_reference
502     operator()(ValueType::size_type i,
503     ValueType::size_type j,
504     ValueType::size_type k) const;
505    
506     /**
507     \brief
508 jgs 108 Return a reference to the element at position i,j,k,m in the view.
509     Assumes a rank 4 view.
510 jgs 82 */
511 woo409 757 ESCRIPT_DLL_API
512 jgs 82 ValueType::reference
513     operator()(ValueType::size_type i,
514     ValueType::size_type j,
515     ValueType::size_type k,
516     ValueType::size_type m);
517    
518 woo409 757 ESCRIPT_DLL_API
519 jgs 82 ValueType::const_reference
520     operator()(ValueType::size_type i,
521     ValueType::size_type j,
522     ValueType::size_type k,
523     ValueType::size_type m) const;
524    
525     /**
526 jgs 108 \brief
527     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 jgs 82 */
533 woo409 757 ESCRIPT_DLL_API
534 jgs 117 static
535     ShapeType
536 jgs 108 getResultSliceShape(const RegionType& region);
537 jgs 82
538 jgs 108 /**
539     \brief
540     Determine the region specified by the given python slice object.
541 jgs 82
542 jgs 108 \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 woo409 757 ESCRIPT_DLL_API
593 jgs 117 RegionType
594 jgs 108 getSliceRegion(const boost::python::object& key) const;
595    
596 jgs 82 /**
597     \brief
598 jgs 108 Copy a data slice specified by the given region from the given view
599     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 jgs 82 */
606 woo409 757 ESCRIPT_DLL_API
607 jgs 82 void
608 jgs 108 copySlice(const DataArrayView& other,
609     const RegionLoopRangeType& region);
610 jgs 82
611     /**
612     \brief
613 jgs 108 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 jgs 82 */
625 woo409 757 ESCRIPT_DLL_API
626 jgs 82 void
627 jgs 108 copySlice(ValueType::size_type thisOffset,
628     const DataArrayView& other,
629     ValueType::size_type otherOffset,
630     const RegionLoopRangeType& region);
631 jgs 82
632     /**
633     \brief
634 jgs 108 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 jgs 102
637 jgs 108 \param other - Input -
638     View to copy data from.
639     \param region - Input -
640     Region in this view to copy data to.
641 jgs 102 */
642 woo409 757 ESCRIPT_DLL_API
643 jgs 108 void
644     copySliceFrom(const DataArrayView& other,
645     const RegionLoopRangeType& region);
646 jgs 102
647     /**
648     \brief
649 jgs 108 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 jgs 102
652 jgs 108 \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 jgs 102 */
661 woo409 757 ESCRIPT_DLL_API
662 jgs 108 void
663     copySliceFrom(ValueType::size_type thisOffset,
664     const DataArrayView& other,
665     ValueType::size_type otherOffset,
666     const RegionLoopRangeType& region);
667 jgs 102
668     /**
669     \brief
670 jgs 108 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 jgs 82 */
679     template <class UnaryFunction>
680 woo409 757 ESCRIPT_DLL_API
681 jgs 108 void
682     unaryOp(UnaryFunction operation);
683 jgs 82
684     /**
685     \brief
686 jgs 108 Perform the unary operation on the data point specified by the given
687     offset. Applies the specified operation to each value in the data
688     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 jgs 82 */
697     template <class UnaryFunction>
698 woo409 757 ESCRIPT_DLL_API
699 jgs 108 void
700     unaryOp(ValueType::size_type offset,
701     UnaryFunction operation);
702 jgs 82
703     /**
704     \brief
705 jgs 108 Perform the binary operation on the data points specified by the default
706     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 jgs 82 */
717     template <class BinaryFunction>
718 woo409 757 ESCRIPT_DLL_API
719 jgs 82 void
720 jgs 108 binaryOp(const DataArrayView& right,
721 jgs 82 BinaryFunction operation);
722    
723     /**
724     \brief
725 jgs 108 Perform the binary operation on the data points specified by the given
726     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 jgs 82 */
741     template <class BinaryFunction>
742 woo409 757 ESCRIPT_DLL_API
743 jgs 82 void
744 jgs 108 binaryOp(ValueType::size_type leftOffset,
745     const DataArrayView& right,
746     ValueType::size_type rightOffset,
747 jgs 82 BinaryFunction operation);
748    
749     /**
750     \brief
751 jgs 108 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 jgs 82 */
763     template <class BinaryFunction>
764 woo409 757 ESCRIPT_DLL_API
765 jgs 82 void
766 jgs 108 binaryOp(double right,
767 jgs 82 BinaryFunction operation);
768    
769     /**
770     \brief
771 jgs 108 Perform the binary operation on the data point specified by the given
772     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 jgs 82 */
785     template <class BinaryFunction>
786 woo409 757 ESCRIPT_DLL_API
787 jgs 82 void
788 jgs 108 binaryOp(ValueType::size_type offset,
789     double right,
790 jgs 82 BinaryFunction operation);
791    
792     /**
793     \brief
794 jgs 108 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     Called by escript::algorithm.
800    
801     \param operation - Input -
802     Operation to apply. Must be a pointer to a function.
803     */
804 jgs 147 template <class BinaryFunction>
805 woo409 757 ESCRIPT_DLL_API
806 jgs 108 double
807 jgs 147 reductionOp(BinaryFunction operation,
808     double initial_value) const;
809 jgs 108
810     /**
811     \brief
812 jgs 113 Perform the given data point reduction operation on the data point
813 jgs 108 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     Called by escript::algorithm.
818    
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 jgs 147 template <class BinaryFunction>
825 woo409 757 ESCRIPT_DLL_API
826 jgs 108 double
827     reductionOp(ValueType::size_type offset,
828 jgs 147 BinaryFunction operation,
829     double initial_value) const;
830 jgs 108
831 gross 584 /**
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.
841     \param right - Input - The right hand side.
842     \param result - Output - The result of the operation.
843     */
844 woo409 757 ESCRIPT_DLL_API
845 gross 584 static
846     void
847     matMult(const DataArrayView& left,
848     const DataArrayView& right,
849     DataArrayView& result);
850    
851 jgs 108 /**
852     \brief
853 gross 584 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 woo409 757 ESCRIPT_DLL_API
858 gross 584 static
859     ShapeType
860     determineResultShape(const DataArrayView& left,
861     const DataArrayView& right);
862    
863     /**
864     \brief
865 ksteube 775 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     static
873     inline
874     void
875     symmetric(DataArrayView& in,
876     ValueType::size_type inOffset,
877     DataArrayView& ev,
878     ValueType::size_type evOffset)
879     {
880     if (in.getRank() == 2) {
881     int i0, i1;
882     int s0=in.getShape()[0];
883     int s1=in.getShape()[1];
884     for (i0=0; i0<s0; i0++) {
885     for (i1=0; i1<s1; i1++) {
886     (*(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;
887     }
888     }
889     }
890     if (in.getRank() == 4) {
891     int i0, i1, i2, i3;
892     int s0=in.getShape()[0];
893     int s1=in.getShape()[1];
894     int s2=in.getShape()[2];
895     int s3=in.getShape()[3];
896     for (i0=0; i0<s0; i0++) {
897     for (i1=0; i1<s1; i1++) {
898     for (i2=0; i2<s2; i2++) {
899     for (i3=0; i3<s3; i3++) {
900     (*(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;
901     }
902     }
903     }
904     }
905     }
906     }
907    
908     /**
909     \brief
910     computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
911    
912     \param in - Input - matrix
913     \param inOffset - Input - offset into in
914     \param ev - Output - The nonsymmetric matrix
915     \param inOffset - Input - offset into ev
916     */
917     static
918     inline
919     void
920     nonsymmetric(DataArrayView& in,
921     ValueType::size_type inOffset,
922     DataArrayView& ev,
923     ValueType::size_type evOffset)
924     {
925     if (in.getRank() == 2) {
926     int i0, i1;
927     int s0=in.getShape()[0];
928     int s1=in.getShape()[1];
929     for (i0=0; i0<s0; i0++) {
930     for (i1=0; i1<s1; i1++) {
931     (*(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;
932     }
933     }
934     }
935     if (in.getRank() == 4) {
936     int i0, i1, i2, i3;
937     int s0=in.getShape()[0];
938     int s1=in.getShape()[1];
939     int s2=in.getShape()[2];
940     int s3=in.getShape()[3];
941     for (i0=0; i0<s0; i0++) {
942     for (i1=0; i1<s1; i1++) {
943     for (i2=0; i2<s2; i2++) {
944     for (i3=0; i3<s3; i3++) {
945     (*(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;
946     }
947     }
948     }
949     }
950     }
951     }
952    
953     /**
954     \brief
955     computes the trace of a matrix
956    
957     \param in - Input - matrix
958     \param inOffset - Input - offset into in
959     \param ev - Output - The nonsymmetric matrix
960     \param inOffset - Input - offset into ev
961     */
962     static
963     inline
964     void
965     matrixtrace(DataArrayView& in,
966     ValueType::size_type inOffset,
967     DataArrayView& ev,
968     ValueType::size_type evOffset,
969     int axis_offset)
970     {
971     if (in.getRank() == 2) {
972     int s0=in.getShape()[0]; // Python wrapper limits to square matrix
973     int i;
974     for (i=0; i<s0; i++) {
975     (*(ev.m_data))[evOffset+ev.index()] += (*(in.m_data))[inOffset+in.index(i,i)];
976     }
977     }
978     else if (in.getRank() == 3) {
979     if (axis_offset==0) {
980     int s0=in.getShape()[0];
981     int s2=in.getShape()[2];
982     int i0, i2;
983     for (i0=0; i0<s0; i0++) {
984     for (i2=0; i2<s2; i2++) {
985     (*(ev.m_data))[evOffset+ev.index(i2)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2)];
986     }
987     }
988     }
989     else if (axis_offset==1) {
990     int s0=in.getShape()[0];
991     int s1=in.getShape()[1];
992     int i0, i1;
993     for (i0=0; i0<s0; i0++) {
994     for (i1=0; i1<s1; i1++) {
995     (*(ev.m_data))[evOffset+ev.index(i0)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1)];
996     }
997     }
998     }
999     }
1000     else if (in.getRank() == 4) {
1001     if (axis_offset==0) {
1002     int s0=in.getShape()[0];
1003     int s2=in.getShape()[2];
1004     int s3=in.getShape()[3];
1005     int i0, i2, i3;
1006     for (i0=0; i0<s0; i0++) {
1007     for (i2=0; i2<s2; i2++) {
1008     for (i3=0; i3<s3; i3++) {
1009     (*(ev.m_data))[evOffset+ev.index(i2,i3)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2,i3)];
1010     }
1011     }
1012     }
1013     }
1014     else if (axis_offset==1) {
1015     int s0=in.getShape()[0];
1016     int s1=in.getShape()[1];
1017     int s3=in.getShape()[3];
1018     int i0, i1, i3;
1019     for (i0=0; i0<s0; i0++) {
1020     for (i1=0; i1<s1; i1++) {
1021     for (i3=0; i3<s3; i3++) {
1022     (*(ev.m_data))[evOffset+ev.index(i0,i3)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1,i3)];
1023     }
1024     }
1025     }
1026     }
1027     else if (axis_offset==2) {
1028     int s0=in.getShape()[0];
1029     int s1=in.getShape()[1];
1030     int s2=in.getShape()[2];
1031     int i0, i1, i2;
1032     for (i0=0; i0<s0; i0++) {
1033     for (i1=0; i1<s1; i1++) {
1034     for (i2=0; i2<s2; i2++) {
1035     (*(ev.m_data))[evOffset+ev.index(i0,i1)] += (*(in.m_data))[inOffset+in.index(i0,i1,i2,i2)];
1036     }
1037     }
1038     }
1039     }
1040     }
1041     }
1042    
1043     /**
1044     \brief
1045     Transpose each data point of this Data object around the given axis.
1046    
1047     \param in - Input - matrix
1048     \param inOffset - Input - offset into in
1049     \param ev - Output - The nonsymmetric matrix
1050     \param inOffset - Input - offset into ev
1051     */
1052     static
1053     inline
1054     void
1055     transpose(DataArrayView& in,
1056     ValueType::size_type inOffset,
1057     DataArrayView& ev,
1058     ValueType::size_type evOffset,
1059     int axis_offset)
1060     {
1061     if (in.getRank() == 4) {
1062     int s0=ev.getShape()[0];
1063     int s1=ev.getShape()[1];
1064     int s2=ev.getShape()[2];
1065     int s3=ev.getShape()[3];
1066     int i0, i1, i2, i3;
1067     if (axis_offset==1) {
1068     for (i0=0; i0<s0; i0++) {
1069     for (i1=0; i1<s1; i1++) {
1070     for (i2=0; i2<s2; i2++) {
1071     for (i3=0; i3<s3; i3++) {
1072     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i3,i0,i1,i2)];
1073     }
1074     }
1075     }
1076     }
1077     }
1078     else if (axis_offset==2) {
1079     for (i0=0; i0<s0; i0++) {
1080     for (i1=0; i1<s1; i1++) {
1081     for (i2=0; i2<s2; i2++) {
1082     for (i3=0; i3<s3; i3++) {
1083     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)];
1084     }
1085     }
1086     }
1087     }
1088     }
1089     else if (axis_offset==3) {
1090     for (i0=0; i0<s0; i0++) {
1091     for (i1=0; i1<s1; i1++) {
1092     for (i2=0; i2<s2; i2++) {
1093     for (i3=0; i3<s3; i3++) {
1094     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i2,i3,i0)];
1095     }
1096     }
1097     }
1098     }
1099     }
1100     else {
1101     for (i0=0; i0<s0; i0++) {
1102     for (i1=0; i1<s1; i1++) {
1103     for (i2=0; i2<s2; i2++) {
1104     for (i3=0; i3<s3; i3++) {
1105     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)];
1106     }
1107     }
1108     }
1109     }
1110     }
1111     }
1112     else if (in.getRank() == 3) {
1113     int s0=ev.getShape()[0];
1114     int s1=ev.getShape()[1];
1115     int s2=ev.getShape()[2];
1116     int i0, i1, i2;
1117     if (axis_offset==1) {
1118     for (i0=0; i0<s0; i0++) {
1119     for (i1=0; i1<s1; i1++) {
1120     for (i2=0; i2<s2; i2++) {
1121     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i2,i0,i1)];
1122     }
1123     }
1124     }
1125     }
1126     else if (axis_offset==2) {
1127     for (i0=0; i0<s0; i0++) {
1128     for (i1=0; i1<s1; i1++) {
1129     for (i2=0; i2<s2; i2++) {
1130     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i2,i0)];
1131     }
1132     }
1133     }
1134     }
1135     else {
1136     // Copy the matrix unchanged
1137     for (i0=0; i0<s0; i0++) {
1138     for (i1=0; i1<s1; i1++) {
1139     for (i2=0; i2<s2; i2++) {
1140     (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2)];
1141     }
1142     }
1143     }
1144     }
1145     }
1146     else if (in.getRank() == 2) {
1147     int s0=ev.getShape()[0];
1148     int s1=ev.getShape()[1];
1149     int i0, i1;
1150     if (axis_offset==1) {
1151     for (i0=0; i0<s0; i0++) {
1152     for (i1=0; i1<s1; i1++) {
1153     (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1154     }
1155     }
1156     }
1157     else {
1158     for (i0=0; i0<s0; i0++) {
1159     for (i1=0; i1<s1; i1++) {
1160     (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i0,i1)];
1161     }
1162     }
1163     }
1164     }
1165     else if (in.getRank() == 1) {
1166     int s0=ev.getShape()[0];
1167     int i0;
1168     for (i0=0; i0<s0; i0++) {
1169     (*(ev.m_data))[evOffset+ev.index(i0)] = (*(in.m_data))[inOffset+in.index(i0)];
1170     }
1171     }
1172     else if (in.getRank() == 0) {
1173     (*(ev.m_data))[evOffset+ev.index()] = (*(in.m_data))[inOffset+in.index()];
1174     }
1175     else {
1176     throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
1177     }
1178     }
1179    
1180     /**
1181     \brief
1182 gross 576 solves a local eigenvalue problem
1183    
1184 gross 580 \param in - Input - matrix
1185     \param inOffset - Input - offset into in
1186 gross 576 \param ev - Output - The eigenvalues
1187 gross 580 \param inOffset - Input - offset into ev
1188 gross 576 */
1189 woo409 757 ESCRIPT_DLL_API
1190 gross 580 static
1191 gross 576 inline
1192     void
1193 gross 584 eigenvalues(DataArrayView& in,
1194     ValueType::size_type inOffset,
1195     DataArrayView& ev,
1196     ValueType::size_type evOffset)
1197     {
1198 gross 580 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1199     double ev0,ev1,ev2;
1200     int s=in.getShape()[0];
1201     if (s==1) {
1202     in00=(*(in.m_data))[inOffset+in.index(0,0)];
1203     eigenvalues1(in00,&ev0);
1204     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1205 gross 576
1206 gross 580 } else if (s==2) {
1207     in00=(*(in.m_data))[inOffset+in.index(0,0)];
1208     in10=(*(in.m_data))[inOffset+in.index(1,0)];
1209     in01=(*(in.m_data))[inOffset+in.index(0,1)];
1210     in11=(*(in.m_data))[inOffset+in.index(1,1)];
1211     eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
1212     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1213     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1214 gross 576
1215 gross 580 } else if (s==3) {
1216     in00=(*(in.m_data))[inOffset+in.index(0,0)];
1217     in10=(*(in.m_data))[inOffset+in.index(1,0)];
1218     in20=(*(in.m_data))[inOffset+in.index(2,0)];
1219     in01=(*(in.m_data))[inOffset+in.index(0,1)];
1220     in11=(*(in.m_data))[inOffset+in.index(1,1)];
1221     in21=(*(in.m_data))[inOffset+in.index(2,1)];
1222     in02=(*(in.m_data))[inOffset+in.index(0,2)];
1223     in12=(*(in.m_data))[inOffset+in.index(1,2)];
1224     in22=(*(in.m_data))[inOffset+in.index(2,2)];
1225     eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1226     &ev0,&ev1,&ev2);
1227     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1228     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1229     (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1230    
1231     }
1232 gross 584 }
1233 gross 580
1234 gross 576 /**
1235     \brief
1236 gross 580 solves a local eigenvalue problem
1237    
1238     \param in - Input - matrix
1239     \param inOffset - Input - offset into in
1240     \param ev - Output - The eigenvalues
1241     \param evOffset - Input - offset into ev
1242     \param V - Output - The eigenvectors
1243     \param VOffset - Input - offset into V
1244     \param tol - Input - eigenvalues with relative difference tol are treated as equal
1245     */
1246 woo409 757 ESCRIPT_DLL_API
1247 gross 580 static
1248     inline
1249     void
1250 gross 584 eigenvalues_and_eigenvectors(DataArrayView& in,
1251     ValueType::size_type inOffset,
1252     DataArrayView& ev,
1253     ValueType::size_type evOffset,
1254     DataArrayView& V,
1255     ValueType::size_type VOffset,
1256     const double tol=1.e-13)
1257 gross 580 {
1258 gross 583 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1259     double V00,V10,V20,V01,V11,V21,V02,V12,V22;
1260 gross 580 double ev0,ev1,ev2;
1261     int s=in.getShape()[0];
1262     if (s==1) {
1263 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1264     eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
1265     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1266     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1267 gross 580 } else if (s==2) {
1268 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1269     in10=(*(in.m_data))[inOffset+in.index(1,0)];
1270     in01=(*(in.m_data))[inOffset+in.index(0,1)];
1271     in11=(*(in.m_data))[inOffset+in.index(1,1)];
1272     eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
1273     &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
1274     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1275     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1276     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1277     (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1278     (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1279     (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1280 gross 580 } else if (s==3) {
1281 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1282     in10=(*(in.m_data))[inOffset+in.index(1,0)];
1283     in20=(*(in.m_data))[inOffset+in.index(2,0)];
1284     in01=(*(in.m_data))[inOffset+in.index(0,1)];
1285     in11=(*(in.m_data))[inOffset+in.index(1,1)];
1286     in21=(*(in.m_data))[inOffset+in.index(2,1)];
1287     in02=(*(in.m_data))[inOffset+in.index(0,2)];
1288     in12=(*(in.m_data))[inOffset+in.index(1,2)];
1289     in22=(*(in.m_data))[inOffset+in.index(2,2)];
1290     eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1291     &ev0,&ev1,&ev2,
1292     &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
1293     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1294     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1295     (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1296     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1297     (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1298     (*(V.m_data))[inOffset+V.index(2,0)]=V20;
1299     (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1300     (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1301     (*(V.m_data))[inOffset+V.index(2,1)]=V21;
1302     (*(V.m_data))[inOffset+V.index(0,2)]=V02;
1303     (*(V.m_data))[inOffset+V.index(1,2)]=V12;
1304     (*(V.m_data))[inOffset+V.index(2,2)]=V22;
1305 gross 580
1306     }
1307     }
1308 jgs 82 protected:
1309    
1310     private:
1311    
1312     //
1313 jgs 108 // The maximum rank allowed for the shape of any view.
1314 jgs 82 static const int m_maxRank=4;
1315    
1316     //
1317 jgs 108 // The data values for the view.
1318     // NOTE: This points to data external to the view.
1319 jgs 117 // This is just a pointer to an array of ValueType.
1320 jgs 82 ValueType* m_data;
1321    
1322     //
1323     // The offset into the data array used by different views.
1324 jgs 108 // This is simply an integer specifying a position in the data array
1325     // pointed to by m_data.
1326 jgs 82 ValueType::size_type m_offset;
1327    
1328     //
1329     // The shape of the data.
1330 jgs 108 // This is simply an STL vector specifying the lengths of each dimension
1331     // of the shape as ints.
1332 jgs 82 ShapeType m_shape;
1333    
1334     //
1335     // The number of values needed for the array.
1336 jgs 108 // This can be derived from m_shape by multiplying the size of each dimension, but
1337     // is stored here for convenience.
1338 jgs 82 int m_noValues;
1339    
1340     };
1341    
1342 woo409 757 ESCRIPT_DLL_API bool operator==(const DataArrayView& left, const DataArrayView& right);
1343     ESCRIPT_DLL_API bool operator!=(const DataArrayView& left, const DataArrayView& right);
1344 jgs 108
1345 jgs 102 /**
1346     \brief
1347 jgs 108 Modify region to copy from in order to
1348     deal with the case where one range in the region contains identical indexes,
1349     eg: <<1,1><0,3><0,3>>
1350     This situation implies we want to copy from an object with rank greater than that of this
1351     object. eg: we want to copy the values from a two dimensional slice out of a three
1352     dimensional object into a two dimensional object.
1353     We do this by taking a slice from the other object where one dimension of
1354     the slice region is of size 1. So in the above example, we modify the above
1355     region like so: <<1,2><0,3><0,3>> and take this slice.
1356     */
1357     DataArrayView::RegionLoopRangeType
1358     getSliceRegionLoopRange(const DataArrayView::RegionType& region);
1359    
1360     /**
1361     \brief
1362 jgs 102 Calculate the slice range from the given python key object
1363     Used by DataArrayView::getSliceRegion.
1364     Returns the python slice object key as a pair of ints where the first
1365     member is the start and the second member is the end. the presence of a possible
1366     step attribute with value different from one will throw an exception
1367    
1368     /param key - Input - key object specifying slice range.
1369     */
1370     std::pair<int,int>
1371     getSliceRange(const boost::python::object& key,
1372     const int shape);
1373    
1374 jgs 108 /**
1375     Inline function definitions.
1376     */
1377    
1378     template <class UnaryFunction>
1379 jgs 82 inline
1380 jgs 108 void
1381     DataArrayView::unaryOp(UnaryFunction operation)
1382 jgs 82 {
1383 jgs 108 unaryOp(m_offset,operation);
1384 jgs 82 }
1385    
1386 jgs 108 template <class UnaryFunction>
1387 jgs 82 inline
1388 jgs 108 void
1389     DataArrayView::unaryOp(ValueType::size_type offset,
1390     UnaryFunction operation)
1391 jgs 82 {
1392 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1393     "Error - Couldn't perform unaryOp due to insufficient storage.");
1394     for (ValueType::size_type i=0;i<noValues();i++) {
1395     (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1396     }
1397 jgs 82 }
1398    
1399 jgs 108 template <class BinaryFunction>
1400 jgs 82 inline
1401 jgs 108 void
1402     DataArrayView::binaryOp(const DataArrayView& right,
1403     BinaryFunction operation)
1404 jgs 82 {
1405 jgs 108 binaryOp(m_offset,right,right.getOffset(),operation);
1406 jgs 82 }
1407    
1408 jgs 108 template <class BinaryFunction>
1409 jgs 82 inline
1410     void
1411 jgs 108 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1412     const DataArrayView& right,
1413     ValueType::size_type rightOffset,
1414     BinaryFunction operation)
1415 jgs 82 {
1416 jgs 108 EsysAssert(getShape()==right.getShape(),
1417     "Error - Couldn't perform binaryOp due to shape mismatch,");
1418     EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1419     "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1420     EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1421     "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1422     for (ValueType::size_type i=0;i<noValues();i++) {
1423     (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1424 jgs 82 }
1425     }
1426    
1427 jgs 108 template <class BinaryFunction>
1428 jgs 82 inline
1429     void
1430 jgs 108 DataArrayView::binaryOp(double right,
1431     BinaryFunction operation)
1432 jgs 82 {
1433 jgs 108 binaryOp(m_offset,right,operation);
1434 jgs 82 }
1435    
1436 jgs 108 template <class BinaryFunction>
1437 jgs 82 inline
1438 jgs 108 void
1439     DataArrayView::binaryOp(ValueType::size_type offset,
1440     double right,
1441     BinaryFunction operation)
1442 jgs 102 {
1443 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1444 jgs 113 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1445 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1446     (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1447 jgs 102 }
1448     }
1449    
1450 jgs 147 template <class BinaryFunction>
1451 jgs 102 inline
1452     double
1453 jgs 147 DataArrayView::reductionOp(BinaryFunction operation,
1454     double initial_value) const
1455 jgs 82 {
1456 jgs 147 return reductionOp(m_offset,operation,initial_value);
1457 jgs 82 }
1458    
1459 jgs 147 template <class BinaryFunction>
1460 jgs 82 inline
1461 jgs 108 double
1462     DataArrayView::reductionOp(ValueType::size_type offset,
1463 jgs 147 BinaryFunction operation,
1464     double initial_value) const
1465 jgs 82 {
1466 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1467     "Error - Couldn't perform reductionOp due to insufficient storage.");
1468 jgs 147 double current_value=initial_value;
1469 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1470 jgs 147 current_value=operation(current_value,(*m_data)[offset+i]);
1471 jgs 82 }
1472 jgs 147 return current_value;
1473 jgs 82 }
1474    
1475     inline
1476 jgs 108 DataArrayView::ValueType::size_type
1477     DataArrayView::relIndex() const
1478 jgs 82 {
1479 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1480     return 0;
1481 jgs 82 }
1482    
1483     inline
1484 jgs 108 DataArrayView::ValueType::size_type
1485     DataArrayView::index() const
1486 jgs 82 {
1487 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1488     return (m_offset);
1489 jgs 82 }
1490    
1491     inline
1492 jgs 108 DataArrayView::ValueType::size_type
1493     DataArrayView::relIndex(ValueType::size_type i) const
1494 jgs 82 {
1495 jgs 108 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1496     EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1497     return i;
1498 jgs 82 }
1499    
1500     inline
1501     DataArrayView::ValueType::size_type
1502     DataArrayView::index(ValueType::size_type i) const
1503     {
1504 jgs 108 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1505     EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1506     return (m_offset+i);
1507 jgs 82 }
1508    
1509     inline
1510     DataArrayView::ValueType::size_type
1511     DataArrayView::relIndex(ValueType::size_type i,
1512     ValueType::size_type j) const
1513     {
1514 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1515 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1516 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1517 jgs 82 return temp;
1518     }
1519    
1520     inline
1521     DataArrayView::ValueType::size_type
1522     DataArrayView::index(ValueType::size_type i,
1523     ValueType::size_type j) const
1524     {
1525 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1526 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1527 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1528     return (m_offset+temp);
1529 jgs 82 }
1530    
1531     inline
1532     DataArrayView::ValueType::size_type
1533     DataArrayView::relIndex(ValueType::size_type i,
1534     ValueType::size_type j,
1535     ValueType::size_type k) const
1536     {
1537 jgs 108 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1538     ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1539     EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1540     return temp;
1541 jgs 82 }
1542    
1543     inline
1544     DataArrayView::ValueType::size_type
1545     DataArrayView::index(ValueType::size_type i,
1546     ValueType::size_type j,
1547     ValueType::size_type k) const
1548     {
1549 jgs 108 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1550     ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1551     EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1552     return (m_offset+temp);
1553 jgs 82 }
1554    
1555     inline
1556     DataArrayView::ValueType::size_type
1557     DataArrayView::relIndex(ValueType::size_type i,
1558     ValueType::size_type j,
1559     ValueType::size_type k,
1560     ValueType::size_type m) const
1561     {
1562 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1563 jgs 82 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];
1564 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1565 jgs 82 return temp;
1566     }
1567    
1568     inline
1569     DataArrayView::ValueType::size_type
1570     DataArrayView::index(ValueType::size_type i,
1571     ValueType::size_type j,
1572     ValueType::size_type k,
1573     ValueType::size_type m) const
1574     {
1575 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1576 jgs 82 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];
1577 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1578     return (m_offset+temp);
1579 jgs 82 }
1580    
1581     inline
1582     DataArrayView::ValueType::reference
1583 jgs 108 DataArrayView::operator()()
1584 jgs 82 {
1585 jgs 108 return (*m_data)[index()];
1586 jgs 82 }
1587    
1588     inline
1589     DataArrayView::ValueType::const_reference
1590 jgs 108 DataArrayView::operator()() const
1591 jgs 82 {
1592 jgs 108 return (*m_data)[index()];
1593 jgs 82 }
1594    
1595     inline
1596     DataArrayView::ValueType::reference
1597 jgs 108 DataArrayView::operator()(ValueType::size_type i)
1598 jgs 82 {
1599 jgs 108 return (*m_data)[index(i)];
1600 jgs 82 }
1601    
1602     inline
1603     DataArrayView::ValueType::const_reference
1604 jgs 108 DataArrayView::operator()(ValueType::size_type i) const
1605 jgs 82 {
1606 jgs 108 return (*m_data)[index(i)];
1607 jgs 82 }
1608    
1609     inline
1610     DataArrayView::ValueType::reference
1611     DataArrayView::operator()(ValueType::size_type i,
1612     ValueType::size_type j)
1613     {
1614 jgs 108 return (*m_data)[index(i,j)];
1615 jgs 82 }
1616    
1617     inline
1618     DataArrayView::ValueType::const_reference
1619     DataArrayView::operator()(ValueType::size_type i,
1620     ValueType::size_type j) const
1621     {
1622 jgs 108 return (*m_data)[index(i,j)];
1623 jgs 82 }
1624    
1625     inline
1626     DataArrayView::ValueType::reference
1627 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1628     ValueType::size_type j,
1629     ValueType::size_type k)
1630 jgs 82 {
1631 jgs 108 return (*m_data)[index(i,j,k)];
1632 jgs 82 }
1633    
1634     inline
1635     DataArrayView::ValueType::const_reference
1636 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1637     ValueType::size_type j,
1638     ValueType::size_type k) const
1639 jgs 82 {
1640 jgs 108 return (*m_data)[index(i,j,k)];
1641 jgs 82 }
1642    
1643     inline
1644     DataArrayView::ValueType::reference
1645 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1646     ValueType::size_type j,
1647     ValueType::size_type k,
1648     ValueType::size_type m)
1649 jgs 82 {
1650 jgs 108 return (*m_data)[index(i,j,k,m)];
1651 jgs 82 }
1652    
1653     inline
1654     DataArrayView::ValueType::const_reference
1655 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1656     ValueType::size_type j,
1657     ValueType::size_type k,
1658     ValueType::size_type m) const
1659 jgs 82 {
1660 jgs 108 return (*m_data)[index(i,j,k,m)];
1661 jgs 82 }
1662    
1663     } // end of namespace
1664    
1665     #endif

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26