/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataMaths.h
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/DataMaths.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1734 - (show annotations)
Thu Aug 28 06:11:56 2008 UTC (11 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 32226 byte(s)
Added operations to Datatypes:
checkOffset
copySlice
copySliceFrom

Fixed some error reporting using EsysAssert.

Added two new c++ test suites:
DataTypesTest
DataMathsTest

Note that the test suite does not compile with dodebug=yes. There is an issue with linking one of the exception functions. I'm going to leave this 
until I have finished the rest of the work, perhaps Ken's scons changes will fix it.


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

  ViewVC Help
Powered by ViewVC 1.1.26