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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1030 - (show annotations)
Wed Mar 14 05:14:44 2007 UTC (12 years, 6 months ago) by phornby
File MIME type: text/plain
File size: 54353 byte(s)
Assemble_CopyElementData.c - resolve with my local edits.
Mesh_saveVTK.c - remove debug print.
DataArrayView.h - remove dllim/export declarations on template functions.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26