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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 757 - (show annotations)
Mon Jun 26 13:12:56 2006 UTC (13 years, 5 months ago) by woo409
File MIME type: text/plain
File size: 40489 byte(s)
+ Merge of intelc_win32 branch (revision 741:755) with trunk. Tested on iVEC altix (run_tests and py_tests all pass)

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 "DataVector.h"
22 #include "LocalOps.h"
23
24 #include <boost/python/numeric.hpp>
25 #include <boost/python/object.hpp>
26
27 #include <vector>
28
29 namespace escript {
30
31 /**
32 \brief
33 DataArrayView provides a view of external data, configured according
34 to the given shape information and offset.
35
36 Description:
37 DataArrayView provides a view of data allocated externally. The
38 external data should provide sufficient values so that the dimensions
39 specified for the view shape will be satisfied. The viewer can update
40 values within the external data but cannot resize the external data.
41
42 The view provided represents a single n-dimensional data-point
43 comprised of values taken from the given data array, starting at the
44 specified offset and extending for as many values as are necessary to
45 satisfy the given shape. The default offset can be changed, or different
46 offsets specified, in order to provide views of other data-points in
47 the underlying data array.
48 */
49
50 class DataArrayView {
51
52 ESCRIPT_DLL_API friend bool operator==(const DataArrayView& left, const DataArrayView& right);
53 ESCRIPT_DLL_API friend bool operator!=(const DataArrayView& left, const DataArrayView& right);
54
55 public:
56
57 //
58 // Some basic types which define the data values and view shapes.
59 typedef DataVector ValueType;
60 //typedef std::vector<double> ValueType;
61 typedef std::vector<int> ShapeType;
62 typedef std::vector<std::pair<int, int> > RegionType;
63 typedef std::vector<std::pair<int, int> > RegionLoopRangeType;
64
65 /**
66 \brief
67 Default constructor for DataArrayView.
68
69 Description:
70 Default constructor for DataArrayView.
71 Creates an empty view with no associated data array.
72
73 This is essentially useless but here for completeness.
74 */
75 ESCRIPT_DLL_API
76 DataArrayView();
77
78 /**
79 \brief
80 Constructor for DataArrayView.
81
82 Description:
83 Constructor for DataArrayView.
84
85 \param data - Input -
86 Array holding data to be viewed. This must be large enough
87 to supply sufficient values for the specified shape starting at
88 the given offset.
89 \param viewShape - Input -
90 The shape of the view to create.
91 \param offset - Input -
92 Offset into the data at which view should start.
93 */
94 ESCRIPT_DLL_API
95 DataArrayView(ValueType& data,
96 const ShapeType& viewShape,
97 int offset=0);
98
99 /**
100 \brief
101 Copy constructor for DataArrayView.
102
103 Description:
104 Copy constructor for DataArrayView.
105
106 \param other - Input -
107 DataArrayView to copy.
108
109 NOTE: The copy references the same data array.
110 */
111 ESCRIPT_DLL_API
112 DataArrayView(const DataArrayView& other);
113
114 /**
115 \brief
116 Copy from a numarray into the data array viewed by this.
117 This must have same shape as the given value - will throw if not.
118 */
119 ESCRIPT_DLL_API
120 void
121 copy(const boost::python::numeric::array& value);
122
123 /**
124 \brief
125 Copy data from another DataArrayView into the data array viewed by this
126 starting at the default offset in both views.
127 The shapes of both views must be the same - will throw if not.
128 NB: views may or may not reference same underlying data array!
129 */
130 ESCRIPT_DLL_API
131 void
132 copy(const DataArrayView& other);
133
134 /**
135 \brief
136 Copy data from another DataArrayView into this view starting at the
137 given offset in this view and the default offset in the other view.
138 The shapes of both views must be the same - will throw if not.
139 NB: views may or may not reference same underlying data array!
140 */
141 ESCRIPT_DLL_API
142 void
143 copy(ValueType::size_type offset,
144 const DataArrayView& other);
145
146 /**
147 \brief
148 Copy data from another DataArrayView into this view starting at the
149 given offsets in each view.
150 The shapes of both views must be compatible - will throw if not.
151 NB: views may or may not reference same underlying data array!
152
153 \param thisOffset - Input -
154 Offset into this view's data array to copy to.
155 \param other - Input -
156 View to copy data from.
157 \param otherOffset - Input -
158 Offset into the other view's data array to copy from.
159 */
160 ESCRIPT_DLL_API
161 void
162 copy(ValueType::size_type thisOffset,
163 const DataArrayView& other,
164 ValueType::size_type otherOffset);
165
166 /**
167 \brief
168 Copy the given single value over the entire view starting at the given
169 offset in the view's data array.
170
171 \param offset - Input -
172 Offset into this view's data array to copy to.
173 \param value - Input -
174 Value to copy.
175 */
176 ESCRIPT_DLL_API
177 void
178 copy(ValueType::size_type offset,
179 ValueType::value_type value);
180
181 /**
182 \brief
183 Check if view is empty. ie: does not point to any actual data.
184 */
185 ESCRIPT_DLL_API
186 bool
187 isEmpty() const;
188
189 /**
190 \brief
191 Return this view's offset into the viewed data array.
192 */
193 ESCRIPT_DLL_API
194 ValueType::size_type
195 getOffset() const;
196
197 /**
198 \brief
199 Set the offset into the data array at which this view is to start.
200 This could be used to step through the underlying data array by incrementing
201 the offset by noValues successively. Thus this view would provide a moving
202 window on the underlying data with the given shape.
203 */
204 ESCRIPT_DLL_API
205 void
206 setOffset(ValueType::size_type offset);
207
208 /**
209 \brief
210 Increment the offset by the number of values in the shape, if possible. Thus
211 moving the view onto the next data point of the given shape in the underlying
212 data array.
213 */
214 ESCRIPT_DLL_API
215 void
216 incrOffset();
217
218 /**
219 \brief
220 Check the (given) offset will not result in two few values being available in
221 the underlying data array for this view's shape.
222 */
223 ESCRIPT_DLL_API
224 bool
225 checkOffset() const;
226
227 ESCRIPT_DLL_API
228 bool
229 checkOffset(ValueType::size_type offset) const;
230
231 /**
232 \brief
233 Return the rank of the shape of this view.
234 */
235 ESCRIPT_DLL_API
236 int
237 getRank() const;
238
239 /**
240 \brief
241 Return the number of values for the shape of this view.
242 */
243 ESCRIPT_DLL_API
244 int
245 noValues() const;
246
247 /**
248 \brief
249 Calculate the number of values for the given shape or region.
250 This is purely a utility method and has no bearing on this view.
251 */
252 ESCRIPT_DLL_API
253 static
254 int
255 noValues(const ShapeType& shape);
256
257 ESCRIPT_DLL_API
258 static
259 int
260 noValues(const RegionLoopRangeType& region);
261
262 /**
263 \brief
264 Return a reference to the underlying data array.
265 */
266 ESCRIPT_DLL_API
267 ValueType&
268 getData() const;
269
270 /**
271 \brief
272 Return a reference to the data value with the given
273 index in this view. This takes into account the offset.
274 */
275 ESCRIPT_DLL_API
276 ValueType::reference
277 getData(ValueType::size_type i) const;
278
279 /**
280 \brief
281 Return the shape of this view.
282 */
283 ESCRIPT_DLL_API
284 const
285 ShapeType&
286 getShape() const;
287
288 /**
289 \brief
290 Return true if the given shape is the same as this view's shape.
291 */
292 ESCRIPT_DLL_API
293 bool
294 checkShape(const ShapeType& other) const;
295
296 /**
297 \brief
298 Create a shape error message. Normally used when there is a shape
299 mismatch between this shape and the other shape.
300
301 \param messagePrefix - Input -
302 First part of the error message.
303 \param other - Input -
304 The other shape.
305 */
306 ESCRIPT_DLL_API
307 std::string
308 createShapeErrorMessage(const std::string& messagePrefix,
309 const ShapeType& other) const;
310
311 /**
312 \brief
313 Return the viewed data as a formatted string.
314 Not recommended for very large objects!
315
316 \param suffix - Input -
317 Suffix appended to index display.
318 */
319 ESCRIPT_DLL_API
320 std::string
321 toString(const std::string& suffix=std::string("")) const;
322
323 /**
324 \brief
325 Return the given shape as a string.
326 This is purely a utility method and has no bearing on this view.
327
328 \param shape - Input.
329 */
330 ESCRIPT_DLL_API
331 static
332 std::string
333 shapeToString(const ShapeType& shape);
334
335 /**
336 \brief
337 Return the 1 dimensional index into the data array of the only element
338 in the view, *ignoring the offset*.
339 Assumes a rank 0 view.
340 */
341 ESCRIPT_DLL_API
342 ValueType::size_type
343 relIndex() const;
344
345 /**
346 \brief
347 Return the 1 dimensional index into the data array of the element at
348 position i in the view, *ignoring the offset*.
349 Assumes a rank 1 view.
350 */
351 ESCRIPT_DLL_API
352 ValueType::size_type
353 relIndex(ValueType::size_type i) const;
354
355 /**
356 \brief
357 Return the 1 dimensional index into the data array of the element at
358 position i,j in the view, *ignoring the offset*.
359 Assumes a rank 2 view.
360 */
361 ESCRIPT_DLL_API
362 ValueType::size_type
363 relIndex(ValueType::size_type i,
364 ValueType::size_type j) const;
365
366 /**
367 \brief
368 Return the 1 dimensional index into the data array of the element at
369 position i,j,k in the view, *ignoring the offset*.
370 Assumes a rank 3 view.
371 */
372 ESCRIPT_DLL_API
373 ValueType::size_type
374 relIndex(ValueType::size_type i,
375 ValueType::size_type j,
376 ValueType::size_type k) const;
377
378 /**
379 \brief
380 Return the 1 dimensional index into the data array of the element at
381 position i,j,k,m in the view, *ignoring the offset*.
382 Assumes a rank 4 view.
383 */
384 ESCRIPT_DLL_API
385 ValueType::size_type
386 relIndex(ValueType::size_type i,
387 ValueType::size_type j,
388 ValueType::size_type k,
389 ValueType::size_type m) const;
390
391 /**
392 \brief
393 Return the 1 dimensional index into the data array of the only element
394 in the view.
395 Assumes a rank 0 view.
396 */
397 ESCRIPT_DLL_API
398 ValueType::size_type
399 index() const;
400
401 /**
402 \brief
403 Return the 1 dimensional index into the data array of the element at
404 position i in the view.
405 Assumes a rank 1 view.
406 */
407 ESCRIPT_DLL_API
408 ValueType::size_type
409 index(ValueType::size_type i) const;
410
411 /**
412 \brief
413 Return the 1 dimensional index into the data array of the element at
414 position i,j in the view.
415 Assumes a rank 2 view.
416 */
417 ESCRIPT_DLL_API
418 ValueType::size_type
419 index(ValueType::size_type i,
420 ValueType::size_type j) const;
421
422 /**
423 \brief
424 Return the 1 dimensional index into the data array of the element at
425 position i,j,k in the view.
426 Assumes a rank 3 view.
427 */
428 ESCRIPT_DLL_API
429 ValueType::size_type
430 index(ValueType::size_type i,
431 ValueType::size_type j,
432 ValueType::size_type k) const;
433
434 /**
435 \brief
436 Return the 1 dimensional index into the data array of the element at
437 position i,j,k,m in the view.
438 Assumes a rank 4 view.
439 */
440 ESCRIPT_DLL_API
441 ValueType::size_type
442 index(ValueType::size_type i,
443 ValueType::size_type j,
444 ValueType::size_type k,
445 ValueType::size_type m) const;
446
447 /**
448 \brief
449 Return a reference for the only element in the view.
450 Assumes a rank 0 view.
451 */
452 ESCRIPT_DLL_API
453 ValueType::reference
454 operator()();
455
456 ESCRIPT_DLL_API
457 ValueType::const_reference
458 operator()() const;
459
460 /**
461 \brief
462 Return a reference to the element at position i in the view.
463 Assumes a rank 1 view.
464 */
465 ESCRIPT_DLL_API
466 ValueType::reference
467 operator()(ValueType::size_type i);
468
469 ESCRIPT_DLL_API
470 ValueType::const_reference
471 operator()(ValueType::size_type i) const;
472
473 /**
474 \brief
475 Return a reference to the element at position i,j in the view.
476 Assumes a rank 2 view.
477 */
478 ESCRIPT_DLL_API
479 ValueType::reference
480 operator()(ValueType::size_type i,
481 ValueType::size_type j);
482
483 ESCRIPT_DLL_API
484 ValueType::const_reference
485 operator()(ValueType::size_type i,
486 ValueType::size_type j) const;
487
488 /**
489 \brief
490 Return a reference to the element at position i,j,k in the view.
491 Assumes a rank 3 view.
492 */
493 ESCRIPT_DLL_API
494 ValueType::reference
495 operator()(ValueType::size_type i,
496 ValueType::size_type j,
497 ValueType::size_type k);
498
499 ESCRIPT_DLL_API
500 ValueType::const_reference
501 operator()(ValueType::size_type i,
502 ValueType::size_type j,
503 ValueType::size_type k) const;
504
505 /**
506 \brief
507 Return a reference to the element at position i,j,k,m in the view.
508 Assumes a rank 4 view.
509 */
510 ESCRIPT_DLL_API
511 ValueType::reference
512 operator()(ValueType::size_type i,
513 ValueType::size_type j,
514 ValueType::size_type k,
515 ValueType::size_type m);
516
517 ESCRIPT_DLL_API
518 ValueType::const_reference
519 operator()(ValueType::size_type i,
520 ValueType::size_type j,
521 ValueType::size_type k,
522 ValueType::size_type m) const;
523
524 /**
525 \brief
526 Determine the shape of the specified slice region.
527 This is purely a utility method and has no bearing on this view.
528
529 \param region - Input -
530 Slice region.
531 */
532 ESCRIPT_DLL_API
533 static
534 ShapeType
535 getResultSliceShape(const RegionType& region);
536
537 /**
538 \brief
539 Determine the region specified by the given python slice object.
540
541 \param key - Input -
542 python slice object specifying region to be returned.
543
544 The slice object is a tuple of n python slice specifiers, where
545 n <= the rank of this Data object. Each slice specifier specifies the
546 range of indexes to be sliced from the corresponding dimension. The
547 first specifier corresponds to the first dimension, the second to the
548 second and so on. Where n < the rank, the remaining dimensions are
549 sliced across the full range of their indicies.
550
551 Each slice specifier is of the form "a:b", which specifies a slice
552 from index a, up to but not including index b. Where index a is ommitted
553 a is assumed to be 0. Where index b is ommitted, b is assumed to be the
554 length of this dimension. Where both are ommitted (eg: ":") the slice is
555 assumed to encompass that entire dimension.
556
557 Where one of the slice specifiers is a single integer, eg: [1], we
558 want to generate a rank-1 dimension object, as opposed to eg: [1,2]
559 which implies we want to take a rank dimensional object with one
560 dimension of size 1.
561
562 The return value is a vector of pairs with length equal to the rank of
563 this object. Each pair corresponds to the range of indicies from the
564 corresponding dimension to be sliced from, as specified in the input
565 slice object.
566
567 Examples:
568
569 For a rank 1 object of shape(5):
570
571 getSliceRegion(:) => < <0,5> >
572 getSliceRegion(2:3) => < <2,3> >
573 getSliceRegion(:3) => < <0,3> >
574 getSliceRegion(2:) => < <2,5> >
575
576 For a rank 2 object of shape(4,5):
577
578 getSliceRegion(2:3) => < <2,3> <0,5> >
579 getSliceRegion(2) => < <2,3> <0,5> >
580 NB: but return object requested will have rank 1, shape(5), with
581 values taken from index 2 of this object's first dimension.
582
583 For a rank 3 object of shape (2,4,6):
584
585 getSliceRegion(0:2,0:4,0:6) => < <0,2> <0,4> <0,6> >
586 getSliceRegion(:,:,:) => < <0,2> <0,4> <0,6> >
587 getSliceRegion(0:1) => < <0,1> <0,4> <0,6> >
588 getSliceRegion(:1,0:2) => < <0,1> <0,2> <0,6> >
589
590 */
591 ESCRIPT_DLL_API
592 RegionType
593 getSliceRegion(const boost::python::object& key) const;
594
595 /**
596 \brief
597 Copy a data slice specified by the given region from the given view
598 into this view, using the default offsets in both views.
599
600 \param other - Input -
601 View to copy data from.
602 \param region - Input -
603 Region in other view to copy data from.
604 */
605 ESCRIPT_DLL_API
606 void
607 copySlice(const DataArrayView& other,
608 const RegionLoopRangeType& region);
609
610 /**
611 \brief
612 Copy a data slice specified by the given region and offset from the
613 given view into this view at the given offset.
614
615 \param thisOffset - Input -
616 Copy the slice to this offset in this view.
617 \param other - Input -
618 View to copy data from.
619 \param otherOffset - Input -
620 Copy the slice from this offset in the given view.
621 \param region - Input -
622 Region in other view to copy data from.
623 */
624 ESCRIPT_DLL_API
625 void
626 copySlice(ValueType::size_type thisOffset,
627 const DataArrayView& other,
628 ValueType::size_type otherOffset,
629 const RegionLoopRangeType& region);
630
631 /**
632 \brief
633 Copy data into a slice specified by the given region in this view from
634 the given view, using the default offsets in both views.
635
636 \param other - Input -
637 View to copy data from.
638 \param region - Input -
639 Region in this view to copy data to.
640 */
641 ESCRIPT_DLL_API
642 void
643 copySliceFrom(const DataArrayView& other,
644 const RegionLoopRangeType& region);
645
646 /**
647 \brief
648 Copy data into a slice specified by the given region and offset in
649 this view from the given view at the given offset.
650
651 \param thisOffset - Input -
652 Copy the slice to this offset in this view.
653 \param other - Input -
654 View to copy data from.
655 \param otherOffset - Input -
656 Copy the slice from this offset in the given view.
657 \param region - Input -
658 Region in this view to copy data to.
659 */
660 ESCRIPT_DLL_API
661 void
662 copySliceFrom(ValueType::size_type thisOffset,
663 const DataArrayView& other,
664 ValueType::size_type otherOffset,
665 const RegionLoopRangeType& region);
666
667 /**
668 \brief
669 Perform the unary operation on the data point specified by the view's
670 default offset. Applies the specified operation to each value in the data
671 point.
672
673 Called by escript::unaryOp.
674
675 \param operation - Input -
676 Operation to apply. Operation must be a pointer to a function.
677 */
678 template <class UnaryFunction>
679 ESCRIPT_DLL_API
680 void
681 unaryOp(UnaryFunction operation);
682
683 /**
684 \brief
685 Perform the unary operation on the data point specified by the given
686 offset. Applies the specified operation to each value in the data
687 point. Operation must be a pointer to a function.
688
689 Called by escript::unaryOp.
690
691 \param offset - Input -
692 Apply the operation to data starting at this offset in this view.
693 \param operation - Input -
694 Operation to apply. Must be a pointer to a function.
695 */
696 template <class UnaryFunction>
697 ESCRIPT_DLL_API
698 void
699 unaryOp(ValueType::size_type offset,
700 UnaryFunction operation);
701
702 /**
703 \brief
704 Perform the binary operation on the data points specified by the default
705 offsets in this view and in view "right". Applies the specified operation
706 to corresponding values in both data points. Operation must be a pointer
707 to a function.
708
709 Called by escript::binaryOp.
710
711 \param right - Input -
712 View to act as RHS in given binary operation.
713 \param operation - Input -
714 Operation to apply. Must be a pointer to a function.
715 */
716 template <class BinaryFunction>
717 ESCRIPT_DLL_API
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 ESCRIPT_DLL_API
742 void
743 binaryOp(ValueType::size_type leftOffset,
744 const DataArrayView& right,
745 ValueType::size_type rightOffset,
746 BinaryFunction operation);
747
748 /**
749 \brief
750 Perform the binary operation on the data point specified by the default
751 offset in this view using the scalar value "right". Applies the specified
752 operation to values in the data point. Operation must be a pointer to
753 a function.
754
755 Called by escript::binaryOp.
756
757 \param right - Input -
758 Value to act as RHS in given binary operation.
759 \param operation - Input -
760 Operation to apply. Must be a pointer to a function.
761 */
762 template <class BinaryFunction>
763 ESCRIPT_DLL_API
764 void
765 binaryOp(double right,
766 BinaryFunction operation);
767
768 /**
769 \brief
770 Perform the binary operation on the data point specified by the given
771 offset in this view using the scalar value "right". Applies the specified
772 operation to values in the data point. Operation must be a pointer
773 to a function.
774
775 Called by escript::binaryOp.
776
777 \param offset - Input -
778 Apply the operation to data starting at this offset in this view.
779 \param right - Input -
780 Value to act as RHS in given binary operation.
781 \param operation - Input -
782 Operation to apply. Must be a pointer to a function.
783 */
784 template <class BinaryFunction>
785 ESCRIPT_DLL_API
786 void
787 binaryOp(ValueType::size_type offset,
788 double right,
789 BinaryFunction operation);
790
791 /**
792 \brief
793 Perform the given data point reduction operation on the data point
794 specified by the default offset into the view. Reduces all elements of
795 the data point using the given operation, returning the result as a
796 scalar. Operation must be a pointer to a function.
797
798 Called by escript::algorithm.
799
800 \param operation - Input -
801 Operation to apply. Must be a pointer to a function.
802 */
803 template <class BinaryFunction>
804 ESCRIPT_DLL_API
805 double
806 reductionOp(BinaryFunction operation,
807 double initial_value) const;
808
809 /**
810 \brief
811 Perform the given data point reduction operation on the data point
812 specified by the given offset into the view. Reduces all elements of
813 the data point using the given operation, returning the result as a
814 scalar. Operation must be a pointer to a function.
815
816 Called by escript::algorithm.
817
818 \param offset - Input -
819 Apply the operation to data starting at this offset in this view.
820 \param operation - Input -
821 Operation to apply. Must be a pointer to a function.
822 */
823 template <class BinaryFunction>
824 ESCRIPT_DLL_API
825 double
826 reductionOp(ValueType::size_type offset,
827 BinaryFunction operation,
828 double initial_value) const;
829
830 /**
831 \brief
832 Perform a matrix multiply of the given views.
833 This is purely a utility method and has no bearing on this view.
834
835 NB: Only multiplies together the two given views, ie: two data-points,
836 would need to call this over all data-points to multiply the entire
837 Data objects involved.
838
839 \param left - Input - The left hand side.
840 \param right - Input - The right hand side.
841 \param result - Output - The result of the operation.
842 */
843 ESCRIPT_DLL_API
844 static
845 void
846 matMult(const DataArrayView& left,
847 const DataArrayView& right,
848 DataArrayView& result);
849
850 /**
851 \brief
852 Determine the shape of the result array for a matrix multiplication
853 of the given views.
854 This is purely a utility method and has no bearing on this view.
855 */
856 ESCRIPT_DLL_API
857 static
858 ShapeType
859 determineResultShape(const DataArrayView& left,
860 const DataArrayView& right);
861
862 /**
863 \brief
864 solves a local eigenvalue problem
865
866 \param in - Input - matrix
867 \param inOffset - Input - offset into in
868 \param ev - Output - The eigenvalues
869 \param inOffset - Input - offset into ev
870 */
871 ESCRIPT_DLL_API
872 static
873 inline
874 void
875 eigenvalues(DataArrayView& in,
876 ValueType::size_type inOffset,
877 DataArrayView& ev,
878 ValueType::size_type evOffset)
879 {
880 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
881 double ev0,ev1,ev2;
882 int s=in.getShape()[0];
883 if (s==1) {
884 in00=(*(in.m_data))[inOffset+in.index(0,0)];
885 eigenvalues1(in00,&ev0);
886 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
887
888 } else if (s==2) {
889 in00=(*(in.m_data))[inOffset+in.index(0,0)];
890 in10=(*(in.m_data))[inOffset+in.index(1,0)];
891 in01=(*(in.m_data))[inOffset+in.index(0,1)];
892 in11=(*(in.m_data))[inOffset+in.index(1,1)];
893 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
894 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
895 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
896
897 } else if (s==3) {
898 in00=(*(in.m_data))[inOffset+in.index(0,0)];
899 in10=(*(in.m_data))[inOffset+in.index(1,0)];
900 in20=(*(in.m_data))[inOffset+in.index(2,0)];
901 in01=(*(in.m_data))[inOffset+in.index(0,1)];
902 in11=(*(in.m_data))[inOffset+in.index(1,1)];
903 in21=(*(in.m_data))[inOffset+in.index(2,1)];
904 in02=(*(in.m_data))[inOffset+in.index(0,2)];
905 in12=(*(in.m_data))[inOffset+in.index(1,2)];
906 in22=(*(in.m_data))[inOffset+in.index(2,2)];
907 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
908 &ev0,&ev1,&ev2);
909 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
910 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
911 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
912
913 }
914 }
915
916 /**
917 \brief
918 solves a local eigenvalue problem
919
920 \param in - Input - matrix
921 \param inOffset - Input - offset into in
922 \param ev - Output - The eigenvalues
923 \param evOffset - Input - offset into ev
924 \param V - Output - The eigenvectors
925 \param VOffset - Input - offset into V
926 \param tol - Input - eigenvalues with relative difference tol are treated as equal
927 */
928 ESCRIPT_DLL_API
929 static
930 inline
931 void
932 eigenvalues_and_eigenvectors(DataArrayView& in,
933 ValueType::size_type inOffset,
934 DataArrayView& ev,
935 ValueType::size_type evOffset,
936 DataArrayView& V,
937 ValueType::size_type VOffset,
938 const double tol=1.e-13)
939 {
940 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
941 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
942 double ev0,ev1,ev2;
943 int s=in.getShape()[0];
944 if (s==1) {
945 in00=(*(in.m_data))[inOffset+in.index(0,0)];
946 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
947 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
948 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
949 } else if (s==2) {
950 in00=(*(in.m_data))[inOffset+in.index(0,0)];
951 in10=(*(in.m_data))[inOffset+in.index(1,0)];
952 in01=(*(in.m_data))[inOffset+in.index(0,1)];
953 in11=(*(in.m_data))[inOffset+in.index(1,1)];
954 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
955 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
956 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
957 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
958 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
959 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
960 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
961 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
962 } else if (s==3) {
963 in00=(*(in.m_data))[inOffset+in.index(0,0)];
964 in10=(*(in.m_data))[inOffset+in.index(1,0)];
965 in20=(*(in.m_data))[inOffset+in.index(2,0)];
966 in01=(*(in.m_data))[inOffset+in.index(0,1)];
967 in11=(*(in.m_data))[inOffset+in.index(1,1)];
968 in21=(*(in.m_data))[inOffset+in.index(2,1)];
969 in02=(*(in.m_data))[inOffset+in.index(0,2)];
970 in12=(*(in.m_data))[inOffset+in.index(1,2)];
971 in22=(*(in.m_data))[inOffset+in.index(2,2)];
972 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
973 &ev0,&ev1,&ev2,
974 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
975 (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
976 (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
977 (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
978 (*(V.m_data))[inOffset+V.index(0,0)]=V00;
979 (*(V.m_data))[inOffset+V.index(1,0)]=V10;
980 (*(V.m_data))[inOffset+V.index(2,0)]=V20;
981 (*(V.m_data))[inOffset+V.index(0,1)]=V01;
982 (*(V.m_data))[inOffset+V.index(1,1)]=V11;
983 (*(V.m_data))[inOffset+V.index(2,1)]=V21;
984 (*(V.m_data))[inOffset+V.index(0,2)]=V02;
985 (*(V.m_data))[inOffset+V.index(1,2)]=V12;
986 (*(V.m_data))[inOffset+V.index(2,2)]=V22;
987
988 }
989 }
990 protected:
991
992 private:
993
994 //
995 // The maximum rank allowed for the shape of any view.
996 static const int m_maxRank=4;
997
998 //
999 // The data values for the view.
1000 // NOTE: This points to data external to the view.
1001 // This is just a pointer to an array of ValueType.
1002 ValueType* m_data;
1003
1004 //
1005 // The offset into the data array used by different views.
1006 // This is simply an integer specifying a position in the data array
1007 // pointed to by m_data.
1008 ValueType::size_type m_offset;
1009
1010 //
1011 // The shape of the data.
1012 // This is simply an STL vector specifying the lengths of each dimension
1013 // of the shape as ints.
1014 ShapeType m_shape;
1015
1016 //
1017 // The number of values needed for the array.
1018 // This can be derived from m_shape by multiplying the size of each dimension, but
1019 // is stored here for convenience.
1020 int m_noValues;
1021
1022 };
1023
1024 ESCRIPT_DLL_API bool operator==(const DataArrayView& left, const DataArrayView& right);
1025 ESCRIPT_DLL_API bool operator!=(const DataArrayView& left, const DataArrayView& right);
1026
1027 /**
1028 \brief
1029 Modify region to copy from in order to
1030 deal with the case where one range in the region contains identical indexes,
1031 eg: <<1,1><0,3><0,3>>
1032 This situation implies we want to copy from an object with rank greater than that of this
1033 object. eg: we want to copy the values from a two dimensional slice out of a three
1034 dimensional object into a two dimensional object.
1035 We do this by taking a slice from the other object where one dimension of
1036 the slice region is of size 1. So in the above example, we modify the above
1037 region like so: <<1,2><0,3><0,3>> and take this slice.
1038 */
1039 DataArrayView::RegionLoopRangeType
1040 getSliceRegionLoopRange(const DataArrayView::RegionType& region);
1041
1042 /**
1043 \brief
1044 Calculate the slice range from the given python key object
1045 Used by DataArrayView::getSliceRegion.
1046 Returns the python slice object key as a pair of ints where the first
1047 member is the start and the second member is the end. the presence of a possible
1048 step attribute with value different from one will throw an exception
1049
1050 /param key - Input - key object specifying slice range.
1051 */
1052 std::pair<int,int>
1053 getSliceRange(const boost::python::object& key,
1054 const int shape);
1055
1056 /**
1057 Inline function definitions.
1058 */
1059
1060 template <class UnaryFunction>
1061 inline
1062 void
1063 DataArrayView::unaryOp(UnaryFunction operation)
1064 {
1065 unaryOp(m_offset,operation);
1066 }
1067
1068 template <class UnaryFunction>
1069 inline
1070 void
1071 DataArrayView::unaryOp(ValueType::size_type offset,
1072 UnaryFunction operation)
1073 {
1074 EsysAssert((!isEmpty()&&checkOffset(offset)),
1075 "Error - Couldn't perform unaryOp due to insufficient storage.");
1076 for (ValueType::size_type i=0;i<noValues();i++) {
1077 (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1078 }
1079 }
1080
1081 template <class BinaryFunction>
1082 inline
1083 void
1084 DataArrayView::binaryOp(const DataArrayView& right,
1085 BinaryFunction operation)
1086 {
1087 binaryOp(m_offset,right,right.getOffset(),operation);
1088 }
1089
1090 template <class BinaryFunction>
1091 inline
1092 void
1093 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1094 const DataArrayView& right,
1095 ValueType::size_type rightOffset,
1096 BinaryFunction operation)
1097 {
1098 EsysAssert(getShape()==right.getShape(),
1099 "Error - Couldn't perform binaryOp due to shape mismatch,");
1100 EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1101 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1102 EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1103 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1104 for (ValueType::size_type i=0;i<noValues();i++) {
1105 (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1106 }
1107 }
1108
1109 template <class BinaryFunction>
1110 inline
1111 void
1112 DataArrayView::binaryOp(double right,
1113 BinaryFunction operation)
1114 {
1115 binaryOp(m_offset,right,operation);
1116 }
1117
1118 template <class BinaryFunction>
1119 inline
1120 void
1121 DataArrayView::binaryOp(ValueType::size_type offset,
1122 double right,
1123 BinaryFunction operation)
1124 {
1125 EsysAssert((!isEmpty()&&checkOffset(offset)),
1126 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1127 for (ValueType::size_type i=0;i<noValues();i++) {
1128 (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1129 }
1130 }
1131
1132 template <class BinaryFunction>
1133 inline
1134 double
1135 DataArrayView::reductionOp(BinaryFunction operation,
1136 double initial_value) const
1137 {
1138 return reductionOp(m_offset,operation,initial_value);
1139 }
1140
1141 template <class BinaryFunction>
1142 inline
1143 double
1144 DataArrayView::reductionOp(ValueType::size_type offset,
1145 BinaryFunction operation,
1146 double initial_value) const
1147 {
1148 EsysAssert((!isEmpty()&&checkOffset(offset)),
1149 "Error - Couldn't perform reductionOp due to insufficient storage.");
1150 double current_value=initial_value;
1151 for (ValueType::size_type i=0;i<noValues();i++) {
1152 current_value=operation(current_value,(*m_data)[offset+i]);
1153 }
1154 return current_value;
1155 }
1156
1157 inline
1158 DataArrayView::ValueType::size_type
1159 DataArrayView::relIndex() const
1160 {
1161 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1162 return 0;
1163 }
1164
1165 inline
1166 DataArrayView::ValueType::size_type
1167 DataArrayView::index() const
1168 {
1169 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1170 return (m_offset);
1171 }
1172
1173 inline
1174 DataArrayView::ValueType::size_type
1175 DataArrayView::relIndex(ValueType::size_type i) const
1176 {
1177 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1178 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1179 return i;
1180 }
1181
1182 inline
1183 DataArrayView::ValueType::size_type
1184 DataArrayView::index(ValueType::size_type i) const
1185 {
1186 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1187 EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1188 return (m_offset+i);
1189 }
1190
1191 inline
1192 DataArrayView::ValueType::size_type
1193 DataArrayView::relIndex(ValueType::size_type i,
1194 ValueType::size_type j) const
1195 {
1196 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1197 ValueType::size_type temp=i+j*m_shape[0];
1198 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1199 return temp;
1200 }
1201
1202 inline
1203 DataArrayView::ValueType::size_type
1204 DataArrayView::index(ValueType::size_type i,
1205 ValueType::size_type j) const
1206 {
1207 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1208 ValueType::size_type temp=i+j*m_shape[0];
1209 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1210 return (m_offset+temp);
1211 }
1212
1213 inline
1214 DataArrayView::ValueType::size_type
1215 DataArrayView::relIndex(ValueType::size_type i,
1216 ValueType::size_type j,
1217 ValueType::size_type k) const
1218 {
1219 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1220 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1221 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1222 return temp;
1223 }
1224
1225 inline
1226 DataArrayView::ValueType::size_type
1227 DataArrayView::index(ValueType::size_type i,
1228 ValueType::size_type j,
1229 ValueType::size_type k) const
1230 {
1231 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1232 ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1233 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1234 return (m_offset+temp);
1235 }
1236
1237 inline
1238 DataArrayView::ValueType::size_type
1239 DataArrayView::relIndex(ValueType::size_type i,
1240 ValueType::size_type j,
1241 ValueType::size_type k,
1242 ValueType::size_type m) const
1243 {
1244 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1245 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];
1246 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1247 return temp;
1248 }
1249
1250 inline
1251 DataArrayView::ValueType::size_type
1252 DataArrayView::index(ValueType::size_type i,
1253 ValueType::size_type j,
1254 ValueType::size_type k,
1255 ValueType::size_type m) const
1256 {
1257 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1258 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];
1259 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1260 return (m_offset+temp);
1261 }
1262
1263 inline
1264 DataArrayView::ValueType::reference
1265 DataArrayView::operator()()
1266 {
1267 return (*m_data)[index()];
1268 }
1269
1270 inline
1271 DataArrayView::ValueType::const_reference
1272 DataArrayView::operator()() const
1273 {
1274 return (*m_data)[index()];
1275 }
1276
1277 inline
1278 DataArrayView::ValueType::reference
1279 DataArrayView::operator()(ValueType::size_type i)
1280 {
1281 return (*m_data)[index(i)];
1282 }
1283
1284 inline
1285 DataArrayView::ValueType::const_reference
1286 DataArrayView::operator()(ValueType::size_type i) const
1287 {
1288 return (*m_data)[index(i)];
1289 }
1290
1291 inline
1292 DataArrayView::ValueType::reference
1293 DataArrayView::operator()(ValueType::size_type i,
1294 ValueType::size_type j)
1295 {
1296 return (*m_data)[index(i,j)];
1297 }
1298
1299 inline
1300 DataArrayView::ValueType::const_reference
1301 DataArrayView::operator()(ValueType::size_type i,
1302 ValueType::size_type j) const
1303 {
1304 return (*m_data)[index(i,j)];
1305 }
1306
1307 inline
1308 DataArrayView::ValueType::reference
1309 DataArrayView::operator()(ValueType::size_type i,
1310 ValueType::size_type j,
1311 ValueType::size_type k)
1312 {
1313 return (*m_data)[index(i,j,k)];
1314 }
1315
1316 inline
1317 DataArrayView::ValueType::const_reference
1318 DataArrayView::operator()(ValueType::size_type i,
1319 ValueType::size_type j,
1320 ValueType::size_type k) const
1321 {
1322 return (*m_data)[index(i,j,k)];
1323 }
1324
1325 inline
1326 DataArrayView::ValueType::reference
1327 DataArrayView::operator()(ValueType::size_type i,
1328 ValueType::size_type j,
1329 ValueType::size_type k,
1330 ValueType::size_type m)
1331 {
1332 return (*m_data)[index(i,j,k,m)];
1333 }
1334
1335 inline
1336 DataArrayView::ValueType::const_reference
1337 DataArrayView::operator()(ValueType::size_type i,
1338 ValueType::size_type j,
1339 ValueType::size_type k,
1340 ValueType::size_type m) const
1341 {
1342 return (*m_data)[index(i,j,k,m)];
1343 }
1344
1345 } // end of namespace
1346
1347 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26