/[escript]/branches/windows_from_1456_trunk_1522_merged_in/escript/src/DataArrayView.h
ViewVC logotype

Contents of /branches/windows_from_1456_trunk_1522_merged_in/escript/src/DataArrayView.h

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26