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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (show annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years, 8 months ago) by robwdcock
File MIME type: text/plain
File size: 39267 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26