/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataArrayView.h
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/DataArrayView.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1714 - (show annotations)
Thu Aug 21 00:01:55 2008 UTC (12 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 57620 byte(s)
Branch commit

Moved getSliceRegion() and getSliceRange() into DataTypes
Data.cpp - modified not to rely on operator() from DataArrayView
         - Used more const& to avoid copies


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26