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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 950 - (hide annotations)
Tue Feb 6 07:01:11 2007 UTC (12 years, 8 months ago) by gross
File MIME type: text/plain
File size: 54501 byte(s)
escript data objects can now be saved to netCDF files, see http://www.unidata.ucar.edu/software/netcdf/.
Currently only constant data are implemented with expanded and tagged data to follow.
There are two new functions to dump a data object

   s=Data(...)
   s.dump(<filename>)

and to recover it

   s=load(<filename>, domain)

Notice that the function space of s is recovered but domain is still need. 

dump and load will replace archive and extract.

The installation needs now the netCDF installed. 


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26