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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2271 - (show annotations)
Mon Feb 16 05:08:29 2009 UTC (10 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 31569 byte(s)
Merging version 2269 to trunk

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 for (int j=0;j<DataTypes::noValues(evShape);++j)
300 {
301 ev[evOffset+j]=0;
302 }
303 if (DataTypes::getRank(inShape) == 2) {
304 int s0=inShape[0]; // Python wrapper limits to square matrix
305 int i;
306 for (i=0; i<s0; i++) {
307 ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
308 }
309 }
310 else if (DataTypes::getRank(inShape) == 3) {
311 if (axis_offset==0) {
312 int s0=inShape[0];
313 int s2=inShape[2];
314 int i0, i2;
315 for (i0=0; i0<s0; i0++) {
316 for (i2=0; i2<s2; i2++) {
317 ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
318 }
319 }
320 }
321 else if (axis_offset==1) {
322 int s0=inShape[0];
323 int s1=inShape[1];
324 int i0, i1;
325 for (i0=0; i0<s0; i0++) {
326 for (i1=0; i1<s1; i1++) {
327 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
328 }
329 }
330 }
331 }
332 else if (DataTypes::getRank(inShape) == 4) {
333 if (axis_offset==0) {
334 int s0=inShape[0];
335 int s2=inShape[2];
336 int s3=inShape[3];
337 int i0, i2, i3;
338 for (i0=0; i0<s0; i0++) {
339 for (i2=0; i2<s2; i2++) {
340 for (i3=0; i3<s3; i3++) {
341 ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
342 }
343 }
344 }
345 }
346 else if (axis_offset==1) {
347 int s0=inShape[0];
348 int s1=inShape[1];
349 int s3=inShape[3];
350 int i0, i1, i3;
351 for (i0=0; i0<s0; i0++) {
352 for (i1=0; i1<s1; i1++) {
353 for (i3=0; i3<s3; i3++) {
354 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
355 }
356 }
357 }
358 }
359 else if (axis_offset==2) {
360 int s0=inShape[0];
361 int s1=inShape[1];
362 int s2=inShape[2];
363 int i0, i1, i2;
364 for (i0=0; i0<s0; i0++) {
365 for (i1=0; i1<s1; i1++) {
366 for (i2=0; i2<s2; i2++) {
367 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
368 }
369 }
370 }
371 }
372 }
373 }
374
375 /**
376 \brief
377 Transpose each data point of this Data object around the given axis.
378
379 \param in - vector containing the input matrix
380 \param inShape - shape of the input matrix
381 \param inOffset - the beginning of the input matrix within the vector "in"
382 \param ev - vector to store the output matrix
383 \param evShape - expected shape of the output matrix
384 \param evOffset - starting location for storing the output matrix in vector ev
385 \param axis_offset
386 */
387 ESCRIPT_DLL_API
388 inline
389 void
390 transpose(const DataTypes::ValueType& in,
391 const DataTypes::ShapeType& inShape,
392 DataTypes::ValueType::size_type inOffset,
393 DataTypes::ValueType& ev,
394 const DataTypes::ShapeType& evShape,
395 DataTypes::ValueType::size_type evOffset,
396 int axis_offset)
397 {
398 int inRank=DataTypes::getRank(inShape);
399 if ( inRank== 4) {
400 int s0=evShape[0];
401 int s1=evShape[1];
402 int s2=evShape[2];
403 int s3=evShape[3];
404 int i0, i1, i2, i3;
405 if (axis_offset==1) {
406 for (i0=0; i0<s0; i0++) {
407 for (i1=0; i1<s1; i1++) {
408 for (i2=0; i2<s2; i2++) {
409 for (i3=0; i3<s3; i3++) {
410 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
411 }
412 }
413 }
414 }
415 }
416 else if (axis_offset==2) {
417 for (i0=0; i0<s0; i0++) {
418 for (i1=0; i1<s1; i1++) {
419 for (i2=0; i2<s2; i2++) {
420 for (i3=0; i3<s3; i3++) {
421 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
422 }
423 }
424 }
425 }
426 }
427 else if (axis_offset==3) {
428 for (i0=0; i0<s0; i0++) {
429 for (i1=0; i1<s1; i1++) {
430 for (i2=0; i2<s2; i2++) {
431 for (i3=0; i3<s3; i3++) {
432 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
433 }
434 }
435 }
436 }
437 }
438 else {
439 for (i0=0; i0<s0; i0++) {
440 for (i1=0; i1<s1; i1++) {
441 for (i2=0; i2<s2; i2++) {
442 for (i3=0; i3<s3; i3++) {
443 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
444 }
445 }
446 }
447 }
448 }
449 }
450 else if (inRank == 3) {
451 int s0=evShape[0];
452 int s1=evShape[1];
453 int s2=evShape[2];
454 int i0, i1, i2;
455 if (axis_offset==1) {
456 for (i0=0; i0<s0; i0++) {
457 for (i1=0; i1<s1; i1++) {
458 for (i2=0; i2<s2; i2++) {
459 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
460 }
461 }
462 }
463 }
464 else if (axis_offset==2) {
465 for (i0=0; i0<s0; i0++) {
466 for (i1=0; i1<s1; i1++) {
467 for (i2=0; i2<s2; i2++) {
468 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
469 }
470 }
471 }
472 }
473 else {
474 // Copy the matrix unchanged
475 for (i0=0; i0<s0; i0++) {
476 for (i1=0; i1<s1; i1++) {
477 for (i2=0; i2<s2; i2++) {
478 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
479 }
480 }
481 }
482 }
483 }
484 else if (inRank == 2) {
485 int s0=evShape[0];
486 int s1=evShape[1];
487 int i0, i1;
488 if (axis_offset==1) {
489 for (i0=0; i0<s0; i0++) {
490 for (i1=0; i1<s1; i1++) {
491 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
492 }
493 }
494 }
495 else {
496 for (i0=0; i0<s0; i0++) {
497 for (i1=0; i1<s1; i1++) {
498 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
499 }
500 }
501 }
502 }
503 else if (inRank == 1) {
504 int s0=evShape[0];
505 int i0;
506 for (i0=0; i0<s0; i0++) {
507 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
508 }
509 }
510 else if (inRank == 0) {
511 ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
512 }
513 else {
514 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
515 }
516 }
517
518 /**
519 \brief
520 swaps the components axis0 and axis1.
521
522 \param in - vector containing the input matrix
523 \param inShape - shape of the input matrix
524 \param inOffset - the beginning of the input matrix within the vector "in"
525 \param ev - vector to store the output matrix
526 \param evShape - expected shape of the output matrix
527 \param evOffset - starting location for storing the output matrix in vector ev
528 \param axis0 - axis index
529 \param axis1 - axis index
530 */
531 ESCRIPT_DLL_API
532 inline
533 void
534 swapaxes(const DataTypes::ValueType& in,
535 const DataTypes::ShapeType& inShape,
536 DataTypes::ValueType::size_type inOffset,
537 DataTypes::ValueType& ev,
538 const DataTypes::ShapeType& evShape,
539 DataTypes::ValueType::size_type evOffset,
540 int axis0,
541 int axis1)
542 {
543 int inRank=DataTypes::getRank(inShape);
544 if (inRank == 4) {
545 int s0=evShape[0];
546 int s1=evShape[1];
547 int s2=evShape[2];
548 int s3=evShape[3];
549 int i0, i1, i2, i3;
550 if (axis0==0) {
551 if (axis1==1) {
552 for (i0=0; i0<s0; i0++) {
553 for (i1=0; i1<s1; i1++) {
554 for (i2=0; i2<s2; i2++) {
555 for (i3=0; i3<s3; i3++) {
556 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
557 }
558 }
559 }
560 }
561 } else if (axis1==2) {
562 for (i0=0; i0<s0; i0++) {
563 for (i1=0; i1<s1; i1++) {
564 for (i2=0; i2<s2; i2++) {
565 for (i3=0; i3<s3; i3++) {
566 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
567 }
568 }
569 }
570 }
571
572 } else if (axis1==3) {
573 for (i0=0; i0<s0; i0++) {
574 for (i1=0; i1<s1; i1++) {
575 for (i2=0; i2<s2; i2++) {
576 for (i3=0; i3<s3; i3++) {
577 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
578 }
579 }
580 }
581 }
582 }
583 } else if (axis0==1) {
584 if (axis1==2) {
585 for (i0=0; i0<s0; i0++) {
586 for (i1=0; i1<s1; i1++) {
587 for (i2=0; i2<s2; i2++) {
588 for (i3=0; i3<s3; i3++) {
589 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
590 }
591 }
592 }
593 }
594 } else if (axis1==3) {
595 for (i0=0; i0<s0; i0++) {
596 for (i1=0; i1<s1; i1++) {
597 for (i2=0; i2<s2; i2++) {
598 for (i3=0; i3<s3; i3++) {
599 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
600 }
601 }
602 }
603 }
604 }
605 } else if (axis0==2) {
606 if (axis1==3) {
607 for (i0=0; i0<s0; i0++) {
608 for (i1=0; i1<s1; i1++) {
609 for (i2=0; i2<s2; i2++) {
610 for (i3=0; i3<s3; i3++) {
611 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
612 }
613 }
614 }
615 }
616 }
617 }
618
619 } else if ( inRank == 3) {
620 int s0=evShape[0];
621 int s1=evShape[1];
622 int s2=evShape[2];
623 int i0, i1, i2;
624 if (axis0==0) {
625 if (axis1==1) {
626 for (i0=0; i0<s0; i0++) {
627 for (i1=0; i1<s1; i1++) {
628 for (i2=0; i2<s2; i2++) {
629 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
630 }
631 }
632 }
633 } else if (axis1==2) {
634 for (i0=0; i0<s0; i0++) {
635 for (i1=0; i1<s1; i1++) {
636 for (i2=0; i2<s2; i2++) {
637 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
638 }
639 }
640 }
641 }
642 } else if (axis0==1) {
643 if (axis1==2) {
644 for (i0=0; i0<s0; i0++) {
645 for (i1=0; i1<s1; i1++) {
646 for (i2=0; i2<s2; i2++) {
647 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
648 }
649 }
650 }
651 }
652 }
653 } else if ( inRank == 2) {
654 int s0=evShape[0];
655 int s1=evShape[1];
656 int i0, i1;
657 if (axis0==0) {
658 if (axis1==1) {
659 for (i0=0; i0<s0; i0++) {
660 for (i1=0; i1<s1; i1++) {
661 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
662 }
663 }
664 }
665 }
666 } else {
667 throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
668 }
669 }
670
671 /**
672 \brief
673 solves a local eigenvalue problem
674
675 \param in - vector containing the input matrix
676 \param inShape - shape of the input matrix
677 \param inOffset - the beginning of the input matrix within the vector "in"
678 \param ev - vector to store the eigenvalues
679 \param evShape - expected shape of the eigenvalues
680 \param evOffset - starting location for storing the eigenvalues in vector ev
681 */
682 ESCRIPT_DLL_API
683 inline
684 void
685 eigenvalues(const DataTypes::ValueType& in,
686 const DataTypes::ShapeType& inShape,
687 DataTypes::ValueType::size_type inOffset,
688 DataTypes::ValueType& ev,
689 const DataTypes::ShapeType& evShape,
690 DataTypes::ValueType::size_type evOffset)
691 {
692 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
693 double ev0,ev1,ev2;
694 int s=inShape[0];
695 if (s==1) {
696 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
697 eigenvalues1(in00,&ev0);
698 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
699
700 } else if (s==2) {
701 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
702 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
703 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
704 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
705 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
706 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
707 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
708
709 } else if (s==3) {
710 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
711 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
712 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
713 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
714 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
715 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
716 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
717 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
718 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
719 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
720 &ev0,&ev1,&ev2);
721 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
722 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
723 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
724
725 }
726 }
727
728 /**
729 \brief
730 solves a local eigenvalue problem
731
732 \param in - vector containing the input matrix
733 \param inShape - shape of the input matrix
734 \param inOffset - the beginning of the input matrix within the vector "in"
735 \param ev - vector to store the eigenvalues
736 \param evShape - expected shape of the eigenvalues
737 \param evOffset - starting location for storing the eigenvalues in ev
738 \param V - vector to store the eigenvectors
739 \param VShape - expected shape of the eigenvectors
740 \param VOffset - starting location for storing the eigenvectors in V
741 \param tol - Input - eigenvalues with relative difference tol are treated as equal
742 */
743 ESCRIPT_DLL_API
744 inline
745 void
746 eigenvalues_and_eigenvectors(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
747 DataTypes::ValueType::size_type inOffset,
748 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
749 DataTypes::ValueType::size_type evOffset,
750 DataTypes::ValueType& V, const DataTypes::ShapeType& VShape,
751 DataTypes::ValueType::size_type VOffset,
752 const double tol=1.e-13)
753 {
754 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
755 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
756 double ev0,ev1,ev2;
757 int s=inShape[0];
758 if (s==1) {
759 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
760 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
761 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
762 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
763 } else if (s==2) {
764 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
765 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
766 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
767 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
768 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
769 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
770 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
771 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
772 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
773 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
774 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
775 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
776 } else if (s==3) {
777 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
778 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
779 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
780 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
781 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
782 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
783 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
784 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
785 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
786 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
787 &ev0,&ev1,&ev2,
788 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
789 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
790 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
791 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
792 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
793 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
794 V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
795 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
796 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
797 V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
798 V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
799 V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
800 V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
801
802 }
803 }
804
805
806 /**
807 Inline function definitions.
808 */
809
810 inline
811 bool
812 checkOffset(const DataTypes::ValueType& data,
813 const DataTypes::ShapeType& shape,
814 DataTypes::ValueType::size_type offset)
815 {
816 return (data.size() >= (offset+DataTypes::noValues(shape)));
817 }
818
819 template <class UnaryFunction>
820 inline
821 void
822 unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
823 DataTypes::ValueType::size_type offset,
824 UnaryFunction operation)
825 {
826 EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
827 "Error - Couldn't perform unaryOp due to insufficient storage.");
828 DataTypes::ValueType::size_type nVals=DataTypes::noValues(shape);
829 for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
830 data[offset+i]=operation(data[offset+i]);
831 }
832 }
833
834
835 template <class BinaryFunction>
836 inline
837 void
838 binaryOp(DataTypes::ValueType& left,
839 const DataTypes::ShapeType& leftShape,
840 DataTypes::ValueType::size_type leftOffset,
841 const DataTypes::ValueType& right,
842 const DataTypes::ShapeType& rightShape,
843 DataTypes::ValueType::size_type rightOffset,
844 BinaryFunction operation)
845 {
846 EsysAssert(leftShape==rightShape,
847 "Error - Couldn't perform binaryOp due to shape mismatch,");
848 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
849 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
850 EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
851 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
852 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
853 left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
854 }
855 }
856
857 template <class BinaryFunction>
858 inline
859 void
860 binaryOp(DataTypes::ValueType& left,
861 const DataTypes::ShapeType& leftShape,
862 DataTypes::ValueType::size_type offset,
863 double right,
864 BinaryFunction operation)
865 {
866 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
867 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
868 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
869 left[offset+i]=operation(left[offset+i],right);
870 }
871 }
872
873 template <class BinaryFunction>
874 inline
875 double
876 reductionOp(const DataTypes::ValueType& left,
877 const DataTypes::ShapeType& leftShape,
878 DataTypes::ValueType::size_type offset,
879 BinaryFunction operation,
880 double initial_value)
881 {
882 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
883 "Error - Couldn't perform reductionOp due to insufficient storage.");
884 double current_value=initial_value;
885 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
886 current_value=operation(current_value,left[offset+i]);
887 }
888 return current_value;
889 }
890
891
892
893 } // end namespace DataMath
894 } // end namespace escript
895 #endif
896

  ViewVC Help
Powered by ViewVC 1.1.26