17 #ifndef escript_DataMaths_20080822_H
18 #define escript_DataMaths_20080822_H
62 template <
class UnaryFunction>
66 UnaryFunction operation);
82 template <
class BinaryFunction>
90 BinaryFunction operation);
108 template <
class BinaryFunction>
114 BinaryFunction operation);
132 template <
class BinaryFunction>
137 BinaryFunction operation,
138 double initial_value);
207 for (i0=0; i0<s0; i0++) {
208 for (i1=0; i1<s1; i1++) {
219 for (i0=0; i0<s0; i0++) {
220 for (i1=0; i1<s1; i1++) {
221 for (i2=0; i2<s2; i2++) {
222 for (i3=0; i3<s3; i3++) {
223 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;
256 for (i0=0; i0<s0; i0++) {
257 for (i1=0; i1<s1; i1++) {
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;
309 for (i=0; i<s0; i++) {
314 if (axis_offset==0) {
318 for (i0=0; i0<s0; i0++) {
319 for (i2=0; i2<s2; i2++) {
324 else if (axis_offset==1) {
328 for (i0=0; i0<s0; i0++) {
329 for (i1=0; i1<s1; i1++) {
336 if (axis_offset==0) {
341 for (i0=0; i0<s0; i0++) {
342 for (i2=0; i2<s2; i2++) {
343 for (i3=0; i3<s3; i3++) {
349 else if (axis_offset==1) {
354 for (i0=0; i0<s0; i0++) {
355 for (i1=0; i1<s1; i1++) {
356 for (i3=0; i3<s3; i3++) {
362 else if (axis_offset==2) {
367 for (i0=0; i0<s0; i0++) {
368 for (i1=0; i1<s1; i1++) {
369 for (i2=0; i2<s2; i2++) {
408 if (axis_offset==1) {
409 for (i0=0; i0<s0; i0++) {
410 for (i1=0; i1<s1; i1++) {
411 for (i2=0; i2<s2; i2++) {
412 for (i3=0; i3<s3; i3++) {
419 else if (axis_offset==2) {
420 for (i0=0; i0<s0; i0++) {
421 for (i1=0; i1<s1; i1++) {
422 for (i2=0; i2<s2; i2++) {
423 for (i3=0; i3<s3; i3++) {
430 else if (axis_offset==3) {
431 for (i0=0; i0<s0; i0++) {
432 for (i1=0; i1<s1; i1++) {
433 for (i2=0; i2<s2; i2++) {
434 for (i3=0; i3<s3; i3++) {
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++) {
453 else if (inRank == 3) {
458 if (axis_offset==1) {
459 for (i0=0; i0<s0; i0++) {
460 for (i1=0; i1<s1; i1++) {
461 for (i2=0; i2<s2; i2++) {
467 else if (axis_offset==2) {
468 for (i0=0; i0<s0; i0++) {
469 for (i1=0; i1<s1; i1++) {
470 for (i2=0; i2<s2; i2++) {
478 for (i0=0; i0<s0; i0++) {
479 for (i1=0; i1<s1; i1++) {
480 for (i2=0; i2<s2; i2++) {
487 else if (inRank == 2) {
491 if (axis_offset==1) {
492 for (i0=0; i0<s0; i0++) {
493 for (i1=0; i1<s1; i1++) {
499 for (i0=0; i0<s0; i0++) {
500 for (i1=0; i1<s1; i1++) {
506 else if (inRank == 1) {
509 for (i0=0; i0<s0; i0++) {
513 else if (inRank == 0) {
514 ev[evOffset] = in[inOffset];
517 throw DataException(
"Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
555 for (i0=0; i0<s0; i0++) {
556 for (i1=0; i1<s1; i1++) {
557 for (i2=0; i2<s2; i2++) {
558 for (i3=0; i3<s3; i3++) {
564 }
else if (axis1==2) {
565 for (i0=0; i0<s0; i0++) {
566 for (i1=0; i1<s1; i1++) {
567 for (i2=0; i2<s2; i2++) {
568 for (i3=0; i3<s3; i3++) {
575 }
else if (axis1==3) {
576 for (i0=0; i0<s0; i0++) {
577 for (i1=0; i1<s1; i1++) {
578 for (i2=0; i2<s2; i2++) {
579 for (i3=0; i3<s3; i3++) {
586 }
else if (axis0==1) {
588 for (i0=0; i0<s0; i0++) {
589 for (i1=0; i1<s1; i1++) {
590 for (i2=0; i2<s2; i2++) {
591 for (i3=0; i3<s3; i3++) {
597 }
else if (axis1==3) {
598 for (i0=0; i0<s0; i0++) {
599 for (i1=0; i1<s1; i1++) {
600 for (i2=0; i2<s2; i2++) {
601 for (i3=0; i3<s3; i3++) {
608 }
else if (axis0==2) {
610 for (i0=0; i0<s0; i0++) {
611 for (i1=0; i1<s1; i1++) {
612 for (i2=0; i2<s2; i2++) {
613 for (i3=0; i3<s3; i3++) {
622 }
else if ( inRank == 3) {
629 for (i0=0; i0<s0; i0++) {
630 for (i1=0; i1<s1; i1++) {
631 for (i2=0; i2<s2; i2++) {
636 }
else if (axis1==2) {
637 for (i0=0; i0<s0; i0++) {
638 for (i1=0; i1<s1; i1++) {
639 for (i2=0; i2<s2; i2++) {
645 }
else if (axis0==1) {
647 for (i0=0; i0<s0; i0++) {
648 for (i1=0; i1<s1; i1++) {
649 for (i2=0; i2<s2; i2++) {
656 }
else if ( inRank == 2) {
662 for (i0=0; i0<s0; i0++) {
663 for (i1=0; i1<s1; i1++) {
670 throw DataException(
"Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
695 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
722 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
755 const double tol=1.e-13)
757 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
758 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
772 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
791 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
822 template <
class UnaryFunction>
827 UnaryFunction operation)
830 "Error - Couldn't perform unaryOp due to insufficient storage.");
833 data[offset+i]=operation(data[offset+i]);
838 template <
class BinaryFunction>
847 BinaryFunction operation)
850 "Error - Couldn't perform binaryOp due to shape mismatch,");
852 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
854 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
856 left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
860 template <
class BinaryFunction>
867 BinaryFunction operation)
870 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
872 left[offset+i]=operation(left[offset+i],right);
876 template <
class BinaryFunction>
882 BinaryFunction operation,
883 double initial_value)
886 "Error - Couldn't perform reductionOp due to insufficient storage.");
887 double current_value=initial_value;
889 current_value=operation(current_value,left[offset+i]);
891 return current_value;
938 for (
size_t z=inOffset;z<inOffset+count;++z)