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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 33249 byte(s)
Don't panic.
Updating copyright stamps

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

  ViewVC Help
Powered by ViewVC 1.1.26