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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 950 - (show annotations)
Tue Feb 6 07:01:11 2007 UTC (11 years, 11 months ago) by gross
File MIME type: text/plain
File size: 54501 byte(s)
escript data objects can now be saved to netCDF files, see http://www.unidata.ucar.edu/software/netcdf/.
Currently only constant data are implemented with expanded and tagged data to follow.
There are two new functions to dump a data object

   s=Data(...)
   s.dump(<filename>)

and to recover it

   s=load(<filename>, domain)

Notice that the function space of s is recovered but domain is still need. 

dump and load will replace archive and extract.

The installation needs now the netCDF installed. 


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26