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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26