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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 682 - (hide 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 jgs 102 // $Id$
2 jgs 108
3 jgs 82 /*
4 elspeth 615 ************************************************************
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 jgs 82 */
14 jgs 474
15 jgs 82 #if !defined escript_DataArrayView_20040323_H
16     #define escript_DataArrayView_20040323_H
17    
18 robwdcock 682 #include "esysUtils/EsysAssert.h"
19 jgs 82
20 jgs 474 #include "DataVector.h"
21 gross 576 #include "LocalOps.h"
22 jgs 474
23 jgs 82 #include <boost/python/numeric.hpp>
24     #include <boost/python/object.hpp>
25    
26 jgs 108 #include <vector>
27 jgs 102
28 jgs 82 namespace escript {
29    
30     /**
31     \brief
32 jgs 108 DataArrayView provides a view of external data, configured according
33     to the given shape information and offset.
34 jgs 82
35     Description:
36     DataArrayView provides a view of data allocated externally. The
37     external data should provide sufficient values so that the dimensions
38 jgs 108 specified for the view shape will be satisfied. The viewer can update
39 jgs 82 values within the external data but cannot resize the external data.
40 jgs 108
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 jgs 117 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 jgs 82 */
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 jgs 117 //
57     // Some basic types which define the data values and view shapes.
58     typedef DataVector ValueType;
59     //typedef std::vector<double> ValueType;
60 jgs 108 typedef std::vector<int> ShapeType;
61 jgs 82 typedef std::vector<std::pair<int, int> > RegionType;
62 jgs 108 typedef std::vector<std::pair<int, int> > RegionLoopRangeType;
63 jgs 82
64     /**
65     \brief
66     Default constructor for DataArrayView.
67    
68     Description:
69 jgs 108 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 jgs 82 */
74     DataArrayView();
75    
76     /**
77     \brief
78     Constructor for DataArrayView.
79    
80     Description:
81     Constructor for DataArrayView.
82    
83 jgs 108 \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 jgs 82 */
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 jgs 108 \param other - Input -
104     DataArrayView to copy.
105 jgs 82
106 jgs 108 NOTE: The copy references the same data array.
107 jgs 82 */
108     DataArrayView(const DataArrayView& other);
109    
110     /**
111     \brief
112 jgs 108 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 jgs 82 */
115     void
116     copy(const boost::python::numeric::array& value);
117    
118     /**
119     \brief
120 jgs 108 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 jgs 82 */
125     void
126     copy(const DataArrayView& other);
127    
128     /**
129     \brief
130 jgs 108 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 jgs 82 */
135     void
136     copy(ValueType::size_type offset,
137     const DataArrayView& other);
138    
139     /**
140     \brief
141 jgs 108 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 jgs 82
146 jgs 108 \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 jgs 82 */
153     void
154     copy(ValueType::size_type thisOffset,
155     const DataArrayView& other,
156     ValueType::size_type otherOffset);
157    
158     /**
159     \brief
160 jgs 108 Copy the given single value over the entire view starting at the given
161     offset in the view's data array.
162 jgs 82
163 jgs 108 \param offset - Input -
164     Offset into this view's data array to copy to.
165     \param value - Input -
166     Value to copy.
167 jgs 82 */
168     void
169     copy(ValueType::size_type offset,
170     ValueType::value_type value);
171    
172     /**
173 jgs 108 \brief
174     Check if view is empty. ie: does not point to any actual data.
175     */
176     bool
177     isEmpty() const;
178    
179     /**
180 jgs 82 \brief
181 jgs 108 Return this view's offset into the viewed data array.
182 jgs 82 */
183     ValueType::size_type
184     getOffset() const;
185    
186     /**
187 jgs 108 \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 jgs 82 */
193 jgs 108 void
194     setOffset(ValueType::size_type offset);
195 jgs 82
196     /**
197     \brief
198 jgs 108 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 jgs 82 */
202 jgs 108 void
203     incrOffset();
204 jgs 82
205     /**
206     \brief
207 jgs 108 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 jgs 82 */
210 jgs 108 bool
211     checkOffset() const;
212 jgs 117
213 jgs 108 bool
214     checkOffset(ValueType::size_type offset) const;
215 jgs 82
216     /**
217 jgs 108 \brief
218     Return the rank of the shape of this view.
219     */
220     int
221     getRank() const;
222    
223     /**
224 jgs 82 \brief
225 jgs 108 Return the number of values for the shape of this view.
226 jgs 82 */
227     int
228     noValues() const;
229    
230     /**
231     \brief
232 jgs 108 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 jgs 82 */
235 jgs 117 static
236     int
237 jgs 108 noValues(const ShapeType& shape);
238 jgs 117
239     static
240     int
241 jgs 108 noValues(const RegionLoopRangeType& region);
242 jgs 82
243     /**
244     \brief
245 jgs 108 Return a reference to the underlying data array.
246 jgs 82 */
247     ValueType&
248     getData() const;
249    
250     /**
251     \brief
252 jgs 108 Return a reference to the data value with the given
253     index in this view. This takes into account the offset.
254 jgs 82 */
255     ValueType::reference
256     getData(ValueType::size_type i) const;
257    
258     /**
259     \brief
260 jgs 108 Return the shape of this view.
261 jgs 82 */
262 jgs 117 const
263     ShapeType&
264 jgs 108 getShape() const;
265 jgs 82
266     /**
267     \brief
268 jgs 108 Return true if the given shape is the same as this view's shape.
269 jgs 82 */
270 jgs 108 bool
271 jgs 117 checkShape(const ShapeType& other) const;
272 jgs 82
273     /**
274     \brief
275 jgs 108 Create a shape error message. Normally used when there is a shape
276 jgs 117 mismatch between this shape and the other shape.
277 jgs 82
278 jgs 108 \param messagePrefix - Input -
279     First part of the error message.
280     \param other - Input -
281     The other shape.
282 jgs 82 */
283 jgs 108 std::string
284     createShapeErrorMessage(const std::string& messagePrefix,
285 jgs 117 const ShapeType& other) const;
286 jgs 82
287     /**
288     \brief
289 jgs 108 Return the viewed data as a formatted string.
290     Not recommended for very large objects!
291 jgs 82
292 jgs 108 \param suffix - Input -
293     Suffix appended to index display.
294 jgs 82 */
295 jgs 108 std::string
296     toString(const std::string& suffix=std::string("")) const;
297 jgs 82
298     /**
299     \brief
300 jgs 108 Return the given shape as a string.
301     This is purely a utility method and has no bearing on this view.
302 jgs 82
303 jgs 108 \param shape - Input.
304 jgs 82 */
305 jgs 117 static
306     std::string
307     shapeToString(const ShapeType& shape);
308 jgs 82
309 jgs 108 /**
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 jgs 82 */
315 jgs 108 ValueType::size_type
316     relIndex() const;
317 jgs 82
318     /**
319 jgs 108 \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 jgs 82 */
324 jgs 108 ValueType::size_type
325     relIndex(ValueType::size_type i) const;
326 jgs 82
327     /**
328 jgs 108 \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 jgs 82 */
333 jgs 108 ValueType::size_type
334     relIndex(ValueType::size_type i,
335     ValueType::size_type j) const;
336 jgs 82
337     /**
338 jgs 108 \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 jgs 82 */
343 jgs 108 ValueType::size_type
344     relIndex(ValueType::size_type i,
345     ValueType::size_type j,
346     ValueType::size_type k) const;
347 jgs 82
348     /**
349     \brief
350 jgs 108 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 jgs 82 */
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 jgs 108 Return the 1 dimensional index into the data array of the only element
363     in the view.
364     Assumes a rank 0 view.
365 jgs 82 */
366     ValueType::size_type
367 jgs 108 index() const;
368 jgs 82
369     /**
370     \brief
371 jgs 108 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 jgs 82 */
375     ValueType::size_type
376 jgs 108 index(ValueType::size_type i) const;
377 jgs 82
378     /**
379     \brief
380 jgs 108 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 jgs 82 */
384     ValueType::size_type
385     index(ValueType::size_type i,
386 jgs 108 ValueType::size_type j) const;
387 jgs 82
388     /**
389     \brief
390 jgs 108 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 jgs 82 */
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 jgs 108 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 jgs 82 */
405     ValueType::size_type
406     index(ValueType::size_type i,
407 jgs 108 ValueType::size_type j,
408     ValueType::size_type k,
409     ValueType::size_type m) const;
410 jgs 82
411     /**
412     \brief
413 jgs 108 Return a reference for the only element in the view.
414     Assumes a rank 0 view.
415 jgs 82 */
416 jgs 108 ValueType::reference
417     operator()();
418 jgs 82
419 jgs 108 ValueType::const_reference
420     operator()() const;
421    
422 jgs 82 /**
423     \brief
424 jgs 108 Return a reference to the element at position i in the view.
425     Assumes a rank 1 view.
426 jgs 82 */
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 jgs 108 Return a reference to the element at position i,j in the view.
436     Assumes a rank 2 view.
437 jgs 82 */
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 jgs 108 Return a reference to the element at position i,j,k in the view.
449     Assumes a rank 3 view.
450 jgs 82 */
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 jgs 108 Return a reference to the element at position i,j,k,m in the view.
464     Assumes a rank 4 view.
465 jgs 82 */
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 jgs 108 \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 jgs 82 */
486 jgs 117 static
487     ShapeType
488 jgs 108 getResultSliceShape(const RegionType& region);
489 jgs 82
490 jgs 108 /**
491     \brief
492     Determine the region specified by the given python slice object.
493 jgs 82
494 jgs 108 \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 jgs 117 RegionType
545 jgs 108 getSliceRegion(const boost::python::object& key) const;
546    
547 jgs 82 /**
548     \brief
549 jgs 108 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 jgs 82 */
557     void
558 jgs 108 copySlice(const DataArrayView& other,
559     const RegionLoopRangeType& region);
560 jgs 82
561     /**
562     \brief
563 jgs 108 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 jgs 82 */
575     void
576 jgs 108 copySlice(ValueType::size_type thisOffset,
577     const DataArrayView& other,
578     ValueType::size_type otherOffset,
579     const RegionLoopRangeType& region);
580 jgs 82
581     /**
582     \brief
583 jgs 108 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 jgs 102
586 jgs 108 \param other - Input -
587     View to copy data from.
588     \param region - Input -
589     Region in this view to copy data to.
590 jgs 102 */
591 jgs 108 void
592     copySliceFrom(const DataArrayView& other,
593     const RegionLoopRangeType& region);
594 jgs 102
595     /**
596     \brief
597 jgs 108 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 jgs 102
600 jgs 108 \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 jgs 102 */
609 jgs 108 void
610     copySliceFrom(ValueType::size_type thisOffset,
611     const DataArrayView& other,
612     ValueType::size_type otherOffset,
613     const RegionLoopRangeType& region);
614 jgs 102
615     /**
616     \brief
617 jgs 108 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 jgs 82 */
626     template <class UnaryFunction>
627 jgs 108 void
628     unaryOp(UnaryFunction operation);
629 jgs 82
630     /**
631     \brief
632 jgs 108 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 jgs 82 */
643     template <class UnaryFunction>
644 jgs 108 void
645     unaryOp(ValueType::size_type offset,
646     UnaryFunction operation);
647 jgs 82
648     /**
649     \brief
650 jgs 108 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 jgs 82 */
662     template <class BinaryFunction>
663     void
664 jgs 108 binaryOp(const DataArrayView& right,
665 jgs 82 BinaryFunction operation);
666    
667     /**
668     \brief
669 jgs 108 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 jgs 82 */
685     template <class BinaryFunction>
686     void
687 jgs 108 binaryOp(ValueType::size_type leftOffset,
688     const DataArrayView& right,
689     ValueType::size_type rightOffset,
690 jgs 82 BinaryFunction operation);
691    
692     /**
693     \brief
694 jgs 108 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 jgs 82 */
706     template <class BinaryFunction>
707     void
708 jgs 108 binaryOp(double right,
709 jgs 82 BinaryFunction operation);
710    
711     /**
712     \brief
713 jgs 108 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 jgs 82 */
727     template <class BinaryFunction>
728     void
729 jgs 108 binaryOp(ValueType::size_type offset,
730     double right,
731 jgs 82 BinaryFunction operation);
732    
733     /**
734     \brief
735 jgs 108 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 jgs 147 template <class BinaryFunction>
746 jgs 108 double
747 jgs 147 reductionOp(BinaryFunction operation,
748     double initial_value) const;
749 jgs 108
750     /**
751     \brief
752 jgs 113 Perform the given data point reduction operation on the data point
753 jgs 108 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 jgs 147 template <class BinaryFunction>
765 jgs 108 double
766     reductionOp(ValueType::size_type offset,
767 jgs 147 BinaryFunction operation,
768     double initial_value) const;
769 jgs 108
770 gross 584 /**
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 jgs 108 /**
790     \brief
791 gross 584 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 gross 576 solves a local eigenvalue problem
803    
804 gross 580 \param in - Input - matrix
805     \param inOffset - Input - offset into in
806 gross 576 \param ev - Output - The eigenvalues
807 gross 580 \param inOffset - Input - offset into ev
808 gross 576 */
809 gross 580 static
810 gross 576 inline
811     void
812 gross 584 eigenvalues(DataArrayView& in,
813     ValueType::size_type inOffset,
814     DataArrayView& ev,
815     ValueType::size_type evOffset)
816     {
817 gross 580 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 gross 576
825 gross 580 } 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 gross 576
834 gross 580 } 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 gross 584 }
852 gross 580
853 gross 576 /**
854     \brief
855 gross 580 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 gross 584 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 gross 580 {
876 gross 583 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
877     double V00,V10,V20,V01,V11,V21,V02,V12,V22;
878 gross 580 double ev0,ev1,ev2;
879     int s=in.getShape()[0];
880     if (s==1) {
881 gross 583 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 gross 580 } else if (s==2) {
886 gross 583 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 gross 580 } else if (s==3) {
899 gross 583 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 gross 580
924     }
925     }
926 jgs 82 protected:
927    
928     private:
929    
930     //
931 jgs 108 // The maximum rank allowed for the shape of any view.
932 jgs 82 static const int m_maxRank=4;
933    
934     //
935 jgs 108 // The data values for the view.
936     // NOTE: This points to data external to the view.
937 jgs 117 // This is just a pointer to an array of ValueType.
938 jgs 82 ValueType* m_data;
939    
940     //
941     // The offset into the data array used by different views.
942 jgs 108 // This is simply an integer specifying a position in the data array
943     // pointed to by m_data.
944 jgs 82 ValueType::size_type m_offset;
945    
946     //
947     // The shape of the data.
948 jgs 108 // This is simply an STL vector specifying the lengths of each dimension
949     // of the shape as ints.
950 jgs 82 ShapeType m_shape;
951    
952     //
953     // The number of values needed for the array.
954 jgs 108 // This can be derived from m_shape by multiplying the size of each dimension, but
955     // is stored here for convenience.
956 jgs 82 int m_noValues;
957    
958     };
959    
960 jgs 108 bool operator==(const DataArrayView& left, const DataArrayView& right);
961     bool operator!=(const DataArrayView& left, const DataArrayView& right);
962    
963 jgs 102 /**
964     \brief
965 jgs 108 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 jgs 102 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 jgs 108 /**
993     Inline function definitions.
994     */
995    
996     template <class UnaryFunction>
997 jgs 82 inline
998 jgs 108 void
999     DataArrayView::unaryOp(UnaryFunction operation)
1000 jgs 82 {
1001 jgs 108 unaryOp(m_offset,operation);
1002 jgs 82 }
1003    
1004 jgs 108 template <class UnaryFunction>
1005 jgs 82 inline
1006 jgs 108 void
1007     DataArrayView::unaryOp(ValueType::size_type offset,
1008     UnaryFunction operation)
1009 jgs 82 {
1010 jgs 108 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 jgs 82 }
1016    
1017 jgs 108 template <class BinaryFunction>
1018 jgs 82 inline
1019 jgs 108 void
1020     DataArrayView::binaryOp(const DataArrayView& right,
1021     BinaryFunction operation)
1022 jgs 82 {
1023 jgs 108 binaryOp(m_offset,right,right.getOffset(),operation);
1024 jgs 82 }
1025    
1026 jgs 108 template <class BinaryFunction>
1027 jgs 82 inline
1028     void
1029 jgs 108 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1030     const DataArrayView& right,
1031     ValueType::size_type rightOffset,
1032     BinaryFunction operation)
1033 jgs 82 {
1034 jgs 108 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 jgs 82 }
1043     }
1044    
1045 jgs 108 template <class BinaryFunction>
1046 jgs 82 inline
1047     void
1048 jgs 108 DataArrayView::binaryOp(double right,
1049     BinaryFunction operation)
1050 jgs 82 {
1051 jgs 108 binaryOp(m_offset,right,operation);
1052 jgs 82 }
1053    
1054 jgs 108 template <class BinaryFunction>
1055 jgs 82 inline
1056 jgs 108 void
1057     DataArrayView::binaryOp(ValueType::size_type offset,
1058     double right,
1059     BinaryFunction operation)
1060 jgs 102 {
1061 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1062 jgs 113 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1063 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1064     (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1065 jgs 102 }
1066     }
1067    
1068 jgs 147 template <class BinaryFunction>
1069 jgs 102 inline
1070     double
1071 jgs 147 DataArrayView::reductionOp(BinaryFunction operation,
1072     double initial_value) const
1073 jgs 82 {
1074 jgs 147 return reductionOp(m_offset,operation,initial_value);
1075 jgs 82 }
1076    
1077 jgs 147 template <class BinaryFunction>
1078 jgs 82 inline
1079 jgs 108 double
1080     DataArrayView::reductionOp(ValueType::size_type offset,
1081 jgs 147 BinaryFunction operation,
1082     double initial_value) const
1083 jgs 82 {
1084 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1085     "Error - Couldn't perform reductionOp due to insufficient storage.");
1086 jgs 147 double current_value=initial_value;
1087 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1088 jgs 147 current_value=operation(current_value,(*m_data)[offset+i]);
1089 jgs 82 }
1090 jgs 147 return current_value;
1091 jgs 82 }
1092    
1093     inline
1094 jgs 108 DataArrayView::ValueType::size_type
1095     DataArrayView::relIndex() const
1096 jgs 82 {
1097 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1098     return 0;
1099 jgs 82 }
1100    
1101     inline
1102 jgs 108 DataArrayView::ValueType::size_type
1103     DataArrayView::index() const
1104 jgs 82 {
1105 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1106     return (m_offset);
1107 jgs 82 }
1108    
1109     inline
1110 jgs 108 DataArrayView::ValueType::size_type
1111     DataArrayView::relIndex(ValueType::size_type i) const
1112 jgs 82 {
1113 jgs 108 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 jgs 82 }
1117    
1118     inline
1119     DataArrayView::ValueType::size_type
1120     DataArrayView::index(ValueType::size_type i) const
1121     {
1122 jgs 108 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 jgs 82 }
1126    
1127     inline
1128     DataArrayView::ValueType::size_type
1129     DataArrayView::relIndex(ValueType::size_type i,
1130     ValueType::size_type j) const
1131     {
1132 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1133 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1134 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1135 jgs 82 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 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1144 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1145 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1146     return (m_offset+temp);
1147 jgs 82 }
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 jgs 108 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 jgs 82 }
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 jgs 108 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 jgs 82 }
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 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1181 jgs 82 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 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1183 jgs 82 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 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1194 jgs 82 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 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1196     return (m_offset+temp);
1197 jgs 82 }
1198    
1199     inline
1200     DataArrayView::ValueType::reference
1201 jgs 108 DataArrayView::operator()()
1202 jgs 82 {
1203 jgs 108 return (*m_data)[index()];
1204 jgs 82 }
1205    
1206     inline
1207     DataArrayView::ValueType::const_reference
1208 jgs 108 DataArrayView::operator()() const
1209 jgs 82 {
1210 jgs 108 return (*m_data)[index()];
1211 jgs 82 }
1212    
1213     inline
1214     DataArrayView::ValueType::reference
1215 jgs 108 DataArrayView::operator()(ValueType::size_type i)
1216 jgs 82 {
1217 jgs 108 return (*m_data)[index(i)];
1218 jgs 82 }
1219    
1220     inline
1221     DataArrayView::ValueType::const_reference
1222 jgs 108 DataArrayView::operator()(ValueType::size_type i) const
1223 jgs 82 {
1224 jgs 108 return (*m_data)[index(i)];
1225 jgs 82 }
1226    
1227     inline
1228     DataArrayView::ValueType::reference
1229     DataArrayView::operator()(ValueType::size_type i,
1230     ValueType::size_type j)
1231     {
1232 jgs 108 return (*m_data)[index(i,j)];
1233 jgs 82 }
1234    
1235     inline
1236     DataArrayView::ValueType::const_reference
1237     DataArrayView::operator()(ValueType::size_type i,
1238     ValueType::size_type j) const
1239     {
1240 jgs 108 return (*m_data)[index(i,j)];
1241 jgs 82 }
1242    
1243     inline
1244     DataArrayView::ValueType::reference
1245 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1246     ValueType::size_type j,
1247     ValueType::size_type k)
1248 jgs 82 {
1249 jgs 108 return (*m_data)[index(i,j,k)];
1250 jgs 82 }
1251    
1252     inline
1253     DataArrayView::ValueType::const_reference
1254 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1255     ValueType::size_type j,
1256     ValueType::size_type k) const
1257 jgs 82 {
1258 jgs 108 return (*m_data)[index(i,j,k)];
1259 jgs 82 }
1260    
1261     inline
1262     DataArrayView::ValueType::reference
1263 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1264     ValueType::size_type j,
1265     ValueType::size_type k,
1266     ValueType::size_type m)
1267 jgs 82 {
1268 jgs 108 return (*m_data)[index(i,j,k,m)];
1269 jgs 82 }
1270    
1271     inline
1272     DataArrayView::ValueType::const_reference
1273 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1274     ValueType::size_type j,
1275     ValueType::size_type k,
1276     ValueType::size_type m) const
1277 jgs 82 {
1278 jgs 108 return (*m_data)[index(i,j,k,m)];
1279 jgs 82 }
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