/[escript]/trunk-mpi-branch/escript/src/DataArrayView.h
ViewVC logotype

Contents of /trunk-mpi-branch/escript/src/DataArrayView.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1306 - (show annotations)
Tue Sep 18 05:51:09 2007 UTC (11 years, 11 months ago) by ksteube
File MIME type: text/plain
File size: 54218 byte(s)
New Copyright in each .c .h .cpp and .py file

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26