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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1833 - (hide annotations)
Thu Oct 2 09:47:24 2008 UTC (11 years, 6 months ago) by phornby
File MIME type: text/plain
File size: 34389 byte(s)
One cannot give a static function external linkage
(eigenvalues_and_eigenvectors and swapaxes).
I'm not even sure the external linkage is needed for inline functions.
1 jfenwick 1724
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
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 jfenwick 1724
14    
15     #ifndef escript_DataMaths_20080822_H
16 jfenwick 1726 #define escript_DataMaths_20080822_H
17 jfenwick 1724 #include "DataAbstract.h"
18 jfenwick 1747 #include "DataException.h"
19     #include "LocalOps.h"
20 jfenwick 1724
21     namespace escript
22     {
23     namespace DataMaths
24     {
25    
26 jfenwick 1779 /**
27     \namespace escript::DataMaths
28     \brief Contains maths operations performed on data vectors.
29 jfenwick 1724
30 jfenwick 1779 In order to properly identify the datapoints, in most cases, the vector, shape and offset of the point must all be supplied.
31     Note that vector in this context refers to a data vector storing datapoints not a mathematical vector. (However, datapoints within the data vector could represent scalars, vectors, matricies, ...).
32     */
33 jfenwick 1724
34     // /**
35     // \brief
36     // Perform the unary operation on the data point specified by the view's
37     // default offset. Applies the specified operation to each value in the data
38     // point.
39     //
40     // Called by escript::unaryOp.
41     //
42     // \param operation - Input -
43     // Operation to apply. Operation must be a pointer to a function.
44     // */
45     // template <class UnaryFunction>
46     // void
47     // unaryOp(DataAbstract& data, UnaryFunction operation);
48    
49     // I'm going to try not to have the above function. It relies on the value of the offset already
50     // being known. I don't want that, offsets need to be explicit.
51    
52    
53    
54     /**
55     \brief
56     Perform the unary operation on the data point specified by the given
57     offset. Applies the specified operation to each value in the data
58     point. Operation must be a pointer to a function.
59    
60     Called by escript::unaryOp.
61    
62 jfenwick 1779 \param data - vector containing the datapoint
63     \param shape - shape of the point
64     \param offset - offset of the point within data
65 jfenwick 1724 \param operation - Input -
66     Operation to apply. Must be a pointer to a function.
67     */
68     template <class UnaryFunction>
69     void
70     unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
71     DataTypes::ValueType::size_type offset,
72     UnaryFunction operation);
73    
74     // /**
75     // \brief
76     // Perform the binary operation on the data points specified by the default
77     // offsets in this view and in view "right". Applies the specified operation
78     // to corresponding values in both data points. Operation must be a pointer
79     // to a function.
80     //
81     // Called by escript::binaryOp.
82     //
83     // \param right - Input -
84     // View to act as RHS in given binary operation.
85     // \param operation - Input -
86     // Operation to apply. Must be a pointer to a function.
87     // */
88     // template <class BinaryFunction>
89     // void
90     // binaryOp(DataAbstract& left, const DataAbstract& right,
91     // BinaryFunction operation);
92    
93     // trying to avoid having this one as well. Again, implicit offset
94    
95     /**
96     \brief
97     Perform the binary operation on the data points specified by the given
98 jfenwick 1779 offsets in the "left" and "right" vectors. Applies the specified operation
99 jfenwick 1724 to corresponding values in both data points. Operation must be a pointer
100     to a function.
101    
102     Called by escript::binaryOp.
103 jfenwick 1779 \param left,right - vectors containing the datapoints
104     \param leftShape,rightShape - shapes of datapoints in the vectors
105     \param leftOffset,rightOffset - beginnings of datapoints in the vectors
106 jfenwick 1724 \param operation - Input -
107     Operation to apply. Must be a pointer to a function.
108     */
109     template <class BinaryFunction>
110     void
111 jfenwick 1779 binaryOp(DataTypes::ValueType& left,
112     const DataTypes::ShapeType& leftShape,
113     DataTypes::ValueType::size_type leftOffset,
114     const DataTypes::ValueType& right,
115     const DataTypes::ShapeType& rightShape,
116 jfenwick 1724 DataTypes::ValueType::size_type rightOffset,
117     BinaryFunction operation);
118    
119     // /**
120     // \brief
121     // Perform the binary operation on the data point specified by the default
122     // offset in this view using the scalar value "right". Applies the specified
123     // operation to values in the data point. Operation must be a pointer to
124     // a function.
125     //
126     // Called by escript::binaryOp.
127     //
128     // \param right - Input -
129     // Value to act as RHS in given binary operation.
130     // \param operation - Input -
131     // Operation to apply. Must be a pointer to a function.
132     // */
133     // template <class BinaryFunction>
134     // void
135     // binaryOp(DataAbstract& left, double right,
136     // BinaryFunction operation);
137    
138     // Implicit offset again
139    
140    
141     /**
142     \brief
143     Perform the binary operation on the data point specified by the given
144 jfenwick 1779 offset in the vector using the scalar value "right". Applies the specified
145 jfenwick 1724 operation to values in the data point. Operation must be a pointer
146     to a function.
147    
148     Called by escript::binaryOp.
149    
150 jfenwick 1779 \param left - vector containing the datapoints
151     \param shape - shape of datapoint in the vector
152     \param offset - beginning of datapoint in the vector
153     \param right - scalar value for the right hand side of the operation
154 jfenwick 1724 \param operation - Input -
155     Operation to apply. Must be a pointer to a function.
156     */
157     template <class BinaryFunction>
158     void
159 jfenwick 1779 binaryOp(DataTypes::ValueType& left,
160     const DataTypes::ShapeType& shape,
161 jfenwick 1724 DataTypes::ValueType::size_type offset,
162     double right,
163     BinaryFunction operation);
164    
165     // /**
166     // \brief
167     // Perform the given data point reduction operation on the data point
168     // specified by the default offset into the view. Reduces all elements of
169     // the data point using the given operation, returning the result as a
170     // scalar. Operation must be a pointer to a function.
171     //
172     // Called by escript::algorithm.
173     //
174     // \param operation - Input -
175     // Operation to apply. Must be a pointer to a function.
176     // */
177     // template <class BinaryFunction>
178     // double
179     // reductionOp(const DataAbstract& left, BinaryFunction operation,
180     // double initial_value);
181    
182     // implicit offset
183    
184     /**
185     \brief
186     Perform the given data point reduction operation on the data point
187     specified by the given offset into the view. Reduces all elements of
188     the data point using the given operation, returning the result as a
189     scalar. Operation must be a pointer to a function.
190    
191     Called by escript::algorithm.
192    
193 jfenwick 1779 \param left - vector containing the datapoint
194     \param shape - shape of datapoints in the vector
195     \param offset - beginning of datapoint in the vector
196 jfenwick 1724 \param operation - Input -
197     Operation to apply. Must be a pointer to a function.
198 jfenwick 1779 \param initial_value
199 jfenwick 1724 */
200     template <class BinaryFunction>
201     double
202 jfenwick 1779 reductionOp(const DataTypes::ValueType& left,
203     const DataTypes::ShapeType& shape,
204 jfenwick 1724 DataTypes::ValueType::size_type offset,
205     BinaryFunction operation,
206     double initial_value);
207    
208     /**
209     \brief
210     Perform a matrix multiply of the given views.
211    
212 jfenwick 1779 NB: Only multiplies together the two given datapoints,
213 jfenwick 1724 would need to call this over all data-points to multiply the entire
214     Data objects involved.
215    
216 jfenwick 1779 \param left,right - vectors containing the datapoints
217     \param leftShape,rightShape - shapes of datapoints in the vectors
218     \param leftOffset,rightOffset - beginnings of datapoints in the vectors
219     \param result - Vector to store the resulting datapoint in
220     \param resultShape - expected shape of the resulting datapoint
221 jfenwick 1724 */
222     ESCRIPT_DLL_API
223     void
224 jfenwick 1779 matMult(const DataTypes::ValueType& left,
225     const DataTypes::ShapeType& leftShape,
226     DataTypes::ValueType::size_type leftOffset,
227     const DataTypes::ValueType& right,
228     const DataTypes::ShapeType& rightShape,
229     DataTypes::ValueType::size_type rightOffset,
230 jfenwick 1734 DataTypes::ValueType& result,
231     const DataTypes::ShapeType& resultShape);
232 jfenwick 1724 // Hmmmm why is there no offset for the result??
233    
234    
235    
236    
237     /**
238     \brief
239     Determine the shape of the result array for a matrix multiplication
240     of the given views.
241 jfenwick 1779
242     \param left,right - shapes of the left and right matricies
243     \return the shape of the matrix which would result from multiplying left and right
244 jfenwick 1724 */
245     ESCRIPT_DLL_API
246     DataTypes::ShapeType
247     determineResultShape(const DataTypes::ShapeType& left,
248     const DataTypes::ShapeType& right);
249    
250     /**
251     \brief
252     computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
253    
254 jfenwick 1779 \param in - vector containing the matrix A
255     \param inShape - shape of the matrix A
256     \param inOffset - the beginning of A within the vector in
257     \param ev - vector to store the output matrix
258     \param evShape - expected shape of the output matrix
259     \param evOffset - starting location for storing ev in vector ev
260 jfenwick 1724 */
261     ESCRIPT_DLL_API
262     inline
263     void
264 jfenwick 1779 symmetric(const DataTypes::ValueType& in,
265     const DataTypes::ShapeType& inShape,
266 jfenwick 1724 DataTypes::ValueType::size_type inOffset,
267 jfenwick 1779 DataTypes::ValueType& ev,
268     const DataTypes::ShapeType& evShape,
269 jfenwick 1724 DataTypes::ValueType::size_type evOffset)
270     {
271     if (DataTypes::getRank(inShape) == 2) {
272     int i0, i1;
273     int s0=inShape[0];
274     int s1=inShape[1];
275     for (i0=0; i0<s0; i0++) {
276     for (i1=0; i1<s1; i1++) {
277     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
278     }
279     }
280     }
281     else if (DataTypes::getRank(inShape) == 4) {
282     int i0, i1, i2, i3;
283     int s0=inShape[0];
284     int s1=inShape[1];
285     int s2=inShape[2];
286     int s3=inShape[3];
287     for (i0=0; i0<s0; i0++) {
288     for (i1=0; i1<s1; i1++) {
289     for (i2=0; i2<s2; i2++) {
290     for (i3=0; i3<s3; i3++) {
291     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] + in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
292     }
293     }
294     }
295     }
296     }
297     }
298    
299     /**
300     \brief
301     computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
302    
303 jfenwick 1779 \param in - vector containing the matrix A
304     \param inShape - shape of the matrix A
305     \param inOffset - the beginning of A within the vector in
306     \param ev - vector to store the output matrix
307     \param evShape - expected shape of the output matrix
308     \param evOffset - starting location for storing ev in vector ev
309 jfenwick 1724 */
310     ESCRIPT_DLL_API
311     inline
312     void
313 jfenwick 1779 nonsymmetric(const DataTypes::ValueType& in,
314     const DataTypes::ShapeType& inShape,
315     DataTypes::ValueType::size_type inOffset,
316     DataTypes::ValueType& ev,
317     const DataTypes::ShapeType& evShape,
318     DataTypes::ValueType::size_type evOffset)
319 jfenwick 1724 {
320     if (DataTypes::getRank(inShape) == 2) {
321     int i0, i1;
322     int s0=inShape[0];
323     int s1=inShape[1];
324     for (i0=0; i0<s0; i0++) {
325     for (i1=0; i1<s1; i1++) {
326     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
327     }
328     }
329     }
330     else if (DataTypes::getRank(inShape) == 4) {
331     int i0, i1, i2, i3;
332     int s0=inShape[0];
333     int s1=inShape[1];
334     int s2=inShape[2];
335     int s3=inShape[3];
336     for (i0=0; i0<s0; i0++) {
337     for (i1=0; i1<s1; i1++) {
338     for (i2=0; i2<s2; i2++) {
339     for (i3=0; i3<s3; i3++) {
340     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] - in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
341     }
342     }
343     }
344     }
345     }
346     }
347    
348     /**
349     \brief
350     computes the trace of a matrix
351    
352 jfenwick 1779 \param in - vector containing the input matrix
353     \param inShape - shape of the input matrix
354     \param inOffset - the beginning of the input matrix within the vector "in"
355     \param ev - vector to store the output matrix
356     \param evShape - expected shape of the output matrix
357     \param evOffset - starting location for storing the output matrix in vector ev
358     \param axis_offset
359 jfenwick 1724 */
360     inline
361     void
362 jfenwick 1779 trace(const DataTypes::ValueType& in,
363     const DataTypes::ShapeType& inShape,
364 jfenwick 1724 DataTypes::ValueType::size_type inOffset,
365 jfenwick 1779 DataTypes::ValueType& ev,
366     const DataTypes::ShapeType& evShape,
367 jfenwick 1724 DataTypes::ValueType::size_type evOffset,
368     int axis_offset)
369     {
370     if (DataTypes::getRank(inShape) == 2) {
371     int s0=inShape[0]; // Python wrapper limits to square matrix
372     int i;
373     for (i=0; i<s0; i++) {
374     ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
375     }
376     }
377     else if (DataTypes::getRank(inShape) == 3) {
378     if (axis_offset==0) {
379     int s0=inShape[0];
380     int s2=inShape[2];
381     int i0, i2;
382     for (i0=0; i0<s0; i0++) {
383     for (i2=0; i2<s2; i2++) {
384     ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
385     }
386     }
387     }
388     else if (axis_offset==1) {
389     int s0=inShape[0];
390     int s1=inShape[1];
391     int i0, i1;
392     for (i0=0; i0<s0; i0++) {
393     for (i1=0; i1<s1; i1++) {
394     ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
395     }
396     }
397     }
398     }
399     else if (DataTypes::getRank(inShape) == 4) {
400     if (axis_offset==0) {
401     int s0=inShape[0];
402     int s2=inShape[2];
403     int s3=inShape[3];
404     int i0, i2, i3;
405     for (i0=0; i0<s0; i0++) {
406     for (i2=0; i2<s2; i2++) {
407     for (i3=0; i3<s3; i3++) {
408     ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
409     }
410     }
411     }
412     }
413     else if (axis_offset==1) {
414     int s0=inShape[0];
415     int s1=inShape[1];
416     int s3=inShape[3];
417     int i0, i1, i3;
418     for (i0=0; i0<s0; i0++) {
419     for (i1=0; i1<s1; i1++) {
420     for (i3=0; i3<s3; i3++) {
421     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
422     }
423     }
424     }
425     }
426     else if (axis_offset==2) {
427     int s0=inShape[0];
428     int s1=inShape[1];
429     int s2=inShape[2];
430     int i0, i1, i2;
431     for (i0=0; i0<s0; i0++) {
432     for (i1=0; i1<s1; i1++) {
433     for (i2=0; i2<s2; i2++) {
434     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
435     }
436     }
437     }
438     }
439     }
440     }
441    
442     /**
443     \brief
444     Transpose each data point of this Data object around the given axis.
445    
446 jfenwick 1779 \param in - vector containing the input matrix
447     \param inShape - shape of the input matrix
448     \param inOffset - the beginning of the input matrix within the vector "in"
449     \param ev - vector to store the output matrix
450     \param evShape - expected shape of the output matrix
451     \param evOffset - starting location for storing the output matrix in vector ev
452     \param axis_offset
453 jfenwick 1724 */
454     ESCRIPT_DLL_API
455     inline
456     void
457 jfenwick 1779 transpose(const DataTypes::ValueType& in,
458     const DataTypes::ShapeType& inShape,
459 jfenwick 1724 DataTypes::ValueType::size_type inOffset,
460 jfenwick 1779 DataTypes::ValueType& ev,
461     const DataTypes::ShapeType& evShape,
462 jfenwick 1724 DataTypes::ValueType::size_type evOffset,
463     int axis_offset)
464     {
465     int inRank=DataTypes::getRank(inShape);
466     if ( inRank== 4) {
467     int s0=evShape[0];
468     int s1=evShape[1];
469     int s2=evShape[2];
470     int s3=evShape[3];
471     int i0, i1, i2, i3;
472     if (axis_offset==1) {
473     for (i0=0; i0<s0; i0++) {
474     for (i1=0; i1<s1; i1++) {
475     for (i2=0; i2<s2; i2++) {
476     for (i3=0; i3<s3; i3++) {
477     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
478     }
479     }
480     }
481     }
482     }
483     else if (axis_offset==2) {
484     for (i0=0; i0<s0; i0++) {
485     for (i1=0; i1<s1; i1++) {
486     for (i2=0; i2<s2; i2++) {
487     for (i3=0; i3<s3; i3++) {
488     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
489     }
490     }
491     }
492     }
493     }
494     else if (axis_offset==3) {
495     for (i0=0; i0<s0; i0++) {
496     for (i1=0; i1<s1; i1++) {
497     for (i2=0; i2<s2; i2++) {
498     for (i3=0; i3<s3; i3++) {
499     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
500     }
501     }
502     }
503     }
504     }
505     else {
506     for (i0=0; i0<s0; i0++) {
507     for (i1=0; i1<s1; i1++) {
508     for (i2=0; i2<s2; i2++) {
509     for (i3=0; i3<s3; i3++) {
510     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
511     }
512     }
513     }
514     }
515     }
516     }
517     else if (inRank == 3) {
518     int s0=evShape[0];
519     int s1=evShape[1];
520     int s2=evShape[2];
521     int i0, i1, i2;
522     if (axis_offset==1) {
523     for (i0=0; i0<s0; i0++) {
524     for (i1=0; i1<s1; i1++) {
525     for (i2=0; i2<s2; i2++) {
526     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
527     }
528     }
529     }
530     }
531     else if (axis_offset==2) {
532     for (i0=0; i0<s0; i0++) {
533     for (i1=0; i1<s1; i1++) {
534     for (i2=0; i2<s2; i2++) {
535     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
536     }
537     }
538     }
539     }
540     else {
541     // Copy the matrix unchanged
542     for (i0=0; i0<s0; i0++) {
543     for (i1=0; i1<s1; i1++) {
544     for (i2=0; i2<s2; i2++) {
545     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
546     }
547     }
548     }
549     }
550     }
551     else if (inRank == 2) {
552     int s0=evShape[0];
553     int s1=evShape[1];
554     int i0, i1;
555     if (axis_offset==1) {
556     for (i0=0; i0<s0; i0++) {
557     for (i1=0; i1<s1; i1++) {
558     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
559     }
560     }
561     }
562     else {
563     for (i0=0; i0<s0; i0++) {
564     for (i1=0; i1<s1; i1++) {
565     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
566     }
567     }
568     }
569     }
570     else if (inRank == 1) {
571     int s0=evShape[0];
572     int i0;
573     for (i0=0; i0<s0; i0++) {
574     ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
575     }
576     }
577     else if (inRank == 0) {
578     ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
579     }
580     else {
581     throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
582     }
583     }
584    
585     /**
586     \brief
587     swaps the components axis0 and axis1.
588    
589 jfenwick 1779 \param in - vector containing the input matrix
590     \param inShape - shape of the input matrix
591     \param inOffset - the beginning of the input matrix within the vector "in"
592     \param ev - vector to store the output matrix
593     \param evShape - expected shape of the output matrix
594     \param evOffset - starting location for storing the output matrix in vector ev
595 jfenwick 1724 \param axis0 - axis index
596     \param axis1 - axis index
597     */
598     ESCRIPT_DLL_API
599     inline
600     void
601 jfenwick 1779 swapaxes(const DataTypes::ValueType& in,
602     const DataTypes::ShapeType& inShape,
603 jfenwick 1724 DataTypes::ValueType::size_type inOffset,
604 jfenwick 1779 DataTypes::ValueType& ev,
605     const DataTypes::ShapeType& evShape,
606 jfenwick 1724 DataTypes::ValueType::size_type evOffset,
607 jfenwick 1779 int axis0,
608     int axis1)
609 jfenwick 1724 {
610     int inRank=DataTypes::getRank(inShape);
611     if (inRank == 4) {
612     int s0=evShape[0];
613     int s1=evShape[1];
614     int s2=evShape[2];
615     int s3=evShape[3];
616     int i0, i1, i2, i3;
617     if (axis0==0) {
618     if (axis1==1) {
619     for (i0=0; i0<s0; i0++) {
620     for (i1=0; i1<s1; i1++) {
621     for (i2=0; i2<s2; i2++) {
622     for (i3=0; i3<s3; i3++) {
623     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
624     }
625     }
626     }
627     }
628     } else if (axis1==2) {
629     for (i0=0; i0<s0; i0++) {
630     for (i1=0; i1<s1; i1++) {
631     for (i2=0; i2<s2; i2++) {
632     for (i3=0; i3<s3; i3++) {
633     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
634     }
635     }
636     }
637     }
638    
639     } else if (axis1==3) {
640     for (i0=0; i0<s0; i0++) {
641     for (i1=0; i1<s1; i1++) {
642     for (i2=0; i2<s2; i2++) {
643     for (i3=0; i3<s3; i3++) {
644     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
645     }
646     }
647     }
648     }
649     }
650     } else if (axis0==1) {
651     if (axis1==2) {
652     for (i0=0; i0<s0; i0++) {
653     for (i1=0; i1<s1; i1++) {
654     for (i2=0; i2<s2; i2++) {
655     for (i3=0; i3<s3; i3++) {
656     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
657     }
658     }
659     }
660     }
661     } else if (axis1==3) {
662     for (i0=0; i0<s0; i0++) {
663     for (i1=0; i1<s1; i1++) {
664     for (i2=0; i2<s2; i2++) {
665     for (i3=0; i3<s3; i3++) {
666     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
667     }
668     }
669     }
670     }
671     }
672     } else if (axis0==2) {
673     if (axis1==3) {
674     for (i0=0; i0<s0; i0++) {
675     for (i1=0; i1<s1; i1++) {
676     for (i2=0; i2<s2; i2++) {
677     for (i3=0; i3<s3; i3++) {
678     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
679     }
680     }
681     }
682     }
683     }
684     }
685    
686     } else if ( inRank == 3) {
687     int s0=evShape[0];
688     int s1=evShape[1];
689     int s2=evShape[2];
690     int i0, i1, i2;
691     if (axis0==0) {
692     if (axis1==1) {
693     for (i0=0; i0<s0; i0++) {
694     for (i1=0; i1<s1; i1++) {
695     for (i2=0; i2<s2; i2++) {
696     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
697     }
698     }
699     }
700     } else if (axis1==2) {
701     for (i0=0; i0<s0; i0++) {
702     for (i1=0; i1<s1; i1++) {
703     for (i2=0; i2<s2; i2++) {
704     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
705     }
706     }
707     }
708     }
709     } else if (axis0==1) {
710     if (axis1==2) {
711     for (i0=0; i0<s0; i0++) {
712     for (i1=0; i1<s1; i1++) {
713     for (i2=0; i2<s2; i2++) {
714     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
715     }
716     }
717     }
718     }
719     }
720     } else if ( inRank == 2) {
721     int s0=evShape[0];
722     int s1=evShape[1];
723     int i0, i1;
724     if (axis0==0) {
725     if (axis1==1) {
726     for (i0=0; i0<s0; i0++) {
727     for (i1=0; i1<s1; i1++) {
728     ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
729     }
730     }
731     }
732     }
733     } else {
734     throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
735     }
736     }
737    
738     /**
739     \brief
740     solves a local eigenvalue problem
741    
742 jfenwick 1779 \param in - vector containing the input matrix
743     \param inShape - shape of the input matrix
744     \param inOffset - the beginning of the input matrix within the vector "in"
745     \param ev - vector to store the eigenvalues
746     \param evShape - expected shape of the eigenvalues
747     \param evOffset - starting location for storing the eigenvalues in vector ev
748 jfenwick 1724 */
749     ESCRIPT_DLL_API
750     inline
751     void
752 jfenwick 1779 eigenvalues(DataTypes::ValueType& in,
753     const DataTypes::ShapeType& inShape,
754 jfenwick 1724 DataTypes::ValueType::size_type inOffset,
755 jfenwick 1779 DataTypes::ValueType& ev,
756     const DataTypes::ShapeType& evShape,
757 jfenwick 1724 DataTypes::ValueType::size_type evOffset)
758     {
759     double in00,in10,in20,in01,in11,in21,in02,in12,in22;
760     double ev0,ev1,ev2;
761     int s=inShape[0];
762     if (s==1) {
763     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
764     eigenvalues1(in00,&ev0);
765     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
766    
767     } else if (s==2) {
768     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
769     in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
770     in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
771     in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
772     eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
773     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
774     ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
775    
776     } else if (s==3) {
777     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
778     in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
779     in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
780     in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
781     in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
782     in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
783     in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
784     in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
785     in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
786     eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
787     &ev0,&ev1,&ev2);
788     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
789     ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
790     ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
791    
792     }
793     }
794    
795     /**
796     \brief
797     solves a local eigenvalue problem
798    
799 jfenwick 1779 \param in - vector containing the input matrix
800     \param inShape - shape of the input matrix
801     \param inOffset - the beginning of the input matrix within the vector "in"
802     \param ev - vector to store the eigenvalues
803     \param evShape - expected shape of the eigenvalues
804     \param evOffset - starting location for storing the eigenvalues in ev
805     \param V - vector to store the eigenvectors
806     \param VShape - expected shape of the eigenvectors
807     \param VOffset - starting location for storing the eigenvectors in V
808 jfenwick 1724 \param tol - Input - eigenvalues with relative difference tol are treated as equal
809     */
810     ESCRIPT_DLL_API
811     inline
812     void
813     eigenvalues_and_eigenvectors(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
814     DataTypes::ValueType::size_type inOffset,
815     DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
816     DataTypes::ValueType::size_type evOffset,
817     DataTypes::ValueType& V, const DataTypes::ShapeType& VShape,
818     DataTypes::ValueType::size_type VOffset,
819     const double tol=1.e-13)
820     {
821     double in00,in10,in20,in01,in11,in21,in02,in12,in22;
822     double V00,V10,V20,V01,V11,V21,V02,V12,V22;
823     double ev0,ev1,ev2;
824     int s=inShape[0];
825     if (s==1) {
826     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
827     eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
828     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
829     V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
830     } else if (s==2) {
831     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
832     in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
833     in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
834     in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
835     eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
836     &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
837     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
838     ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
839     V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
840     V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
841     V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
842     V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
843     } else if (s==3) {
844     in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
845     in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
846     in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
847     in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
848     in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
849     in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
850     in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
851     in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
852     in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
853     eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
854     &ev0,&ev1,&ev2,
855     &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
856     ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
857     ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
858     ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
859     V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
860     V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
861     V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
862     V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
863     V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
864     V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
865     V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
866     V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
867     V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
868    
869     }
870     }
871    
872    
873     /**
874     Inline function definitions.
875     */
876    
877     // template <class UnaryFunction>
878     // inline
879     // void
880     // DataArrayView::unaryOp(UnaryFunction operation)
881     // {
882     // unaryOp(m_offset,operation);
883     // }
884    
885    
886     inline
887     bool
888     checkOffset(const DataTypes::ValueType& data,
889     const DataTypes::ShapeType& shape,
890     DataTypes::ValueType::size_type offset)
891     {
892     return (data.size() >= (offset+DataTypes::noValues(shape)));
893     }
894    
895     template <class UnaryFunction>
896     inline
897     void
898     unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
899     DataTypes::ValueType::size_type offset,
900     UnaryFunction operation)
901     {
902     EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
903     "Error - Couldn't perform unaryOp due to insufficient storage.");
904     DataTypes::ValueType::size_type nVals=DataTypes::noValues(shape);
905     for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
906     data[offset+i]=operation(data[offset+i]);
907     }
908     }
909    
910     // template <class BinaryFunction>
911     // inline
912     // void
913     // binaryOp(const DataArrayView& right,
914     // BinaryFunction operation)
915     // {
916     // binaryOp(m_offset,right,right.getOffset(),operation);
917     // }
918    
919    
920     template <class BinaryFunction>
921     inline
922     void
923     binaryOp(DataTypes::ValueType& left,
924     const DataTypes::ShapeType& leftShape,
925     DataTypes::ValueType::size_type leftOffset,
926     const DataTypes::ValueType& right,
927     const DataTypes::ShapeType& rightShape,
928     DataTypes::ValueType::size_type rightOffset,
929     BinaryFunction operation)
930     {
931     EsysAssert(leftShape==rightShape,
932     "Error - Couldn't perform binaryOp due to shape mismatch,");
933     EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
934     "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
935     EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
936     "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
937     for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
938     left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
939     }
940     }
941    
942     template <class BinaryFunction>
943     inline
944     void
945     binaryOp(DataTypes::ValueType& left,
946     const DataTypes::ShapeType& leftShape,
947     DataTypes::ValueType::size_type offset,
948     double right,
949     BinaryFunction operation)
950     {
951     EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
952     "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
953     for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
954     left[offset+i]=operation(left[offset+i],right);
955     }
956     }
957    
958     template <class BinaryFunction>
959     inline
960     double
961     reductionOp(const DataTypes::ValueType& left,
962     const DataTypes::ShapeType& leftShape,
963     DataTypes::ValueType::size_type offset,
964     BinaryFunction operation,
965     double initial_value)
966     {
967     EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
968     "Error - Couldn't perform reductionOp due to insufficient storage.");
969     double current_value=initial_value;
970     for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
971     current_value=operation(current_value,left[offset+i]);
972     }
973     return current_value;
974     }
975    
976    
977    
978     } // end namespace DataMath
979     } // end namespace escript
980     #endif
981    

  ViewVC Help
Powered by ViewVC 1.1.26