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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (show annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 4 months ago) by phornby
File MIME type: text/plain
File size: 55396 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26