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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1724 - (show annotations)
Mon Aug 25 05:38:57 2008 UTC (11 years, 2 months ago) by jfenwick
Original Path: branches/arrayview_from_1695_trunk/escript/src/DataMaths.h
File MIME type: text/plain
File size: 32145 byte(s)
Branch commit

Moved createShapeErrorMessage() into DataTypes.h
Modified functions in DataAlgorithm.h to use non-DataArrayView accesses.

Added getVector() to each of DataTagged, DataConstant, DataExpanded - This method returns 
the underlying DataVector by reference/constant reference. Note that this method does not 
exist in DataAbstract so (at the momement) in order to pull the data from something you 
need to know what you are looking at. (Lower level access is still possible via double* 
though).

DataTagged now has a getOffsetForTag method and a getDefaultOffset method.

DataMaths.h - A new file containing the reductionOps etc from DataArrayView (but without 
requiring DAV).
This file requires significant commenting improvements.


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

  ViewVC Help
Powered by ViewVC 1.1.26