/[escript]/branches/trilinos_from_5897/escriptcore/src/DataMaths.h
ViewVC logotype

Contents of /branches/trilinos_from_5897/escriptcore/src/DataMaths.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5963 - (show annotations)
Mon Feb 22 06:59:27 2016 UTC (3 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 40454 byte(s)
sync and fix.

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

  ViewVC Help
Powered by ViewVC 1.1.26