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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (10 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 31477 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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

  ViewVC Help
Powered by ViewVC 1.1.26