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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 775 - (show annotations)
Mon Jul 10 04:00:08 2006 UTC (13 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 49720 byte(s)
Modified the following python methods in escript/py_src/util.py to
call faster C++ methods:
	escript_trace
	escript_transpose
	escript_symmetric
	escript_nonsymmetric

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 static
873 inline
874 void
875 symmetric(DataArrayView& in,
876 ValueType::size_type inOffset,
877 DataArrayView& ev,
878 ValueType::size_type evOffset)
879 {
880 if (in.getRank() == 2) {
881 int i0, i1;
882 int s0=in.getShape()[0];
883 int s1=in.getShape()[1];
884 for (i0=0; i0<s0; i0++) {
885 for (i1=0; i1<s1; i1++) {
886 (*(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;
887 }
888 }
889 }
890 if (in.getRank() == 4) {
891 int i0, i1, i2, i3;
892 int s0=in.getShape()[0];
893 int s1=in.getShape()[1];
894 int s2=in.getShape()[2];
895 int s3=in.getShape()[3];
896 for (i0=0; i0<s0; i0++) {
897 for (i1=0; i1<s1; i1++) {
898 for (i2=0; i2<s2; i2++) {
899 for (i3=0; i3<s3; i3++) {
900 (*(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;
901 }
902 }
903 }
904 }
905 }
906 }
907
908 /**
909 \brief
910 computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
911
912 \param in - Input - matrix
913 \param inOffset - Input - offset into in
914 \param ev - Output - The nonsymmetric matrix
915 \param inOffset - Input - offset into ev
916 */
917 static
918 inline
919 void
920 nonsymmetric(DataArrayView& in,
921 ValueType::size_type inOffset,
922 DataArrayView& ev,
923 ValueType::size_type evOffset)
924 {
925 if (in.getRank() == 2) {
926 int i0, i1;
927 int s0=in.getShape()[0];
928 int s1=in.getShape()[1];
929 for (i0=0; i0<s0; i0++) {
930 for (i1=0; i1<s1; i1++) {
931 (*(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;
932 }
933 }
934 }
935 if (in.getRank() == 4) {
936 int i0, i1, i2, i3;
937 int s0=in.getShape()[0];
938 int s1=in.getShape()[1];
939 int s2=in.getShape()[2];
940 int s3=in.getShape()[3];
941 for (i0=0; i0<s0; i0++) {
942 for (i1=0; i1<s1; i1++) {
943 for (i2=0; i2<s2; i2++) {
944 for (i3=0; i3<s3; i3++) {
945 (*(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;
946 }
947 }
948 }
949 }
950 }
951 }
952
953 /**
954 \brief
955 computes the trace of a matrix
956
957 \param in - Input - matrix
958 \param inOffset - Input - offset into in
959 \param ev - Output - The nonsymmetric matrix
960 \param inOffset - Input - offset into ev
961 */
962 static
963 inline
964 void
965 matrixtrace(DataArrayView& in,
966 ValueType::size_type inOffset,
967 DataArrayView& ev,
968 ValueType::size_type evOffset,
969 int axis_offset)
970 {
971 if (in.getRank() == 2) {
972 int s0=in.getShape()[0]; // Python wrapper limits to square matrix
973 int i;
974 for (i=0; i<s0; i++) {
975 (*(ev.m_data))[evOffset+ev.index()] += (*(in.m_data))[inOffset+in.index(i,i)];
976 }
977 }
978 else if (in.getRank() == 3) {
979 if (axis_offset==0) {
980 int s0=in.getShape()[0];
981 int s2=in.getShape()[2];
982 int i0, i2;
983 for (i0=0; i0<s0; i0++) {
984 for (i2=0; i2<s2; i2++) {
985 (*(ev.m_data))[evOffset+ev.index(i2)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2)];
986 }
987 }
988 }
989 else if (axis_offset==1) {
990 int s0=in.getShape()[0];
991 int s1=in.getShape()[1];
992 int i0, i1;
993 for (i0=0; i0<s0; i0++) {
994 for (i1=0; i1<s1; i1++) {
995 (*(ev.m_data))[evOffset+ev.index(i0)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1)];
996 }
997 }
998 }
999 }
1000 else if (in.getRank() == 4) {
1001 if (axis_offset==0) {
1002 int s0=in.getShape()[0];
1003 int s2=in.getShape()[2];
1004 int s3=in.getShape()[3];
1005 int i0, i2, i3;
1006 for (i0=0; i0<s0; i0++) {
1007 for (i2=0; i2<s2; i2++) {
1008 for (i3=0; i3<s3; i3++) {
1009 (*(ev.m_data))[evOffset+ev.index(i2,i3)] += (*(in.m_data))[inOffset+in.index(i0,i0,i2,i3)];
1010 }
1011 }
1012 }
1013 }
1014 else if (axis_offset==1) {
1015 int s0=in.getShape()[0];
1016 int s1=in.getShape()[1];
1017 int s3=in.getShape()[3];
1018 int i0, i1, i3;
1019 for (i0=0; i0<s0; i0++) {
1020 for (i1=0; i1<s1; i1++) {
1021 for (i3=0; i3<s3; i3++) {
1022 (*(ev.m_data))[evOffset+ev.index(i0,i3)] += (*(in.m_data))[inOffset+in.index(i0,i1,i1,i3)];
1023 }
1024 }
1025 }
1026 }
1027 else if (axis_offset==2) {
1028 int s0=in.getShape()[0];
1029 int s1=in.getShape()[1];
1030 int s2=in.getShape()[2];
1031 int i0, i1, i2;
1032 for (i0=0; i0<s0; i0++) {
1033 for (i1=0; i1<s1; i1++) {
1034 for (i2=0; i2<s2; i2++) {
1035 (*(ev.m_data))[evOffset+ev.index(i0,i1)] += (*(in.m_data))[inOffset+in.index(i0,i1,i2,i2)];
1036 }
1037 }
1038 }
1039 }
1040 }
1041 }
1042
1043 /**
1044 \brief
1045 Transpose each data point of this Data object around the given axis.
1046
1047 \param in - Input - matrix
1048 \param inOffset - Input - offset into in
1049 \param ev - Output - The nonsymmetric matrix
1050 \param inOffset - Input - offset into ev
1051 */
1052 static
1053 inline
1054 void
1055 transpose(DataArrayView& in,
1056 ValueType::size_type inOffset,
1057 DataArrayView& ev,
1058 ValueType::size_type evOffset,
1059 int axis_offset)
1060 {
1061 if (in.getRank() == 4) {
1062 int s0=ev.getShape()[0];
1063 int s1=ev.getShape()[1];
1064 int s2=ev.getShape()[2];
1065 int s3=ev.getShape()[3];
1066 int i0, i1, i2, i3;
1067 if (axis_offset==1) {
1068 for (i0=0; i0<s0; i0++) {
1069 for (i1=0; i1<s1; i1++) {
1070 for (i2=0; i2<s2; i2++) {
1071 for (i3=0; i3<s3; i3++) {
1072 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i3,i0,i1,i2)];
1073 }
1074 }
1075 }
1076 }
1077 }
1078 else if (axis_offset==2) {
1079 for (i0=0; i0<s0; i0++) {
1080 for (i1=0; i1<s1; i1++) {
1081 for (i2=0; i2<s2; i2++) {
1082 for (i3=0; i3<s3; i3++) {
1083 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i2,i3,i0,i1)];
1084 }
1085 }
1086 }
1087 }
1088 }
1089 else if (axis_offset==3) {
1090 for (i0=0; i0<s0; i0++) {
1091 for (i1=0; i1<s1; i1++) {
1092 for (i2=0; i2<s2; i2++) {
1093 for (i3=0; i3<s3; i3++) {
1094 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i1,i2,i3,i0)];
1095 }
1096 }
1097 }
1098 }
1099 }
1100 else {
1101 for (i0=0; i0<s0; i0++) {
1102 for (i1=0; i1<s1; i1++) {
1103 for (i2=0; i2<s2; i2++) {
1104 for (i3=0; i3<s3; i3++) {
1105 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2,i3)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2,i3)];
1106 }
1107 }
1108 }
1109 }
1110 }
1111 }
1112 else if (in.getRank() == 3) {
1113 int s0=ev.getShape()[0];
1114 int s1=ev.getShape()[1];
1115 int s2=ev.getShape()[2];
1116 int i0, i1, i2;
1117 if (axis_offset==1) {
1118 for (i0=0; i0<s0; i0++) {
1119 for (i1=0; i1<s1; i1++) {
1120 for (i2=0; i2<s2; i2++) {
1121 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i2,i0,i1)];
1122 }
1123 }
1124 }
1125 }
1126 else if (axis_offset==2) {
1127 for (i0=0; i0<s0; i0++) {
1128 for (i1=0; i1<s1; i1++) {
1129 for (i2=0; i2<s2; i2++) {
1130 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i1,i2,i0)];
1131 }
1132 }
1133 }
1134 }
1135 else {
1136 // Copy the matrix unchanged
1137 for (i0=0; i0<s0; i0++) {
1138 for (i1=0; i1<s1; i1++) {
1139 for (i2=0; i2<s2; i2++) {
1140 (*(ev.m_data))[evOffset+ev.index(i0,i1,i2)] = (*(in.m_data))[inOffset+in.index(i0,i1,i2)];
1141 }
1142 }
1143 }
1144 }
1145 }
1146 else if (in.getRank() == 2) {
1147 int s0=ev.getShape()[0];
1148 int s1=ev.getShape()[1];
1149 int i0, i1;
1150 if (axis_offset==1) {
1151 for (i0=0; i0<s0; i0++) {
1152 for (i1=0; i1<s1; i1++) {
1153 (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i1,i0)];
1154 }
1155 }
1156 }
1157 else {
1158 for (i0=0; i0<s0; i0++) {
1159 for (i1=0; i1<s1; i1++) {
1160 (*(ev.m_data))[evOffset+ev.index(i0,i1)] = (*(in.m_data))[inOffset+in.index(i0,i1)];
1161 }
1162 }
1163 }
1164 }
1165 else if (in.getRank() == 1) {
1166 int s0=ev.getShape()[0];
1167 int i0;
1168 for (i0=0; i0<s0; i0++) {
1169 (*(ev.m_data))[evOffset+ev.index(i0)] = (*(in.m_data))[inOffset+in.index(i0)];
1170 }
1171 }
1172 else if (in.getRank() == 0) {
1173 (*(ev.m_data))[evOffset+ev.index()] = (*(in.m_data))[inOffset+in.index()];
1174 }
1175 else {
1176 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
1177 }
1178 }
1179
1180 /**
1181 \brief
1182 solves a local eigenvalue problem
1183
1184 \param in - Input - matrix
1185 \param inOffset - Input - offset into in
1186 \param ev - Output - The eigenvalues
1187 \param inOffset - Input - offset into ev
1188 */
1189 ESCRIPT_DLL_API
1190 static
1191 inline
1192 void
1193 eigenvalues(DataArrayView& in,
1194 ValueType::size_type inOffset,
1195 DataArrayView& ev,
1196 ValueType::size_type evOffset)
1197 {
1198 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1199 double ev0,ev1,ev2;
1200 int s=in.getShape()[0];
1201 if (s==1) {
1202 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1203 eigenvalues1(in00,&ev0);
1204 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1205
1206 } else if (s==2) {
1207 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1208 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1209 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1210 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1211 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
1212 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1213 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1214
1215 } else if (s==3) {
1216 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1217 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1218 in20=(*(in.m_data))[inOffset+in.index(2,0)];
1219 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1220 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1221 in21=(*(in.m_data))[inOffset+in.index(2,1)];
1222 in02=(*(in.m_data))[inOffset+in.index(0,2)];
1223 in12=(*(in.m_data))[inOffset+in.index(1,2)];
1224 in22=(*(in.m_data))[inOffset+in.index(2,2)];
1225 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1226 &ev0,&ev1,&ev2);
1227 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1228 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1229 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1230
1231 }
1232 }
1233
1234 /**
1235 \brief
1236 solves a local eigenvalue problem
1237
1238 \param in - Input - matrix
1239 \param inOffset - Input - offset into in
1240 \param ev - Output - The eigenvalues
1241 \param evOffset - Input - offset into ev
1242 \param V - Output - The eigenvectors
1243 \param VOffset - Input - offset into V
1244 \param tol - Input - eigenvalues with relative difference tol are treated as equal
1245 */
1246 ESCRIPT_DLL_API
1247 static
1248 inline
1249 void
1250 eigenvalues_and_eigenvectors(DataArrayView& in,
1251 ValueType::size_type inOffset,
1252 DataArrayView& ev,
1253 ValueType::size_type evOffset,
1254 DataArrayView& V,
1255 ValueType::size_type VOffset,
1256 const double tol=1.e-13)
1257 {
1258 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
1259 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
1260 double ev0,ev1,ev2;
1261 int s=in.getShape()[0];
1262 if (s==1) {
1263 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1264 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
1265 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1266 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1267 } else if (s==2) {
1268 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1269 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1270 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1271 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1272 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
1273 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
1274 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1275 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1276 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1277 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1278 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1279 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1280 } else if (s==3) {
1281 in00=(*(in.m_data))[inOffset+in.index(0,0)];
1282 in10=(*(in.m_data))[inOffset+in.index(1,0)];
1283 in20=(*(in.m_data))[inOffset+in.index(2,0)];
1284 in01=(*(in.m_data))[inOffset+in.index(0,1)];
1285 in11=(*(in.m_data))[inOffset+in.index(1,1)];
1286 in21=(*(in.m_data))[inOffset+in.index(2,1)];
1287 in02=(*(in.m_data))[inOffset+in.index(0,2)];
1288 in12=(*(in.m_data))[inOffset+in.index(1,2)];
1289 in22=(*(in.m_data))[inOffset+in.index(2,2)];
1290 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
1291 &ev0,&ev1,&ev2,
1292 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
1293 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
1294 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
1295 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
1296 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
1297 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
1298 (*(V.m_data))[inOffset+V.index(2,0)]=V20;
1299 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
1300 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
1301 (*(V.m_data))[inOffset+V.index(2,1)]=V21;
1302 (*(V.m_data))[inOffset+V.index(0,2)]=V02;
1303 (*(V.m_data))[inOffset+V.index(1,2)]=V12;
1304 (*(V.m_data))[inOffset+V.index(2,2)]=V22;
1305
1306 }
1307 }
1308 protected:
1309
1310 private:
1311
1312 //
1313 // The maximum rank allowed for the shape of any view.
1314 static const int m_maxRank=4;
1315
1316 //
1317 // The data values for the view.
1318 // NOTE: This points to data external to the view.
1319 // This is just a pointer to an array of ValueType.
1320 ValueType* m_data;
1321
1322 //
1323 // The offset into the data array used by different views.
1324 // This is simply an integer specifying a position in the data array
1325 // pointed to by m_data.
1326 ValueType::size_type m_offset;
1327
1328 //
1329 // The shape of the data.
1330 // This is simply an STL vector specifying the lengths of each dimension
1331 // of the shape as ints.
1332 ShapeType m_shape;
1333
1334 //
1335 // The number of values needed for the array.
1336 // This can be derived from m_shape by multiplying the size of each dimension, but
1337 // is stored here for convenience.
1338 int m_noValues;
1339
1340 };
1341
1342 ESCRIPT_DLL_API bool operator==(const DataArrayView& left, const DataArrayView& right);
1343 ESCRIPT_DLL_API bool operator!=(const DataArrayView& left, const DataArrayView& right);
1344
1345 /**
1346 \brief
1347 Modify region to copy from in order to
1348 deal with the case where one range in the region contains identical indexes,
1349 eg: <<1,1><0,3><0,3>>
1350 This situation implies we want to copy from an object with rank greater than that of this
1351 object. eg: we want to copy the values from a two dimensional slice out of a three
1352 dimensional object into a two dimensional object.
1353 We do this by taking a slice from the other object where one dimension of
1354 the slice region is of size 1. So in the above example, we modify the above
1355 region like so: <<1,2><0,3><0,3>> and take this slice.
1356 */
1357 DataArrayView::RegionLoopRangeType
1358 getSliceRegionLoopRange(const DataArrayView::RegionType& region);
1359
1360 /**
1361 \brief
1362 Calculate the slice range from the given python key object
1363 Used by DataArrayView::getSliceRegion.
1364 Returns the python slice object key as a pair of ints where the first
1365 member is the start and the second member is the end. the presence of a possible
1366 step attribute with value different from one will throw an exception
1367
1368 /param key - Input - key object specifying slice range.
1369 */
1370 std::pair<int,int>
1371 getSliceRange(const boost::python::object& key,
1372 const int shape);
1373
1374 /**
1375 Inline function definitions.
1376 */
1377
1378 template <class UnaryFunction>
1379 inline
1380 void
1381 DataArrayView::unaryOp(UnaryFunction operation)
1382 {
1383 unaryOp(m_offset,operation);
1384 }
1385
1386 template <class UnaryFunction>
1387 inline
1388 void
1389 DataArrayView::unaryOp(ValueType::size_type offset,
1390 UnaryFunction operation)
1391 {
1392 EsysAssert((!isEmpty()&&checkOffset(offset)),
1393 "Error - Couldn't perform unaryOp due to insufficient storage.");
1394 for (ValueType::size_type i=0;i<noValues();i++) {
1395 (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1396 }
1397 }
1398
1399 template <class BinaryFunction>
1400 inline
1401 void
1402 DataArrayView::binaryOp(const DataArrayView& right,
1403 BinaryFunction operation)
1404 {
1405 binaryOp(m_offset,right,right.getOffset(),operation);
1406 }
1407
1408 template <class BinaryFunction>
1409 inline
1410 void
1411 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1412 const DataArrayView& right,
1413 ValueType::size_type rightOffset,
1414 BinaryFunction operation)
1415 {
1416 EsysAssert(getShape()==right.getShape(),
1417 "Error - Couldn't perform binaryOp due to shape mismatch,");
1418 EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1419 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1420 EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1421 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1422 for (ValueType::size_type i=0;i<noValues();i++) {
1423 (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1424 }
1425 }
1426
1427 template <class BinaryFunction>
1428 inline
1429 void
1430 DataArrayView::binaryOp(double right,
1431 BinaryFunction operation)
1432 {
1433 binaryOp(m_offset,right,operation);
1434 }
1435
1436 template <class BinaryFunction>
1437 inline
1438 void
1439 DataArrayView::binaryOp(ValueType::size_type offset,
1440 double right,
1441 BinaryFunction operation)
1442 {
1443 EsysAssert((!isEmpty()&&checkOffset(offset)),
1444 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1445 for (ValueType::size_type i=0;i<noValues();i++) {
1446 (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1447 }
1448 }
1449
1450 template <class BinaryFunction>
1451 inline
1452 double
1453 DataArrayView::reductionOp(BinaryFunction operation,
1454 double initial_value) const
1455 {
1456 return reductionOp(m_offset,operation,initial_value);
1457 }
1458
1459 template <class BinaryFunction>
1460 inline
1461 double
1462 DataArrayView::reductionOp(ValueType::size_type offset,
1463 BinaryFunction operation,
1464 double initial_value) const
1465 {
1466 EsysAssert((!isEmpty()&&checkOffset(offset)),
1467 "Error - Couldn't perform reductionOp due to insufficient storage.");
1468 double current_value=initial_value;
1469 for (ValueType::size_type i=0;i<noValues();i++) {
1470 current_value=operation(current_value,(*m_data)[offset+i]);
1471 }
1472 return current_value;
1473 }
1474
1475 inline
1476 DataArrayView::ValueType::size_type
1477 DataArrayView::relIndex() const
1478 {
1479 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1480 return 0;
1481 }
1482
1483 inline
1484 DataArrayView::ValueType::size_type
1485 DataArrayView::index() const
1486 {
1487 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1488 return (m_offset);
1489 }
1490
1491 inline
1492 DataArrayView::ValueType::size_type
1493 DataArrayView::relIndex(ValueType::size_type i) const
1494 {
1495 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1496 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1497 return i;
1498 }
1499
1500 inline
1501 DataArrayView::ValueType::size_type
1502 DataArrayView::index(ValueType::size_type i) const
1503 {
1504 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1505 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1506 return (m_offset+i);
1507 }
1508
1509 inline
1510 DataArrayView::ValueType::size_type
1511 DataArrayView::relIndex(ValueType::size_type i,
1512 ValueType::size_type j) const
1513 {
1514 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1515 ValueType::size_type temp=i+j*m_shape[0];
1516 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1517 return temp;
1518 }
1519
1520 inline
1521 DataArrayView::ValueType::size_type
1522 DataArrayView::index(ValueType::size_type i,
1523 ValueType::size_type j) const
1524 {
1525 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1526 ValueType::size_type temp=i+j*m_shape[0];
1527 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1528 return (m_offset+temp);
1529 }
1530
1531 inline
1532 DataArrayView::ValueType::size_type
1533 DataArrayView::relIndex(ValueType::size_type i,
1534 ValueType::size_type j,
1535 ValueType::size_type k) const
1536 {
1537 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1538 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1539 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1540 return temp;
1541 }
1542
1543 inline
1544 DataArrayView::ValueType::size_type
1545 DataArrayView::index(ValueType::size_type i,
1546 ValueType::size_type j,
1547 ValueType::size_type k) const
1548 {
1549 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1550 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1551 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1552 return (m_offset+temp);
1553 }
1554
1555 inline
1556 DataArrayView::ValueType::size_type
1557 DataArrayView::relIndex(ValueType::size_type i,
1558 ValueType::size_type j,
1559 ValueType::size_type k,
1560 ValueType::size_type m) const
1561 {
1562 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1563 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];
1564 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1565 return temp;
1566 }
1567
1568 inline
1569 DataArrayView::ValueType::size_type
1570 DataArrayView::index(ValueType::size_type i,
1571 ValueType::size_type j,
1572 ValueType::size_type k,
1573 ValueType::size_type m) const
1574 {
1575 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1576 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];
1577 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1578 return (m_offset+temp);
1579 }
1580
1581 inline
1582 DataArrayView::ValueType::reference
1583 DataArrayView::operator()()
1584 {
1585 return (*m_data)[index()];
1586 }
1587
1588 inline
1589 DataArrayView::ValueType::const_reference
1590 DataArrayView::operator()() const
1591 {
1592 return (*m_data)[index()];
1593 }
1594
1595 inline
1596 DataArrayView::ValueType::reference
1597 DataArrayView::operator()(ValueType::size_type i)
1598 {
1599 return (*m_data)[index(i)];
1600 }
1601
1602 inline
1603 DataArrayView::ValueType::const_reference
1604 DataArrayView::operator()(ValueType::size_type i) const
1605 {
1606 return (*m_data)[index(i)];
1607 }
1608
1609 inline
1610 DataArrayView::ValueType::reference
1611 DataArrayView::operator()(ValueType::size_type i,
1612 ValueType::size_type j)
1613 {
1614 return (*m_data)[index(i,j)];
1615 }
1616
1617 inline
1618 DataArrayView::ValueType::const_reference
1619 DataArrayView::operator()(ValueType::size_type i,
1620 ValueType::size_type j) const
1621 {
1622 return (*m_data)[index(i,j)];
1623 }
1624
1625 inline
1626 DataArrayView::ValueType::reference
1627 DataArrayView::operator()(ValueType::size_type i,
1628 ValueType::size_type j,
1629 ValueType::size_type k)
1630 {
1631 return (*m_data)[index(i,j,k)];
1632 }
1633
1634 inline
1635 DataArrayView::ValueType::const_reference
1636 DataArrayView::operator()(ValueType::size_type i,
1637 ValueType::size_type j,
1638 ValueType::size_type k) const
1639 {
1640 return (*m_data)[index(i,j,k)];
1641 }
1642
1643 inline
1644 DataArrayView::ValueType::reference
1645 DataArrayView::operator()(ValueType::size_type i,
1646 ValueType::size_type j,
1647 ValueType::size_type k,
1648 ValueType::size_type m)
1649 {
1650 return (*m_data)[index(i,j,k,m)];
1651 }
1652
1653 inline
1654 DataArrayView::ValueType::const_reference
1655 DataArrayView::operator()(ValueType::size_type i,
1656 ValueType::size_type j,
1657 ValueType::size_type k,
1658 ValueType::size_type m) const
1659 {
1660 return (*m_data)[index(i,j,k,m)];
1661 }
1662
1663 } // end of namespace
1664
1665 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26