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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 54218 byte(s) 
The MPI branch is hereby closed. All future work should be in trunk.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26