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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 804 - (show annotations)
Thu Aug 10 01:12:16 2006 UTC (13 years, 2 months ago) by gross
File MIME type: text/plain
File size: 54471 byte(s)
the new function swap_axes + tests added. (It replaces swap).


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26