/[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 1726 - (show annotations)
Tue Aug 26 03:33:34 2008 UTC (11 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 32182 byte(s)
BinaryOp and UnaryOp - modified to accept the extra parameters required 
to operate without DataArrayView.  There are still a few parallel 
methods which accept them as params.

Added getVector() members to DataAbstract - these versions will throw.
DataC.cpp - uses new method to get the shape.
Added constant form of getVector() to DataConstant
Fixed the #include protection on DataMaths.h


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 // Hmmmm why is there no offset for the result??
222
223
224
225
226 /**
227 \brief
228 Determine the shape of the result array for a matrix multiplication
229 of the given views.
230 */
231 ESCRIPT_DLL_API
232 DataTypes::ShapeType
233 determineResultShape(const DataTypes::ShapeType& left,
234 const DataTypes::ShapeType& right);
235
236 /**
237 \brief
238 computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
239
240 \param in - Input - matrix
241 \param inOffset - Input - offset into in
242 \param ev - Output - The symmetric matrix
243 \param inOffset - Input - offset into ev
244 */
245 ESCRIPT_DLL_API
246 inline
247 void
248 symmetric(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
249 DataTypes::ValueType::size_type inOffset,
250 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
251 DataTypes::ValueType::size_type evOffset)
252 {
253 if (DataTypes::getRank(inShape) == 2) {
254 int i0, i1;
255 int s0=inShape[0];
256 int s1=inShape[1];
257 for (i0=0; i0<s0; i0++) {
258 for (i1=0; i1<s1; i1++) {
259 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
260 }
261 }
262 }
263 else if (DataTypes::getRank(inShape) == 4) {
264 int i0, i1, i2, i3;
265 int s0=inShape[0];
266 int s1=inShape[1];
267 int s2=inShape[2];
268 int s3=inShape[3];
269 for (i0=0; i0<s0; i0++) {
270 for (i1=0; i1<s1; i1++) {
271 for (i2=0; i2<s2; i2++) {
272 for (i3=0; i3<s3; i3++) {
273 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;
274 }
275 }
276 }
277 }
278 }
279 }
280
281 /**
282 \brief
283 computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
284
285 \param in - Input - matrix
286 \param inOffset - Input - offset into in
287 \param ev - Output - The nonsymmetric matrix
288 \param inOffset - Input - offset into ev
289 */
290 ESCRIPT_DLL_API
291 inline
292 void
293 nonsymmetric(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
294 DataTypes::ValueType::size_type inOffset,
295 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
296 DataTypes::ValueType::size_type evOffset)
297 {
298 if (DataTypes::getRank(inShape) == 2) {
299 int i0, i1;
300 int s0=inShape[0];
301 int s1=inShape[1];
302 for (i0=0; i0<s0; i0++) {
303 for (i1=0; i1<s1; i1++) {
304 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
305 }
306 }
307 }
308 else if (DataTypes::getRank(inShape) == 4) {
309 int i0, i1, i2, i3;
310 int s0=inShape[0];
311 int s1=inShape[1];
312 int s2=inShape[2];
313 int s3=inShape[3];
314 for (i0=0; i0<s0; i0++) {
315 for (i1=0; i1<s1; i1++) {
316 for (i2=0; i2<s2; i2++) {
317 for (i3=0; i3<s3; i3++) {
318 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;
319 }
320 }
321 }
322 }
323 }
324 }
325
326 /**
327 \brief
328 computes the trace of a matrix
329
330 \param in - Input - matrix
331 \param inOffset - Input - offset into in
332 \param ev - Output - The nonsymmetric matrix
333 \param inOffset - Input - offset into ev
334 */
335 inline
336 void
337 trace(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
338 DataTypes::ValueType::size_type inOffset,
339 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
340 DataTypes::ValueType::size_type evOffset,
341 int axis_offset)
342 {
343 if (DataTypes::getRank(inShape) == 2) {
344 int s0=inShape[0]; // Python wrapper limits to square matrix
345 int i;
346 for (i=0; i<s0; i++) {
347 ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
348 }
349 }
350 else if (DataTypes::getRank(inShape) == 3) {
351 if (axis_offset==0) {
352 int s0=inShape[0];
353 int s2=inShape[2];
354 int i0, i2;
355 for (i0=0; i0<s0; i0++) {
356 for (i2=0; i2<s2; i2++) {
357 ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
358 }
359 }
360 }
361 else if (axis_offset==1) {
362 int s0=inShape[0];
363 int s1=inShape[1];
364 int i0, i1;
365 for (i0=0; i0<s0; i0++) {
366 for (i1=0; i1<s1; i1++) {
367 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
368 }
369 }
370 }
371 }
372 else if (DataTypes::getRank(inShape) == 4) {
373 if (axis_offset==0) {
374 int s0=inShape[0];
375 int s2=inShape[2];
376 int s3=inShape[3];
377 int i0, i2, i3;
378 for (i0=0; i0<s0; i0++) {
379 for (i2=0; i2<s2; i2++) {
380 for (i3=0; i3<s3; i3++) {
381 ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
382 }
383 }
384 }
385 }
386 else if (axis_offset==1) {
387 int s0=inShape[0];
388 int s1=inShape[1];
389 int s3=inShape[3];
390 int i0, i1, i3;
391 for (i0=0; i0<s0; i0++) {
392 for (i1=0; i1<s1; i1++) {
393 for (i3=0; i3<s3; i3++) {
394 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
395 }
396 }
397 }
398 }
399 else if (axis_offset==2) {
400 int s0=inShape[0];
401 int s1=inShape[1];
402 int s2=inShape[2];
403 int i0, i1, i2;
404 for (i0=0; i0<s0; i0++) {
405 for (i1=0; i1<s1; i1++) {
406 for (i2=0; i2<s2; i2++) {
407 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
408 }
409 }
410 }
411 }
412 }
413 }
414
415 /**
416 \brief
417 Transpose each data point of this Data object around the given axis.
418
419 \param in - Input - matrix
420 \param inOffset - Input - offset into in
421 \param ev - Output - The nonsymmetric matrix
422 \param inOffset - Input - offset into ev
423 */
424 ESCRIPT_DLL_API
425 inline
426 void
427 transpose(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
428 DataTypes::ValueType::size_type inOffset,
429 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
430 DataTypes::ValueType::size_type evOffset,
431 int axis_offset)
432 {
433 int inRank=DataTypes::getRank(inShape);
434 if ( inRank== 4) {
435 int s0=evShape[0];
436 int s1=evShape[1];
437 int s2=evShape[2];
438 int s3=evShape[3];
439 int i0, i1, i2, i3;
440 if (axis_offset==1) {
441 for (i0=0; i0<s0; i0++) {
442 for (i1=0; i1<s1; i1++) {
443 for (i2=0; i2<s2; i2++) {
444 for (i3=0; i3<s3; i3++) {
445 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
446 }
447 }
448 }
449 }
450 }
451 else if (axis_offset==2) {
452 for (i0=0; i0<s0; i0++) {
453 for (i1=0; i1<s1; i1++) {
454 for (i2=0; i2<s2; i2++) {
455 for (i3=0; i3<s3; i3++) {
456 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
457 }
458 }
459 }
460 }
461 }
462 else if (axis_offset==3) {
463 for (i0=0; i0<s0; i0++) {
464 for (i1=0; i1<s1; i1++) {
465 for (i2=0; i2<s2; i2++) {
466 for (i3=0; i3<s3; i3++) {
467 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
468 }
469 }
470 }
471 }
472 }
473 else {
474 for (i0=0; i0<s0; i0++) {
475 for (i1=0; i1<s1; i1++) {
476 for (i2=0; i2<s2; i2++) {
477 for (i3=0; i3<s3; i3++) {
478 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
479 }
480 }
481 }
482 }
483 }
484 }
485 else if (inRank == 3) {
486 int s0=evShape[0];
487 int s1=evShape[1];
488 int s2=evShape[2];
489 int i0, i1, i2;
490 if (axis_offset==1) {
491 for (i0=0; i0<s0; i0++) {
492 for (i1=0; i1<s1; i1++) {
493 for (i2=0; i2<s2; i2++) {
494 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
495 }
496 }
497 }
498 }
499 else if (axis_offset==2) {
500 for (i0=0; i0<s0; i0++) {
501 for (i1=0; i1<s1; i1++) {
502 for (i2=0; i2<s2; i2++) {
503 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
504 }
505 }
506 }
507 }
508 else {
509 // Copy the matrix unchanged
510 for (i0=0; i0<s0; i0++) {
511 for (i1=0; i1<s1; i1++) {
512 for (i2=0; i2<s2; i2++) {
513 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
514 }
515 }
516 }
517 }
518 }
519 else if (inRank == 2) {
520 int s0=evShape[0];
521 int s1=evShape[1];
522 int i0, i1;
523 if (axis_offset==1) {
524 for (i0=0; i0<s0; i0++) {
525 for (i1=0; i1<s1; i1++) {
526 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
527 }
528 }
529 }
530 else {
531 for (i0=0; i0<s0; i0++) {
532 for (i1=0; i1<s1; i1++) {
533 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
534 }
535 }
536 }
537 }
538 else if (inRank == 1) {
539 int s0=evShape[0];
540 int i0;
541 for (i0=0; i0<s0; i0++) {
542 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
543 }
544 }
545 else if (inRank == 0) {
546 ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
547 }
548 else {
549 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
550 }
551 }
552
553 /**
554 \brief
555 swaps the components axis0 and axis1.
556
557 \param in - Input - matrix
558 \param inOffset - Input - offset into in
559 \param ev - Output - The nonsymmetric matrix
560 \param inOffset - Input - offset into ev
561 \param axis0 - axis index
562 \param axis1 - axis index
563 */
564 ESCRIPT_DLL_API
565 static
566 inline
567 void
568 swapaxes(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
569 DataTypes::ValueType::size_type inOffset,
570 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
571 DataTypes::ValueType::size_type evOffset,
572 int axis0, int axis1)
573 {
574 int inRank=DataTypes::getRank(inShape);
575 if (inRank == 4) {
576 int s0=evShape[0];
577 int s1=evShape[1];
578 int s2=evShape[2];
579 int s3=evShape[3];
580 int i0, i1, i2, i3;
581 if (axis0==0) {
582 if (axis1==1) {
583 for (i0=0; i0<s0; i0++) {
584 for (i1=0; i1<s1; i1++) {
585 for (i2=0; i2<s2; i2++) {
586 for (i3=0; i3<s3; i3++) {
587 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
588 }
589 }
590 }
591 }
592 } else if (axis1==2) {
593 for (i0=0; i0<s0; i0++) {
594 for (i1=0; i1<s1; i1++) {
595 for (i2=0; i2<s2; i2++) {
596 for (i3=0; i3<s3; i3++) {
597 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
598 }
599 }
600 }
601 }
602
603 } else if (axis1==3) {
604 for (i0=0; i0<s0; i0++) {
605 for (i1=0; i1<s1; i1++) {
606 for (i2=0; i2<s2; i2++) {
607 for (i3=0; i3<s3; i3++) {
608 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
609 }
610 }
611 }
612 }
613 }
614 } else if (axis0==1) {
615 if (axis1==2) {
616 for (i0=0; i0<s0; i0++) {
617 for (i1=0; i1<s1; i1++) {
618 for (i2=0; i2<s2; i2++) {
619 for (i3=0; i3<s3; i3++) {
620 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
621 }
622 }
623 }
624 }
625 } else if (axis1==3) {
626 for (i0=0; i0<s0; i0++) {
627 for (i1=0; i1<s1; i1++) {
628 for (i2=0; i2<s2; i2++) {
629 for (i3=0; i3<s3; i3++) {
630 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
631 }
632 }
633 }
634 }
635 }
636 } else if (axis0==2) {
637 if (axis1==3) {
638 for (i0=0; i0<s0; i0++) {
639 for (i1=0; i1<s1; i1++) {
640 for (i2=0; i2<s2; i2++) {
641 for (i3=0; i3<s3; i3++) {
642 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
643 }
644 }
645 }
646 }
647 }
648 }
649
650 } else if ( inRank == 3) {
651 int s0=evShape[0];
652 int s1=evShape[1];
653 int s2=evShape[2];
654 int i0, i1, i2;
655 if (axis0==0) {
656 if (axis1==1) {
657 for (i0=0; i0<s0; i0++) {
658 for (i1=0; i1<s1; i1++) {
659 for (i2=0; i2<s2; i2++) {
660 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
661 }
662 }
663 }
664 } else if (axis1==2) {
665 for (i0=0; i0<s0; i0++) {
666 for (i1=0; i1<s1; i1++) {
667 for (i2=0; i2<s2; i2++) {
668 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
669 }
670 }
671 }
672 }
673 } else if (axis0==1) {
674 if (axis1==2) {
675 for (i0=0; i0<s0; i0++) {
676 for (i1=0; i1<s1; i1++) {
677 for (i2=0; i2<s2; i2++) {
678 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
679 }
680 }
681 }
682 }
683 }
684 } else if ( inRank == 2) {
685 int s0=evShape[0];
686 int s1=evShape[1];
687 int i0, i1;
688 if (axis0==0) {
689 if (axis1==1) {
690 for (i0=0; i0<s0; i0++) {
691 for (i1=0; i1<s1; i1++) {
692 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
693 }
694 }
695 }
696 }
697 } else {
698 throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
699 }
700 }
701
702 /**
703 \brief
704 solves a local eigenvalue problem
705
706 \param in - Input - matrix
707 \param inOffset - Input - offset into in
708 \param ev - Output - The eigenvalues
709 \param inOffset - Input - offset into ev
710 */
711 ESCRIPT_DLL_API
712 inline
713 void
714 eigenvalues(DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
715 DataTypes::ValueType::size_type inOffset,
716 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
717 DataTypes::ValueType::size_type evOffset)
718 {
719 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
720 double ev0,ev1,ev2;
721 int s=inShape[0];
722 if (s==1) {
723 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
724 eigenvalues1(in00,&ev0);
725 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
726
727 } else if (s==2) {
728 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
729 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
730 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
731 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
732 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
733 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
734 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
735
736 } else if (s==3) {
737 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
738 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
739 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
740 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
741 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
742 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
743 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
744 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
745 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
746 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
747 &ev0,&ev1,&ev2);
748 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
749 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
750 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
751
752 }
753 }
754
755 /**
756 \brief
757 solves a local eigenvalue problem
758
759 \param in - Input - matrix
760 \param inOffset - Input - offset into in
761 \param ev - Output - The eigenvalues
762 \param evOffset - Input - offset into ev
763 \param V - Output - The eigenvectors
764 \param VOffset - Input - offset into V
765 \param tol - Input - eigenvalues with relative difference tol are treated as equal
766 */
767 ESCRIPT_DLL_API
768 static
769 inline
770 void
771 eigenvalues_and_eigenvectors(const DataTypes::ValueType& in, const DataTypes::ShapeType& inShape,
772 DataTypes::ValueType::size_type inOffset,
773 DataTypes::ValueType& ev, const DataTypes::ShapeType& evShape,
774 DataTypes::ValueType::size_type evOffset,
775 DataTypes::ValueType& V, const DataTypes::ShapeType& VShape,
776 DataTypes::ValueType::size_type VOffset,
777 const double tol=1.e-13)
778 {
779 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
780 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
781 double ev0,ev1,ev2;
782 int s=inShape[0];
783 if (s==1) {
784 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
785 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
786 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
787 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
788 } else if (s==2) {
789 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
790 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
791 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
792 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
793 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
794 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
795 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
796 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
797 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
798 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
799 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
800 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
801 } else if (s==3) {
802 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
803 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
804 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
805 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
806 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
807 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
808 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
809 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
810 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
811 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
812 &ev0,&ev1,&ev2,
813 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
814 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
815 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
816 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
817 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
818 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
819 V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
820 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
821 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
822 V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
823 V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
824 V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
825 V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
826
827 }
828 }
829
830
831 /**
832 Inline function definitions.
833 */
834
835 // template <class UnaryFunction>
836 // inline
837 // void
838 // DataArrayView::unaryOp(UnaryFunction operation)
839 // {
840 // unaryOp(m_offset,operation);
841 // }
842
843
844 inline
845 bool
846 checkOffset(const DataTypes::ValueType& data,
847 const DataTypes::ShapeType& shape,
848 DataTypes::ValueType::size_type offset)
849 {
850 return (data.size() >= (offset+DataTypes::noValues(shape)));
851 }
852
853 template <class UnaryFunction>
854 inline
855 void
856 unaryOp(DataTypes::ValueType& data, const DataTypes::ShapeType& shape,
857 DataTypes::ValueType::size_type offset,
858 UnaryFunction operation)
859 {
860 EsysAssert((data.size()>0)&&checkOffset(data,shape,offset),
861 "Error - Couldn't perform unaryOp due to insufficient storage.");
862 DataTypes::ValueType::size_type nVals=DataTypes::noValues(shape);
863 for (DataTypes::ValueType::size_type i=0;i<nVals;i++) {
864 data[offset+i]=operation(data[offset+i]);
865 }
866 }
867
868 // template <class BinaryFunction>
869 // inline
870 // void
871 // binaryOp(const DataArrayView& right,
872 // BinaryFunction operation)
873 // {
874 // binaryOp(m_offset,right,right.getOffset(),operation);
875 // }
876
877
878 template <class BinaryFunction>
879 inline
880 void
881 binaryOp(DataTypes::ValueType& left,
882 const DataTypes::ShapeType& leftShape,
883 DataTypes::ValueType::size_type leftOffset,
884 const DataTypes::ValueType& right,
885 const DataTypes::ShapeType& rightShape,
886 DataTypes::ValueType::size_type rightOffset,
887 BinaryFunction operation)
888 {
889 EsysAssert(leftShape==rightShape,
890 "Error - Couldn't perform binaryOp due to shape mismatch,");
891 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape, leftOffset)),
892 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
893 EsysAssert(((right.size()>0)&&checkOffset(right,rightShape,rightOffset)),
894 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
895 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
896 left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
897 }
898 }
899
900 template <class BinaryFunction>
901 inline
902 void
903 binaryOp(DataTypes::ValueType& left,
904 const DataTypes::ShapeType& leftShape,
905 DataTypes::ValueType::size_type offset,
906 double right,
907 BinaryFunction operation)
908 {
909 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
910 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
911 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
912 left[offset+i]=operation(left[offset+i],right);
913 }
914 }
915
916 template <class BinaryFunction>
917 inline
918 double
919 reductionOp(const DataTypes::ValueType& left,
920 const DataTypes::ShapeType& leftShape,
921 DataTypes::ValueType::size_type offset,
922 BinaryFunction operation,
923 double initial_value)
924 {
925 EsysAssert(((left.size()>0)&&checkOffset(left,leftShape,offset)),
926 "Error - Couldn't perform reductionOp due to insufficient storage.");
927 double current_value=initial_value;
928 for (DataTypes::ValueType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
929 current_value=operation(current_value,left[offset+i]);
930 }
931 return current_value;
932 }
933
934
935
936 } // end namespace DataMath
937 } // end namespace escript
938 #endif
939

  ViewVC Help
Powered by ViewVC 1.1.26