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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 583 - (hide annotations)
Wed Mar 8 08:15:34 2006 UTC (16 years, 11 months ago) by gross
File MIME type: text/plain
File size: 39742 byte(s) 
_eigenvalues_and_eigenvector method added of data object. the algorithm has been tested on floats in python but not on data objects.
1 jgs 102 // $Id$
2 jgs 108
3 jgs 82 /*
4     ******************************************************************************
5     * *
6     * COPYRIGHT ACcESS 2004 - All Rights Reserved *
7     * *
8     * This software is the property of ACcESS. No part of this code *
9     * may be copied in any form or by any means without the expressed written *
10     * consent of ACcESS. Copying, use or modification of this software *
11     * by any unauthorised person is illegal unless that person has a software *
12     * license agreement with ACcESS. *
13     * *
14     ******************************************************************************
15     */
16 jgs 474
17 jgs 82 #if !defined escript_DataArrayView_20040323_H
18     #define escript_DataArrayView_20040323_H
19    
20 jgs 468 #include "EsysAssert.h"
21 jgs 82
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     friend bool operator==(const DataArrayView& left, const DataArrayView& right);
54     friend bool operator!=(const DataArrayView& left, const DataArrayView& right);
55    
56     public:
57    
58 jgs 117 //
59     // Some basic types which define the data values and view shapes.
60     typedef DataVector ValueType;
61     //typedef std::vector<double> ValueType;
62 jgs 108 typedef std::vector<int> ShapeType;
63 jgs 82 typedef std::vector<std::pair<int, int> > RegionType;
64 jgs 108 typedef std::vector<std::pair<int, int> > RegionLoopRangeType;
65 jgs 82
66     /**
67     \brief
68     Default constructor for DataArrayView.
69    
70     Description:
71 jgs 108 Default constructor for DataArrayView.
72     Creates an empty view with no associated data array.
73    
74     This is essentially useless but here for completeness.
75 jgs 82 */
76     DataArrayView();
77    
78     /**
79     \brief
80     Constructor for DataArrayView.
81    
82     Description:
83     Constructor for DataArrayView.
84    
85 jgs 108 \param data - Input -
86     Array holding data to be viewed. This must be large enough
87     to supply sufficient values for the specified shape starting at
88     the given offset.
89     \param viewShape - Input -
90     The shape of the view to create.
91     \param offset - Input -
92     Offset into the data at which view should start.
93 jgs 82 */
94     DataArrayView(ValueType& data,
95     const ShapeType& viewShape,
96     int offset=0);
97    
98     /**
99     \brief
100     Copy constructor for DataArrayView.
101    
102     Description:
103     Copy constructor for DataArrayView.
104    
105 jgs 108 \param other - Input -
106     DataArrayView to copy.
107 jgs 82
108 jgs 108 NOTE: The copy references the same data array.
109 jgs 82 */
110     DataArrayView(const DataArrayView& other);
111    
112     /**
113     \brief
114 jgs 108 Copy from a numarray into the data array viewed by this.
115     This must have same shape as the given value - will throw if not.
116 jgs 82 */
117     void
118     copy(const boost::python::numeric::array& value);
119    
120     /**
121     \brief
122 jgs 108 Copy data from another DataArrayView into the data array viewed by this
123     starting at the default offset in both views.
124     The shapes of both views must be the same - will throw if not.
125     NB: views may or may not reference same underlying data array!
126 jgs 82 */
127     void
128     copy(const DataArrayView& other);
129    
130     /**
131     \brief
132 jgs 108 Copy data from another DataArrayView into this view starting at the
133     given offset in this view and the default offset in the other view.
134     The shapes of both views must be the same - will throw if not.
135     NB: views may or may not reference same underlying data array!
136 jgs 82 */
137     void
138     copy(ValueType::size_type offset,
139     const DataArrayView& other);
140    
141     /**
142     \brief
143 jgs 108 Copy data from another DataArrayView into this view starting at the
144     given offsets in each view.
145     The shapes of both views must be compatible - will throw if not.
146     NB: views may or may not reference same underlying data array!
147 jgs 82
148 jgs 108 \param thisOffset - Input -
149     Offset into this view's data array to copy to.
150     \param other - Input -
151     View to copy data from.
152     \param otherOffset - Input -
153     Offset into the other view's data array to copy from.
154 jgs 82 */
155     void
156     copy(ValueType::size_type thisOffset,
157     const DataArrayView& other,
158     ValueType::size_type otherOffset);
159    
160     /**
161     \brief
162 jgs 108 Copy the given single value over the entire view starting at the given
163     offset in the view's data array.
164 jgs 82
165 jgs 108 \param offset - Input -
166     Offset into this view's data array to copy to.
167     \param value - Input -
168     Value to copy.
169 jgs 82 */
170     void
171     copy(ValueType::size_type offset,
172     ValueType::value_type value);
173    
174     /**
175 jgs 108 \brief
176     Check if view is empty. ie: does not point to any actual data.
177     */
178     bool
179     isEmpty() const;
180    
181     /**
182 jgs 82 \brief
183 jgs 108 Return this view's offset into the viewed data array.
184 jgs 82 */
185     ValueType::size_type
186     getOffset() const;
187    
188     /**
189 jgs 108 \brief
190     Set the offset into the data array at which this view is to start.
191     This could be used to step through the underlying data array by incrementing
192     the offset by noValues successively. Thus this view would provide a moving
193     window on the underlying data with the given shape.
194 jgs 82 */
195 jgs 108 void
196     setOffset(ValueType::size_type offset);
197 jgs 82
198     /**
199     \brief
200 jgs 108 Increment the offset by the number of values in the shape, if possible. Thus
201     moving the view onto the next data point of the given shape in the underlying
202     data array.
203 jgs 82 */
204 jgs 108 void
205     incrOffset();
206 jgs 82
207     /**
208     \brief
209 jgs 108 Check the (given) offset will not result in two few values being available in
210     the underlying data array for this view's shape.
211 jgs 82 */
212 jgs 108 bool
213     checkOffset() const;
214 jgs 117
215 jgs 108 bool
216     checkOffset(ValueType::size_type offset) const;
217 jgs 82
218     /**
219 jgs 108 \brief
220     Return the rank of the shape of this view.
221     */
222     int
223     getRank() const;
224    
225     /**
226 jgs 82 \brief
227 jgs 108 Return the number of values for the shape of this view.
228 jgs 82 */
229     int
230     noValues() const;
231    
232     /**
233     \brief
234 jgs 108 Calculate the number of values for the given shape or region.
235     This is purely a utility method and has no bearing on this view.
236 jgs 82 */
237 jgs 117 static
238     int
239 jgs 108 noValues(const ShapeType& shape);
240 jgs 117
241     static
242     int
243 jgs 108 noValues(const RegionLoopRangeType& region);
244 jgs 82
245     /**
246     \brief
247 jgs 108 Return a reference to the underlying data array.
248 jgs 82 */
249     ValueType&
250     getData() const;
251    
252     /**
253     \brief
254 jgs 108 Return a reference to the data value with the given
255     index in this view. This takes into account the offset.
256 jgs 82 */
257     ValueType::reference
258     getData(ValueType::size_type i) const;
259    
260     /**
261     \brief
262 jgs 108 Return the shape of this view.
263 jgs 82 */
264 jgs 117 const
265     ShapeType&
266 jgs 108 getShape() const;
267 jgs 82
268     /**
269     \brief
270 jgs 108 Return true if the given shape is the same as this view's shape.
271 jgs 82 */
272 jgs 108 bool
273 jgs 117 checkShape(const ShapeType& other) const;
274 jgs 82
275     /**
276     \brief
277 jgs 108 Create a shape error message. Normally used when there is a shape
278 jgs 117 mismatch between this shape and the other shape.
279 jgs 82
280 jgs 108 \param messagePrefix - Input -
281     First part of the error message.
282     \param other - Input -
283     The other shape.
284 jgs 82 */
285 jgs 108 std::string
286     createShapeErrorMessage(const std::string& messagePrefix,
287 jgs 117 const ShapeType& other) const;
288 jgs 82
289     /**
290     \brief
291 jgs 108 Return the viewed data as a formatted string.
292     Not recommended for very large objects!
293 jgs 82
294 jgs 108 \param suffix - Input -
295     Suffix appended to index display.
296 jgs 82 */
297 jgs 108 std::string
298     toString(const std::string& suffix=std::string("")) const;
299 jgs 82
300     /**
301     \brief
302 jgs 108 Return the given shape as a string.
303     This is purely a utility method and has no bearing on this view.
304 jgs 82
305 jgs 108 \param shape - Input.
306 jgs 82 */
307 jgs 117 static
308     std::string
309     shapeToString(const ShapeType& shape);
310 jgs 82
311 jgs 108 /**
312     \brief
313     Return the 1 dimensional index into the data array of the only element
314     in the view, *ignoring the offset*.
315     Assumes a rank 0 view.
316 jgs 82 */
317 jgs 108 ValueType::size_type
318     relIndex() const;
319 jgs 82
320     /**
321 jgs 108 \brief
322     Return the 1 dimensional index into the data array of the element at
323     position i in the view, *ignoring the offset*.
324     Assumes a rank 1 view.
325 jgs 82 */
326 jgs 108 ValueType::size_type
327     relIndex(ValueType::size_type i) const;
328 jgs 82
329     /**
330 jgs 108 \brief
331     Return the 1 dimensional index into the data array of the element at
332     position i,j in the view, *ignoring the offset*.
333     Assumes a rank 2 view.
334 jgs 82 */
335 jgs 108 ValueType::size_type
336     relIndex(ValueType::size_type i,
337     ValueType::size_type j) const;
338 jgs 82
339     /**
340 jgs 108 \brief
341     Return the 1 dimensional index into the data array of the element at
342     position i,j,k in the view, *ignoring the offset*.
343     Assumes a rank 3 view.
344 jgs 82 */
345 jgs 108 ValueType::size_type
346     relIndex(ValueType::size_type i,
347     ValueType::size_type j,
348     ValueType::size_type k) const;
349 jgs 82
350     /**
351     \brief
352 jgs 108 Return the 1 dimensional index into the data array of the element at
353     position i,j,k,m in the view, *ignoring the offset*.
354     Assumes a rank 4 view.
355 jgs 82 */
356     ValueType::size_type
357     relIndex(ValueType::size_type i,
358     ValueType::size_type j,
359     ValueType::size_type k,
360     ValueType::size_type m) const;
361    
362     /**
363     \brief
364 jgs 108 Return the 1 dimensional index into the data array of the only element
365     in the view.
366     Assumes a rank 0 view.
367 jgs 82 */
368     ValueType::size_type
369 jgs 108 index() const;
370 jgs 82
371     /**
372     \brief
373 jgs 108 Return the 1 dimensional index into the data array of the element at
374     position i in the view.
375     Assumes a rank 1 view.
376 jgs 82 */
377     ValueType::size_type
378 jgs 108 index(ValueType::size_type i) 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 in the view.
384     Assumes a rank 2 view.
385 jgs 82 */
386     ValueType::size_type
387     index(ValueType::size_type i,
388 jgs 108 ValueType::size_type j) const;
389 jgs 82
390     /**
391     \brief
392 jgs 108 Return the 1 dimensional index into the data array of the element at
393     position i,j,k in the view.
394     Assumes a rank 3 view.
395 jgs 82 */
396     ValueType::size_type
397     index(ValueType::size_type i,
398     ValueType::size_type j,
399     ValueType::size_type k) const;
400    
401     /**
402     \brief
403 jgs 108 Return the 1 dimensional index into the data array of the element at
404     position i,j,k,m in the view.
405     Assumes a rank 4 view.
406 jgs 82 */
407     ValueType::size_type
408     index(ValueType::size_type i,
409 jgs 108 ValueType::size_type j,
410     ValueType::size_type k,
411     ValueType::size_type m) const;
412 jgs 82
413     /**
414     \brief
415 jgs 108 Return a reference for the only element in the view.
416     Assumes a rank 0 view.
417 jgs 82 */
418 jgs 108 ValueType::reference
419     operator()();
420 jgs 82
421 jgs 108 ValueType::const_reference
422     operator()() const;
423    
424 jgs 82 /**
425     \brief
426 jgs 108 Return a reference to the element at position i in the view.
427     Assumes a rank 1 view.
428 jgs 82 */
429     ValueType::reference
430     operator()(ValueType::size_type i);
431    
432     ValueType::const_reference
433     operator()(ValueType::size_type i) const;
434    
435     /**
436     \brief
437 jgs 108 Return a reference to the element at position i,j in the view.
438     Assumes a rank 2 view.
439 jgs 82 */
440     ValueType::reference
441     operator()(ValueType::size_type i,
442     ValueType::size_type j);
443    
444     ValueType::const_reference
445     operator()(ValueType::size_type i,
446     ValueType::size_type j) const;
447    
448     /**
449     \brief
450 jgs 108 Return a reference to the element at position i,j,k in the view.
451     Assumes a rank 3 view.
452 jgs 82 */
453     ValueType::reference
454     operator()(ValueType::size_type i,
455     ValueType::size_type j,
456     ValueType::size_type k);
457    
458     ValueType::const_reference
459     operator()(ValueType::size_type i,
460     ValueType::size_type j,
461     ValueType::size_type k) const;
462    
463     /**
464     \brief
465 jgs 108 Return a reference to the element at position i,j,k,m in the view.
466     Assumes a rank 4 view.
467 jgs 82 */
468     ValueType::reference
469     operator()(ValueType::size_type i,
470     ValueType::size_type j,
471     ValueType::size_type k,
472     ValueType::size_type m);
473    
474     ValueType::const_reference
475     operator()(ValueType::size_type i,
476     ValueType::size_type j,
477     ValueType::size_type k,
478     ValueType::size_type m) const;
479    
480     /**
481 jgs 108 \brief
482     Determine the shape of the specified slice region.
483     This is purely a utility method and has no bearing on this view.
484    
485     \param region - Input -
486     Slice region.
487 jgs 82 */
488 jgs 117 static
489     ShapeType
490 jgs 108 getResultSliceShape(const RegionType& region);
491 jgs 82
492 jgs 108 /**
493     \brief
494     Determine the region specified by the given python slice object.
495 jgs 82
496 jgs 108 \param key - Input -
497     python slice object specifying region to be returned.
498    
499     The slice object is a tuple of n python slice specifiers, where
500     n <= the rank of this Data object. Each slice specifier specifies the
501     range of indexes to be sliced from the corresponding dimension. The
502     first specifier corresponds to the first dimension, the second to the
503     second and so on. Where n < the rank, the remaining dimensions are
504     sliced across the full range of their indicies.
505    
506     Each slice specifier is of the form "a:b", which specifies a slice
507     from index a, up to but not including index b. Where index a is ommitted
508     a is assumed to be 0. Where index b is ommitted, b is assumed to be the
509     length of this dimension. Where both are ommitted (eg: ":") the slice is
510     assumed to encompass that entire dimension.
511    
512     Where one of the slice specifiers is a single integer, eg: [1], we
513     want to generate a rank-1 dimension object, as opposed to eg: [1,2]
514     which implies we want to take a rank dimensional object with one
515     dimension of size 1.
516    
517     The return value is a vector of pairs with length equal to the rank of
518     this object. Each pair corresponds to the range of indicies from the
519     corresponding dimension to be sliced from, as specified in the input
520     slice object.
521    
522     Examples:
523    
524     For a rank 1 object of shape(5):
525    
526     getSliceRegion(:) => < <0,5> >
527     getSliceRegion(2:3) => < <2,3> >
528     getSliceRegion(:3) => < <0,3> >
529     getSliceRegion(2:) => < <2,5> >
530    
531     For a rank 2 object of shape(4,5):
532    
533     getSliceRegion(2:3) => < <2,3> <0,5> >
534     getSliceRegion(2) => < <2,3> <0,5> >
535     NB: but return object requested will have rank 1, shape(5), with
536     values taken from index 2 of this object's first dimension.
537    
538     For a rank 3 object of shape (2,4,6):
539    
540     getSliceRegion(0:2,0:4,0:6) => < <0,2> <0,4> <0,6> >
541     getSliceRegion(:,:,:) => < <0,2> <0,4> <0,6> >
542     getSliceRegion(0:1) => < <0,1> <0,4> <0,6> >
543     getSliceRegion(:1,0:2) => < <0,1> <0,2> <0,6> >
544    
545     */
546 jgs 117 RegionType
547 jgs 108 getSliceRegion(const boost::python::object& key) const;
548    
549 jgs 82 /**
550     \brief
551 jgs 108 Copy a data slice specified by the given region from the given view
552     into this view, using the default offsets in both views.
553    
554     \param other - Input -
555     View to copy data from.
556     \param region - Input -
557     Region in other view to copy data from.
558 jgs 82 */
559     void
560 jgs 108 copySlice(const DataArrayView& other,
561     const RegionLoopRangeType& region);
562 jgs 82
563     /**
564     \brief
565 jgs 108 Copy a data slice specified by the given region and offset from the
566     given view into this view at the given offset.
567    
568     \param thisOffset - Input -
569     Copy the slice to this offset in this view.
570     \param other - Input -
571     View to copy data from.
572     \param otherOffset - Input -
573     Copy the slice from this offset in the given view.
574     \param region - Input -
575     Region in other view to copy data from.
576 jgs 82 */
577     void
578 jgs 108 copySlice(ValueType::size_type thisOffset,
579     const DataArrayView& other,
580     ValueType::size_type otherOffset,
581     const RegionLoopRangeType& region);
582 jgs 82
583     /**
584     \brief
585 jgs 108 Copy data into a slice specified by the given region in this view from
586     the given view, using the default offsets in both views.
587 jgs 102
588 jgs 108 \param other - Input -
589     View to copy data from.
590     \param region - Input -
591     Region in this view to copy data to.
592 jgs 102 */
593 jgs 108 void
594     copySliceFrom(const DataArrayView& other,
595     const RegionLoopRangeType& region);
596 jgs 102
597     /**
598     \brief
599 jgs 108 Copy data into a slice specified by the given region and offset in
600     this view from the given view at the given offset.
601 jgs 102
602 jgs 108 \param thisOffset - Input -
603     Copy the slice to this offset in this view.
604     \param other - Input -
605     View to copy data from.
606     \param otherOffset - Input -
607     Copy the slice from this offset in the given view.
608     \param region - Input -
609     Region in this view to copy data to.
610 jgs 102 */
611 jgs 108 void
612     copySliceFrom(ValueType::size_type thisOffset,
613     const DataArrayView& other,
614     ValueType::size_type otherOffset,
615     const RegionLoopRangeType& region);
616 jgs 102
617     /**
618     \brief
619 jgs 108 Perform the unary operation on the data point specified by the view's
620     default offset. Applies the specified operation to each value in the data
621     point.
622    
623     Called by escript::unaryOp.
624    
625     \param operation - Input -
626     Operation to apply. Operation must be a pointer to a function.
627 jgs 82 */
628     template <class UnaryFunction>
629 jgs 108 void
630     unaryOp(UnaryFunction operation);
631 jgs 82
632     /**
633     \brief
634 jgs 108 Perform the unary operation on the data point specified by the given
635     offset. Applies the specified operation to each value in the data
636     point. Operation must be a pointer to a function.
637    
638     Called by escript::unaryOp.
639    
640     \param offset - Input -
641     Apply the operation to data starting at this offset in this view.
642     \param operation - Input -
643     Operation to apply. Must be a pointer to a function.
644 jgs 82 */
645     template <class UnaryFunction>
646 jgs 108 void
647     unaryOp(ValueType::size_type offset,
648     UnaryFunction operation);
649 jgs 82
650     /**
651     \brief
652 jgs 108 Perform the binary operation on the data points specified by the default
653     offsets in this view and in view "right". Applies the specified operation
654     to corresponding values in both data points. Operation must be a pointer
655     to a function.
656    
657     Called by escript::binaryOp.
658    
659     \param right - Input -
660     View to act as RHS in given binary operation.
661     \param operation - Input -
662     Operation to apply. Must be a pointer to a function.
663 jgs 82 */
664     template <class BinaryFunction>
665     void
666 jgs 108 binaryOp(const DataArrayView& right,
667 jgs 82 BinaryFunction operation);
668    
669     /**
670     \brief
671 jgs 108 Perform the binary operation on the data points specified by the given
672     offsets in this view and in view "right". Applies the specified operation
673     to corresponding values in both data points. Operation must be a pointer
674     to a function.
675    
676     Called by escript::binaryOp.
677    
678     \param leftOffset - Input -
679     Apply the operation to data starting at this offset in this view.
680     \param right - Input -
681     View to act as RHS in given binary operation.
682     \param rightOffset - Input -
683     Apply the operation to data starting at this offset in right.
684     \param operation - Input -
685     Operation to apply. Must be a pointer to a function.
686 jgs 82 */
687     template <class BinaryFunction>
688     void
689 jgs 108 binaryOp(ValueType::size_type leftOffset,
690     const DataArrayView& right,
691     ValueType::size_type rightOffset,
692 jgs 82 BinaryFunction operation);
693    
694     /**
695     \brief
696 jgs 108 Perform the binary operation on the data point specified by the default
697     offset in this view using the scalar value "right". Applies the specified
698     operation to values in the data point. Operation must be a pointer to
699     a function.
700    
701     Called by escript::binaryOp.
702    
703     \param right - Input -
704     Value to act as RHS in given binary operation.
705     \param operation - Input -
706     Operation to apply. Must be a pointer to a function.
707 jgs 82 */
708     template <class BinaryFunction>
709     void
710 jgs 108 binaryOp(double right,
711 jgs 82 BinaryFunction operation);
712    
713     /**
714     \brief
715 jgs 108 Perform the binary operation on the data point specified by the given
716     offset in this view using the scalar value "right". Applies the specified
717     operation to values in the data point. Operation must be a pointer
718     to a function.
719    
720     Called by escript::binaryOp.
721    
722     \param offset - Input -
723     Apply the operation to data starting at this offset in this view.
724     \param right - Input -
725     Value to act as RHS in given binary operation.
726     \param operation - Input -
727     Operation to apply. Must be a pointer to a function.
728 jgs 82 */
729     template <class BinaryFunction>
730     void
731 jgs 108 binaryOp(ValueType::size_type offset,
732     double right,
733 jgs 82 BinaryFunction operation);
734    
735     /**
736     \brief
737 jgs 108 Perform the given data point reduction operation on the data point
738     specified by the default offset into the view. Reduces all elements of
739     the data point using the given operation, returning the result as a
740     scalar. Operation must be a pointer to a function.
741    
742     Called by escript::algorithm.
743    
744     \param operation - Input -
745     Operation to apply. Must be a pointer to a function.
746     */
747 jgs 147 template <class BinaryFunction>
748 jgs 108 double
749 jgs 147 reductionOp(BinaryFunction operation,
750     double initial_value) const;
751 jgs 108
752     /**
753     \brief
754 jgs 113 Perform the given data point reduction operation on the data point
755 jgs 108 specified by the given offset into the view. Reduces all elements of
756     the data point using the given operation, returning the result as a
757     scalar. Operation must be a pointer to a function.
758    
759     Called by escript::algorithm.
760    
761     \param offset - Input -
762     Apply the operation to data starting at this offset in this view.
763     \param operation - Input -
764     Operation to apply. Must be a pointer to a function.
765     */
766 jgs 147 template <class BinaryFunction>
767 jgs 108 double
768     reductionOp(ValueType::size_type offset,
769 jgs 147 BinaryFunction operation,
770     double initial_value) const;
771 jgs 108
772     /**
773     \brief
774 gross 576 solves a local eigenvalue problem
775    
776 gross 580 \param in - Input - matrix
777     \param inOffset - Input - offset into in
778 gross 576 \param ev - Output - The eigenvalues
779 gross 580 \param inOffset - Input - offset into ev
780 gross 576 */
781 gross 580 static
782 gross 576 inline
783     void
784 gross 580 DataArrayView::eigenvalues(DataArrayView& in,
785     ValueType::size_type inOffset,
786     DataArrayView& ev,
787     ValueType::size_type evOffset)
788     {
789     double in00,in10,in20,in01,in11,in21,in02,in12,in22;
790     double ev0,ev1,ev2;
791     int s=in.getShape()[0];
792     if (s==1) {
793     in00=(*(in.m_data))[inOffset+in.index(0,0)];
794     eigenvalues1(in00,&ev0);
795     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
796 gross 576
797 gross 580 } else if (s==2) {
798     in00=(*(in.m_data))[inOffset+in.index(0,0)];
799     in10=(*(in.m_data))[inOffset+in.index(1,0)];
800     in01=(*(in.m_data))[inOffset+in.index(0,1)];
801     in11=(*(in.m_data))[inOffset+in.index(1,1)];
802     eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
803     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
804     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
805 gross 576
806 gross 580 } else if (s==3) {
807     in00=(*(in.m_data))[inOffset+in.index(0,0)];
808     in10=(*(in.m_data))[inOffset+in.index(1,0)];
809     in20=(*(in.m_data))[inOffset+in.index(2,0)];
810     in01=(*(in.m_data))[inOffset+in.index(0,1)];
811     in11=(*(in.m_data))[inOffset+in.index(1,1)];
812     in21=(*(in.m_data))[inOffset+in.index(2,1)];
813     in02=(*(in.m_data))[inOffset+in.index(0,2)];
814     in12=(*(in.m_data))[inOffset+in.index(1,2)];
815     in22=(*(in.m_data))[inOffset+in.index(2,2)];
816     eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
817     &ev0,&ev1,&ev2);
818     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
819     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
820     (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
821    
822     }
823     }
824    
825 gross 576 /**
826     \brief
827 gross 580 solves a local eigenvalue problem
828    
829     \param in - Input - matrix
830     \param inOffset - Input - offset into in
831     \param ev - Output - The eigenvalues
832     \param evOffset - Input - offset into ev
833     \param V - Output - The eigenvectors
834     \param VOffset - Input - offset into V
835     \param tol - Input - eigenvalues with relative difference tol are treated as equal
836     */
837     static
838     inline
839     void
840     DataArrayView::eigenvalues_and_eigenvectors(DataArrayView& in,
841     ValueType::size_type inOffset,
842     DataArrayView& ev,
843     ValueType::size_type evOffset,
844     DataArrayView& V,
845     ValueType::size_type VOffset,
846     const double tol=1.e-13)
847     {
848 gross 583 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
849     double V00,V10,V20,V01,V11,V21,V02,V12,V22;
850 gross 580 double ev0,ev1,ev2;
851     int s=in.getShape()[0];
852     if (s==1) {
853 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
854     eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
855     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
856     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
857 gross 580 } else if (s==2) {
858 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
859     in10=(*(in.m_data))[inOffset+in.index(1,0)];
860     in01=(*(in.m_data))[inOffset+in.index(0,1)];
861     in11=(*(in.m_data))[inOffset+in.index(1,1)];
862     eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
863     &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
864     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
865     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
866     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
867     (*(V.m_data))[inOffset+V.index(1,0)]=V10;
868     (*(V.m_data))[inOffset+V.index(0,1)]=V01;
869     (*(V.m_data))[inOffset+V.index(1,1)]=V11;
870 gross 580 } else if (s==3) {
871 gross 583 in00=(*(in.m_data))[inOffset+in.index(0,0)];
872     in10=(*(in.m_data))[inOffset+in.index(1,0)];
873     in20=(*(in.m_data))[inOffset+in.index(2,0)];
874     in01=(*(in.m_data))[inOffset+in.index(0,1)];
875     in11=(*(in.m_data))[inOffset+in.index(1,1)];
876     in21=(*(in.m_data))[inOffset+in.index(2,1)];
877     in02=(*(in.m_data))[inOffset+in.index(0,2)];
878     in12=(*(in.m_data))[inOffset+in.index(1,2)];
879     in22=(*(in.m_data))[inOffset+in.index(2,2)];
880     eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
881     &ev0,&ev1,&ev2,
882     &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
883     (*(ev.m_data))[evOffset+ev.index(0)]=ev0;
884     (*(ev.m_data))[evOffset+ev.index(1)]=ev1;
885     (*(ev.m_data))[evOffset+ev.index(2)]=ev2;
886     (*(V.m_data))[inOffset+V.index(0,0)]=V00;
887     (*(V.m_data))[inOffset+V.index(1,0)]=V10;
888     (*(V.m_data))[inOffset+V.index(2,0)]=V20;
889     (*(V.m_data))[inOffset+V.index(0,1)]=V01;
890     (*(V.m_data))[inOffset+V.index(1,1)]=V11;
891     (*(V.m_data))[inOffset+V.index(2,1)]=V21;
892     (*(V.m_data))[inOffset+V.index(0,2)]=V02;
893     (*(V.m_data))[inOffset+V.index(1,2)]=V12;
894     (*(V.m_data))[inOffset+V.index(2,2)]=V22;
895 gross 580
896     }
897     }
898     /**
899     \brief
900 jgs 113 Perform a matrix multiply of the given views.
901     This is purely a utility method and has no bearing on this view.
902 jgs 82
903 jgs 126 NB: Only multiplies together the two given views, ie: two data-points,
904     would need to call this over all data-points to multiply the entire
905     Data objects involved.
906    
907 jgs 82 \param left - Input - The left hand side.
908     \param right - Input - The right hand side.
909     \param result - Output - The result of the operation.
910     */
911 jgs 117 static
912     void
913 jgs 82 matMult(const DataArrayView& left,
914     const DataArrayView& right,
915     DataArrayView& result);
916    
917 jgs 108 /**
918     \brief
919     Determine the shape of the result array for a matrix multiplication
920     of the given views.
921     This is purely a utility method and has no bearing on this view.
922     */
923 jgs 117 static
924     ShapeType
925 jgs 108 determineResultShape(const DataArrayView& left,
926     const DataArrayView& right);
927    
928 jgs 82 protected:
929    
930     private:
931    
932     //
933 jgs 108 // The maximum rank allowed for the shape of any view.
934 jgs 82 static const int m_maxRank=4;
935    
936     //
937 jgs 108 // The data values for the view.
938     // NOTE: This points to data external to the view.
939 jgs 117 // This is just a pointer to an array of ValueType.
940 jgs 82 ValueType* m_data;
941    
942     //
943     // The offset into the data array used by different views.
944 jgs 108 // This is simply an integer specifying a position in the data array
945     // pointed to by m_data.
946 jgs 82 ValueType::size_type m_offset;
947    
948     //
949     // The shape of the data.
950 jgs 108 // This is simply an STL vector specifying the lengths of each dimension
951     // of the shape as ints.
952 jgs 82 ShapeType m_shape;
953    
954     //
955     // The number of values needed for the array.
956 jgs 108 // This can be derived from m_shape by multiplying the size of each dimension, but
957     // is stored here for convenience.
958 jgs 82 int m_noValues;
959    
960     };
961    
962 jgs 108 bool operator==(const DataArrayView& left, const DataArrayView& right);
963     bool operator!=(const DataArrayView& left, const DataArrayView& right);
964    
965 jgs 102 /**
966     \brief
967 jgs 108 Modify region to copy from in order to
968     deal with the case where one range in the region contains identical indexes,
969     eg: <<1,1><0,3><0,3>>
970     This situation implies we want to copy from an object with rank greater than that of this
971     object. eg: we want to copy the values from a two dimensional slice out of a three
972     dimensional object into a two dimensional object.
973     We do this by taking a slice from the other object where one dimension of
974     the slice region is of size 1. So in the above example, we modify the above
975     region like so: <<1,2><0,3><0,3>> and take this slice.
976     */
977     DataArrayView::RegionLoopRangeType
978     getSliceRegionLoopRange(const DataArrayView::RegionType& region);
979    
980     /**
981     \brief
982 jgs 102 Calculate the slice range from the given python key object
983     Used by DataArrayView::getSliceRegion.
984     Returns the python slice object key as a pair of ints where the first
985     member is the start and the second member is the end. the presence of a possible
986     step attribute with value different from one will throw an exception
987    
988     /param key - Input - key object specifying slice range.
989     */
990     std::pair<int,int>
991     getSliceRange(const boost::python::object& key,
992     const int shape);
993    
994 jgs 108 /**
995     Inline function definitions.
996     */
997    
998     template <class UnaryFunction>
999 jgs 82 inline
1000 jgs 108 void
1001     DataArrayView::unaryOp(UnaryFunction operation)
1002 jgs 82 {
1003 jgs 108 unaryOp(m_offset,operation);
1004 jgs 82 }
1005    
1006 jgs 108 template <class UnaryFunction>
1007 jgs 82 inline
1008 jgs 108 void
1009     DataArrayView::unaryOp(ValueType::size_type offset,
1010     UnaryFunction operation)
1011 jgs 82 {
1012 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1013     "Error - Couldn't perform unaryOp due to insufficient storage.");
1014     for (ValueType::size_type i=0;i<noValues();i++) {
1015     (*m_data)[offset+i]=operation((*m_data)[offset+i]);
1016     }
1017 jgs 82 }
1018    
1019 jgs 108 template <class BinaryFunction>
1020 jgs 82 inline
1021 jgs 108 void
1022     DataArrayView::binaryOp(const DataArrayView& right,
1023     BinaryFunction operation)
1024 jgs 82 {
1025 jgs 108 binaryOp(m_offset,right,right.getOffset(),operation);
1026 jgs 82 }
1027    
1028 jgs 108 template <class BinaryFunction>
1029 jgs 82 inline
1030     void
1031 jgs 108 DataArrayView::binaryOp(ValueType::size_type leftOffset,
1032     const DataArrayView& right,
1033     ValueType::size_type rightOffset,
1034     BinaryFunction operation)
1035 jgs 82 {
1036 jgs 108 EsysAssert(getShape()==right.getShape(),
1037     "Error - Couldn't perform binaryOp due to shape mismatch,");
1038     EsysAssert((!isEmpty()&&checkOffset(leftOffset)),
1039     "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1040     EsysAssert((!right.isEmpty()&&right.checkOffset(rightOffset)),
1041     "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
1042     for (ValueType::size_type i=0;i<noValues();i++) {
1043     (*m_data)[leftOffset+i]=operation((*m_data)[leftOffset+i],(*right.m_data)[rightOffset+i]);
1044 jgs 82 }
1045     }
1046    
1047 jgs 108 template <class BinaryFunction>
1048 jgs 82 inline
1049     void
1050 jgs 108 DataArrayView::binaryOp(double right,
1051     BinaryFunction operation)
1052 jgs 82 {
1053 jgs 108 binaryOp(m_offset,right,operation);
1054 jgs 82 }
1055    
1056 jgs 108 template <class BinaryFunction>
1057 jgs 82 inline
1058 jgs 108 void
1059     DataArrayView::binaryOp(ValueType::size_type offset,
1060     double right,
1061     BinaryFunction operation)
1062 jgs 102 {
1063 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1064 jgs 113 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
1065 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1066     (*m_data)[offset+i]=operation((*m_data)[offset+i],right);
1067 jgs 102 }
1068     }
1069    
1070 jgs 147 template <class BinaryFunction>
1071 jgs 102 inline
1072     double
1073 jgs 147 DataArrayView::reductionOp(BinaryFunction operation,
1074     double initial_value) const
1075 jgs 82 {
1076 jgs 147 return reductionOp(m_offset,operation,initial_value);
1077 jgs 82 }
1078    
1079 jgs 147 template <class BinaryFunction>
1080 jgs 82 inline
1081 jgs 108 double
1082     DataArrayView::reductionOp(ValueType::size_type offset,
1083 jgs 147 BinaryFunction operation,
1084     double initial_value) const
1085 jgs 82 {
1086 jgs 108 EsysAssert((!isEmpty()&&checkOffset(offset)),
1087     "Error - Couldn't perform reductionOp due to insufficient storage.");
1088 jgs 147 double current_value=initial_value;
1089 jgs 108 for (ValueType::size_type i=0;i<noValues();i++) {
1090 jgs 147 current_value=operation(current_value,(*m_data)[offset+i]);
1091 jgs 82 }
1092 jgs 147 return current_value;
1093 jgs 82 }
1094    
1095     inline
1096 jgs 108 DataArrayView::ValueType::size_type
1097     DataArrayView::relIndex() const
1098 jgs 82 {
1099 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1100     return 0;
1101 jgs 82 }
1102    
1103     inline
1104 jgs 108 DataArrayView::ValueType::size_type
1105     DataArrayView::index() const
1106 jgs 82 {
1107 jgs 108 EsysAssert((getRank()==0),"Incorrect number of indices for the rank of this object.");
1108     return (m_offset);
1109 jgs 82 }
1110    
1111     inline
1112 jgs 108 DataArrayView::ValueType::size_type
1113     DataArrayView::relIndex(ValueType::size_type i) const
1114 jgs 82 {
1115 jgs 108 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1116     EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1117     return i;
1118 jgs 82 }
1119    
1120     inline
1121     DataArrayView::ValueType::size_type
1122     DataArrayView::index(ValueType::size_type i) const
1123     {
1124 jgs 108 EsysAssert((getRank()==1),"Incorrect number of indices for the rank of this object.");
1125     EsysAssert((i < noValues(m_shape)), "Error - Invalid index.");
1126     return (m_offset+i);
1127 jgs 82 }
1128    
1129     inline
1130     DataArrayView::ValueType::size_type
1131     DataArrayView::relIndex(ValueType::size_type i,
1132     ValueType::size_type j) const
1133     {
1134 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1135 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1136 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1137 jgs 82 return temp;
1138     }
1139    
1140     inline
1141     DataArrayView::ValueType::size_type
1142     DataArrayView::index(ValueType::size_type i,
1143     ValueType::size_type j) const
1144     {
1145 jgs 108 EsysAssert((getRank()==2),"Incorrect number of indices for the rank of this object.");
1146 jgs 82 ValueType::size_type temp=i+j*m_shape[0];
1147 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1148     return (m_offset+temp);
1149 jgs 82 }
1150    
1151     inline
1152     DataArrayView::ValueType::size_type
1153     DataArrayView::relIndex(ValueType::size_type i,
1154     ValueType::size_type j,
1155     ValueType::size_type k) const
1156     {
1157 jgs 108 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1158     ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1159     EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1160     return temp;
1161 jgs 82 }
1162    
1163     inline
1164     DataArrayView::ValueType::size_type
1165     DataArrayView::index(ValueType::size_type i,
1166     ValueType::size_type j,
1167     ValueType::size_type k) const
1168     {
1169 jgs 108 EsysAssert((getRank()==3),"Incorrect number of indices for the rank of this object.");
1170     ValueType::size_type temp=i+j*m_shape[0]+k*m_shape[1]*m_shape[0];
1171     EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1172     return (m_offset+temp);
1173 jgs 82 }
1174    
1175     inline
1176     DataArrayView::ValueType::size_type
1177     DataArrayView::relIndex(ValueType::size_type i,
1178     ValueType::size_type j,
1179     ValueType::size_type k,
1180     ValueType::size_type m) const
1181     {
1182 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1183 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];
1184 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1185 jgs 82 return temp;
1186     }
1187    
1188     inline
1189     DataArrayView::ValueType::size_type
1190     DataArrayView::index(ValueType::size_type i,
1191     ValueType::size_type j,
1192     ValueType::size_type k,
1193     ValueType::size_type m) const
1194     {
1195 jgs 108 EsysAssert((getRank()==4),"Incorrect number of indices for the rank of this object.");
1196 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];
1197 jgs 108 EsysAssert((temp < noValues(m_shape)), "Error - Invalid index.");
1198     return (m_offset+temp);
1199 jgs 82 }
1200    
1201     inline
1202     DataArrayView::ValueType::reference
1203 jgs 108 DataArrayView::operator()()
1204 jgs 82 {
1205 jgs 108 return (*m_data)[index()];
1206 jgs 82 }
1207    
1208     inline
1209     DataArrayView::ValueType::const_reference
1210 jgs 108 DataArrayView::operator()() const
1211 jgs 82 {
1212 jgs 108 return (*m_data)[index()];
1213 jgs 82 }
1214    
1215     inline
1216     DataArrayView::ValueType::reference
1217 jgs 108 DataArrayView::operator()(ValueType::size_type i)
1218 jgs 82 {
1219 jgs 108 return (*m_data)[index(i)];
1220 jgs 82 }
1221    
1222     inline
1223     DataArrayView::ValueType::const_reference
1224 jgs 108 DataArrayView::operator()(ValueType::size_type i) const
1225 jgs 82 {
1226 jgs 108 return (*m_data)[index(i)];
1227 jgs 82 }
1228    
1229     inline
1230     DataArrayView::ValueType::reference
1231     DataArrayView::operator()(ValueType::size_type i,
1232     ValueType::size_type j)
1233     {
1234 jgs 108 return (*m_data)[index(i,j)];
1235 jgs 82 }
1236    
1237     inline
1238     DataArrayView::ValueType::const_reference
1239     DataArrayView::operator()(ValueType::size_type i,
1240     ValueType::size_type j) const
1241     {
1242 jgs 108 return (*m_data)[index(i,j)];
1243 jgs 82 }
1244    
1245     inline
1246     DataArrayView::ValueType::reference
1247 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1248     ValueType::size_type j,
1249     ValueType::size_type k)
1250 jgs 82 {
1251 jgs 108 return (*m_data)[index(i,j,k)];
1252 jgs 82 }
1253    
1254     inline
1255     DataArrayView::ValueType::const_reference
1256 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1257     ValueType::size_type j,
1258     ValueType::size_type k) const
1259 jgs 82 {
1260 jgs 108 return (*m_data)[index(i,j,k)];
1261 jgs 82 }
1262    
1263     inline
1264     DataArrayView::ValueType::reference
1265 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1266     ValueType::size_type j,
1267     ValueType::size_type k,
1268     ValueType::size_type m)
1269 jgs 82 {
1270 jgs 108 return (*m_data)[index(i,j,k,m)];
1271 jgs 82 }
1272    
1273     inline
1274     DataArrayView::ValueType::const_reference
1275 jgs 108 DataArrayView::operator()(ValueType::size_type i,
1276     ValueType::size_type j,
1277     ValueType::size_type k,
1278     ValueType::size_type m) const
1279 jgs 82 {
1280 jgs 108 return (*m_data)[index(i,j,k,m)];
1281 jgs 82 }
1282    
1283     } // end of namespace
1284    
1285     #endif

Properties

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