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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 800 - (show annotations)
Tue Aug 8 11:23:18 2006 UTC (13 years, 4 months ago) by gross
File MIME type: text/plain
File size: 52849 byte(s)
new function _swap. Python wrapper + testing is still missing.


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 axis_offset and axis_offset+1
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 */
1192 ESCRIPT_DLL_API
1193 static
1194 inline
1195 void
1196 swap(DataArrayView& in,
1197 ValueType::size_type inOffset,
1198 DataArrayView& ev,
1199 ValueType::size_type evOffset,
1200 int axis_offset)
1201 {
1202 if (in.getRank() == 4) {
1203 int s0=ev.getShape()[0];
1204 int s1=ev.getShape()[1];
1205 int s2=ev.getShape()[2];
1206 int s3=ev.getShape()[3];
1207 int i0, i1, i2, i3;
1208 if (axis_offset==0) {
1209 for (i0=0; i0<s0; i0++) {
1210 for (i1=0; i1<s1; i1++) {
1211 for (i2=0; i2<s2; i2++) {
1212 for (i3=0; i3<s3; i3++) {
1213 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i0,i2,i3)];
1214 }
1215 }
1216 }
1217 }
1218 }
1219 else if (axis_offset==1) {
1220 for (i0=0; i0<s0; i0++) {
1221 for (i1=0; i1<s1; i1++) {
1222 for (i2=0; i2<s2; i2++) {
1223 for (i3=0; i3<s3; i3++) {
1224 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i2,i1,i3)];
1225 }
1226 }
1227 }
1228 }
1229 }
1230 else {
1231 for (i0=0; i0<s0; i0++) {
1232 for (i1=0; i1<s1; i1++) {
1233 for (i2=0; i2<s2; i2++) {
1234 for (i3=0; i3<s3; i3++) {
1235 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i3,i2)];
1236 }
1237 }
1238 }
1239 }
1240 }
1241 }
1242 else if (in.getRank() == 3) {
1243 int s0=ev.getShape()[0];
1244 int s1=ev.getShape()[1];
1245 int s2=ev.getShape()[2];
1246 int i0, i1, i2;
1247 if (axis_offset==0) {
1248 for (i0=0; i0<s0; i0++) {
1249 for (i1=0; i1<s1; i1++) {
1250 for (i2=0; i2<s2; i2++) {
1251 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i0,i2)];
1252 }
1253 }
1254 }
1255 }
1256 else {
1257 for (i0=0; i0<s0; i0++) {
1258 for (i1=0; i1<s1; i1++) {
1259 for (i2=0; i2<s2; i2++) {
1260 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i2,i1)];
1261 }
1262 }
1263 }
1264 }
1265 }
1266 else if (in.getRank() == 2) {
1267 int s0=ev.getShape()[0];
1268 int s1=ev.getShape()[1];
1269 int i0, i1;
1270 for (i0=0; i0<s0; i0++) {
1271 for (i1=0; i1<s1; i1++) {
1272 (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1273 }
1274 }
1275 }
1276 else if (in.getRank() == 1) {
1277 int s0=ev.getShape()[0];
1278 int i0;
1279 for (i0=0; i0<s0; i0++) {
1280 (*(ev.m_data))[evOffset+ev.index(i0)] = (*(in.m_data))[inOffset+in.index(i0)];
1281 }
1282 }
1283 else if (in.getRank() == 0) {
1284 (*(ev.m_data))[evOffset+ev.index()] = (*(in.m_data))[inOffset+in.index()];
1285 }
1286 else {
1287 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
1288 }
1289 }
1290
1291 /**
1292 \brief
1293 solves a local eigenvalue problem
1294
1295 \param in - Input - matrix
1296 \param inOffset - Input - offset into in
1297 \param ev - Output - The eigenvalues
1298 \param inOffset - Input - offset into ev
1299 */
1300 ESCRIPT_DLL_API
1301 static
1302 inline
1303 void
1304 eigenvalues(DataArrayView& in,
1305 ValueType::size_type inOffset,
1306 DataArrayView& ev,
1307 ValueType::size_type evOffset)
1308 {
1309 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1310 double ev0,ev1,ev2;
1311 int s=in.getShape()[0];
1312 if (s==1) {
1313 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1314 eigenvalues1(in00,&ev0);
1315 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1316
1317 } else if (s==2) {
1318 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1319 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1320 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1321 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1322 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
1323 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1324 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1325
1326 } else if (s==3) {
1327 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1328 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1329 in20=(*(in.m_data))[inOffset+in.index(2,0)];
1330 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1331 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1332 in21=(*(in.m_data))[inOffset+in.index(2,1)];
1333 in02=(*(in.m_data))[inOffset+in.index(0,2)];
1334 in12=(*(in.m_data))[inOffset+in.index(1,2)];
1335 in22=(*(in.m_data))[inOffset+in.index(2,2)];
1336 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1337 &ev0,&ev1,&ev2);
1338 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1339 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1340 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1341
1342 }
1343 }
1344
1345 /**
1346 \brief
1347 solves a local eigenvalue problem
1348
1349 \param in - Input - matrix
1350 \param inOffset - Input - offset into in
1351 \param ev - Output - The eigenvalues
1352 \param evOffset - Input - offset into ev
1353 \param V - Output - The eigenvectors
1354 \param VOffset - Input - offset into V
1355 \param tol - Input - eigenvalues with relative difference tol are treated as equal
1356 */
1357 ESCRIPT_DLL_API
1358 static
1359 inline
1360 void
1361 eigenvalues_and_eigenvectors(DataArrayView& in,
1362 ValueType::size_type inOffset,
1363 DataArrayView& ev,
1364 ValueType::size_type evOffset,
1365 DataArrayView& V,
1366 ValueType::size_type VOffset,
1367 const double tol=1.e-13)
1368 {
1369 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1370 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
1371 double ev0,ev1,ev2;
1372 int s=in.getShape()[0];
1373 if (s==1) {
1374 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1375 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
1376 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1377 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1378 } else if (s==2) {
1379 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1380 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1381 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1382 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1383 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
1384 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
1385 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1386 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1387 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1388 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1389 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1390 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1391 } else if (s==3) {
1392 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1393 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1394 in20=(*(in.m_data))[inOffset+in.index(2,0)];
1395 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1396 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1397 in21=(*(in.m_data))[inOffset+in.index(2,1)];
1398 in02=(*(in.m_data))[inOffset+in.index(0,2)];
1399 in12=(*(in.m_data))[inOffset+in.index(1,2)];
1400 in22=(*(in.m_data))[inOffset+in.index(2,2)];
1401 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1402 &ev0,&ev1,&ev2,
1403 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
1404 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1405 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1406 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1407 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1408 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1409 (*(V.m_data))[inOffset+V.index(2,0)]=V20;
1410 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1411 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1412 (*(V.m_data))[inOffset+V.index(2,1)]=V21;
1413 (*(V.m_data))[inOffset+V.index(0,2)]=V02;
1414 (*(V.m_data))[inOffset+V.index(1,2)]=V12;
1415 (*(V.m_data))[inOffset+V.index(2,2)]=V22;
1416
1417 }
1418 }
1419 protected:
1420
1421 private:
1422
1423 //
1424 // The maximum rank allowed for the shape of any view.
1425 static const int m_maxRank=4;
1426
1427 //
1428 // The data values for the view.
1429 // NOTE: This points to data external to the view.
1430 // This is just a pointer to an array of ValueType.
1431 ValueType* m_data;
1432
1433 //
1434 // The offset into the data array used by different views.
1435 // This is simply an integer specifying a position in the data array
1436 // pointed to by m_data.
1437 ValueType::size_type m_offset;
1438
1439 //
1440 // The shape of the data.
1441 // This is simply an STL vector specifying the lengths of each dimension
1442 // of the shape as ints.
1443 ShapeType m_shape;
1444
1445 //
1446 // The number of values needed for the array.
1447 // This can be derived from m_shape by multiplying the size of each dimension, but
1448 // is stored here for convenience.
1449 int m_noValues;
1450
1451 };
1452
1453 ESCRIPT_DLL_API bool operator==(const DataArrayView& left, const DataArrayView& right);
1454 ESCRIPT_DLL_API bool operator!=(const DataArrayView& left, const DataArrayView& right);
1455
1456 /**
1457 \brief
1458 Modify region to copy from in order to
1459 deal with the case where one range in the region contains identical indexes,
1460 eg: <<1,1><0,3><0,3>>
1461 This situation implies we want to copy from an object with rank greater than that of this
1462 object. eg: we want to copy the values from a two dimensional slice out of a three
1463 dimensional object into a two dimensional object.
1464 We do this by taking a slice from the other object where one dimension of
1465 the slice region is of size 1. So in the above example, we modify the above
1466 region like so: <<1,2><0,3><0,3>> and take this slice.
1467 */
1468 DataArrayView::RegionLoopRangeType
1469 getSliceRegionLoopRange(const DataArrayView::RegionType& region);
1470
1471 /**
1472 \brief
1473 Calculate the slice range from the given python key object
1474 Used by DataArrayView::getSliceRegion.
1475 Returns the python slice object key as a pair of ints where the first
1476 member is the start and the second member is the end. the presence of a possible
1477 step attribute with value different from one will throw an exception
1478
1479 /param key - Input - key object specifying slice range.
1480 */
1481 std::pair<int,int>
1482 getSliceRange(const boost::python::object& key,
1483 const int shape);
1484
1485 /**
1486 Inline function definitions.
1487 */
1488
1489 template <class UnaryFunction>
1490 inline
1491 void
1492 DataArrayView::unaryOp(UnaryFunction operation)
1493 {
1494 unaryOp(m_offset,operation);
1495 }
1496
1497 template <class UnaryFunction>
1498 inline
1499 void
1500 DataArrayView::unaryOp(ValueType::size_type offset,
1501 UnaryFunction operation)
1502 {
1503 EsysAssert((!isEmpty()&&checkOffset(offset)),
1504 "Error - Couldn't perform unaryOp due to insufficient storage.");
1505 for (ValueType::size_type i=0;i<noValues();i++) {
1506 (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1507 }
1508 }
1509
1510 template <class BinaryFunction>
1511 inline
1512 void
1513 DataArrayView::binaryOp(const DataArrayView& right,
1514 BinaryFunction operation)
1515 {
1516 binaryOp(m_offset,right,right.getOffset(),operation);
1517 }
1518
1519 template <class BinaryFunction>
1520 inline
1521 void
1522 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1523 const DataArrayView& right,
1524 ValueType::size_type rightOffset,
1525 BinaryFunction operation)
1526 {
1527 EsysAssert(getShape()==right.getShape(),
1528 "Error - Couldn't perform binaryOp due to shape mismatch,");
1529 EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1530 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1531 EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1532 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1533 for (ValueType::size_type i=0;i<noValues();i++) {
1534 (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1535 }
1536 }
1537
1538 template <class BinaryFunction>
1539 inline
1540 void
1541 DataArrayView::binaryOp(double right,
1542 BinaryFunction operation)
1543 {
1544 binaryOp(m_offset,right,operation);
1545 }
1546
1547 template <class BinaryFunction>
1548 inline
1549 void
1550 DataArrayView::binaryOp(ValueType::size_type offset,
1551 double right,
1552 BinaryFunction operation)
1553 {
1554 EsysAssert((!isEmpty()&&checkOffset(offset)),
1555 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1556 for (ValueType::size_type i=0;i<noValues();i++) {
1557 (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1558 }
1559 }
1560
1561 template <class BinaryFunction>
1562 inline
1563 double
1564 DataArrayView::reductionOp(BinaryFunction operation,
1565 double initial_value) const
1566 {
1567 return reductionOp(m_offset,operation,initial_value);
1568 }
1569
1570 template <class BinaryFunction>
1571 inline
1572 double
1573 DataArrayView::reductionOp(ValueType::size_type offset,
1574 BinaryFunction operation,
1575 double initial_value) const
1576 {
1577 EsysAssert((!isEmpty()&&checkOffset(offset)),
1578 "Error - Couldn't perform reductionOp due to insufficient storage.");
1579 double current_value=initial_value;
1580 for (ValueType::size_type i=0;i<noValues();i++) {
1581 current_value=operation(current_value,(*m_data)[offset+i]);
1582 }
1583 return current_value;
1584 }
1585
1586 inline
1587 DataArrayView::ValueType::size_type
1588 DataArrayView::relIndex() const
1589 {
1590 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1591 return 0;
1592 }
1593
1594 inline
1595 DataArrayView::ValueType::size_type
1596 DataArrayView::index() const
1597 {
1598 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1599 return (m_offset);
1600 }
1601
1602 inline
1603 DataArrayView::ValueType::size_type
1604 DataArrayView::relIndex(ValueType::size_type i) const
1605 {
1606 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1607 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1608 return i;
1609 }
1610
1611 inline
1612 DataArrayView::ValueType::size_type
1613 DataArrayView::index(ValueType::size_type i) const
1614 {
1615 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1616 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1617 return (m_offset+i);
1618 }
1619
1620 inline
1621 DataArrayView::ValueType::size_type
1622 DataArrayView::relIndex(ValueType::size_type i,
1623 ValueType::size_type j) const
1624 {
1625 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1626 ValueType::size_type temp=i+j*m_shape[0];
1627 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1628 return temp;
1629 }
1630
1631 inline
1632 DataArrayView::ValueType::size_type
1633 DataArrayView::index(ValueType::size_type i,
1634 ValueType::size_type j) const
1635 {
1636 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1637 ValueType::size_type temp=i+j*m_shape[0];
1638 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1639 return (m_offset+temp);
1640 }
1641
1642 inline
1643 DataArrayView::ValueType::size_type
1644 DataArrayView::relIndex(ValueType::size_type i,
1645 ValueType::size_type j,
1646 ValueType::size_type k) const
1647 {
1648 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1649 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1650 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1651 return temp;
1652 }
1653
1654 inline
1655 DataArrayView::ValueType::size_type
1656 DataArrayView::index(ValueType::size_type i,
1657 ValueType::size_type j,
1658 ValueType::size_type k) const
1659 {
1660 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1661 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1662 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1663 return (m_offset+temp);
1664 }
1665
1666 inline
1667 DataArrayView::ValueType::size_type
1668 DataArrayView::relIndex(ValueType::size_type i,
1669 ValueType::size_type j,
1670 ValueType::size_type k,
1671 ValueType::size_type m) const
1672 {
1673 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1674 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];
1675 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1676 return temp;
1677 }
1678
1679 inline
1680 DataArrayView::ValueType::size_type
1681 DataArrayView::index(ValueType::size_type i,
1682 ValueType::size_type j,
1683 ValueType::size_type k,
1684 ValueType::size_type m) const
1685 {
1686 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1687 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];
1688 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1689 return (m_offset+temp);
1690 }
1691
1692 inline
1693 DataArrayView::ValueType::reference
1694 DataArrayView::operator()()
1695 {
1696 return (*m_data)[index()];
1697 }
1698
1699 inline
1700 DataArrayView::ValueType::const_reference
1701 DataArrayView::operator()() const
1702 {
1703 return (*m_data)[index()];
1704 }
1705
1706 inline
1707 DataArrayView::ValueType::reference
1708 DataArrayView::operator()(ValueType::size_type i)
1709 {
1710 return (*m_data)[index(i)];
1711 }
1712
1713 inline
1714 DataArrayView::ValueType::const_reference
1715 DataArrayView::operator()(ValueType::size_type i) const
1716 {
1717 return (*m_data)[index(i)];
1718 }
1719
1720 inline
1721 DataArrayView::ValueType::reference
1722 DataArrayView::operator()(ValueType::size_type i,
1723 ValueType::size_type j)
1724 {
1725 return (*m_data)[index(i,j)];
1726 }
1727
1728 inline
1729 DataArrayView::ValueType::const_reference
1730 DataArrayView::operator()(ValueType::size_type i,
1731 ValueType::size_type j) const
1732 {
1733 return (*m_data)[index(i,j)];
1734 }
1735
1736 inline
1737 DataArrayView::ValueType::reference
1738 DataArrayView::operator()(ValueType::size_type i,
1739 ValueType::size_type j,
1740 ValueType::size_type k)
1741 {
1742 return (*m_data)[index(i,j,k)];
1743 }
1744
1745 inline
1746 DataArrayView::ValueType::const_reference
1747 DataArrayView::operator()(ValueType::size_type i,
1748 ValueType::size_type j,
1749 ValueType::size_type k) const
1750 {
1751 return (*m_data)[index(i,j,k)];
1752 }
1753
1754 inline
1755 DataArrayView::ValueType::reference
1756 DataArrayView::operator()(ValueType::size_type i,
1757 ValueType::size_type j,
1758 ValueType::size_type k,
1759 ValueType::size_type m)
1760 {
1761 return (*m_data)[index(i,j,k,m)];
1762 }
1763
1764 inline
1765 DataArrayView::ValueType::const_reference
1766 DataArrayView::operator()(ValueType::size_type i,
1767 ValueType::size_type j,
1768 ValueType::size_type k,
1769 ValueType::size_type m) const
1770 {
1771 return (*m_data)[index(i,j,k,m)];
1772 }
1773
1774 } // end of namespace
1775
1776 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26