/[escript]/branches/windows_from_1456_trunk_1522_merged_in/escript/src/DataArrayView.h
ViewVC logotype

Annotation of /branches/windows_from_1456_trunk_1522_merged_in/escript/src/DataArrayView.h

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26