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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 800 - (hide annotations)
Tue Aug 8 11:23:18 2006 UTC (13 years, 4 months ago) by gross
File MIME type: text/plain
File size: 52849 byte(s)
new function _swap. Python wrapper + testing is still missing.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26