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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26