/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataMaths.h
ViewVC logotype

Annotation of /branches/arrayview_from_1695_trunk/escript/src/DataMaths.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1779 - (hide annotations)
Thu Sep 11 01:06:15 2008 UTC (11 years, 2 months ago) by jfenwick
File MIME type: text/plain
File size: 34492 byte(s)
Branch commit:

Updated the doxyfile for deprecated/modified options (eg encoding is set 
to UTF-8).
Also set it not to scan directories called .svn or deprecated.

Removed #includes of DataArrayView.h
Moved DataArrayView into a new deprecated directory

Commented out setTaggedValue in DataExpanded which still referred to 
DataArrayView

Fixed some of the doxygen comments in DataMaths.h, DataTypes.h, 
DataTagged.h

Removed the empty JoelMods.cpp_

Unit tests appear to indicate that this branch is now "no more broken" 
than the version I branched from.



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

  ViewVC Help
Powered by ViewVC 1.1.26