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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26