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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1724 - (hide annotations)
Mon Aug 25 05:38:57 2008 UTC (11 years, 2 months ago) by jfenwick
Original Path: branches/arrayview_from_1695_trunk/escript/src/DataMaths.h
File MIME type: text/plain
File size: 32145 byte(s)
Branch commit

Moved createShapeErrorMessage() into DataTypes.h
Modified functions in DataAlgorithm.h to use non-DataArrayView accesses.

Added getVector() to each of DataTagged, DataConstant, DataExpanded - This method returns 
the underlying DataVector by reference/constant reference. Note that this method does not 
exist in DataAbstract so (at the momement) in order to pull the data from something you 
need to know what you are looking at. (Lower level access is still possible via double* 
though).

DataTagged now has a getOffsetForTag method and a getDefaultOffset method.

DataMaths.h - A new file containing the reductionOps etc from DataArrayView (but without 
requiring DAV).
This file requires significant commenting improvements.


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

  ViewVC Help
Powered by ViewVC 1.1.26