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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (show annotations)
Wed Sep 17 01:45:46 2008 UTC (10 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 34492 byte(s)
Merged noarrayview branch onto trunk.


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

  ViewVC Help
Powered by ViewVC 1.1.26