18 #ifndef escript_DataMaths_20080822_H 19 #define escript_DataMaths_20080822_H 63 template <
class UnaryFunction>
67 UnaryFunction operation);
83 template <
class BinaryFunction>
91 BinaryFunction operation);
109 template <
class BinaryFunction>
115 BinaryFunction operation);
133 template <
class BinaryFunction>
138 BinaryFunction operation,
139 double initial_value);
208 for (i0=0; i0<s0; i0++) {
209 for (i1=0; i1<s1; i1++) {
220 for (i0=0; i0<s0; i0++) {
221 for (i1=0; i1<s1; i1++) {
222 for (i2=0; i2<s2; i2++) {
223 for (i3=0; i3<s3; i3++) {
224 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;
257 for (i0=0; i0<s0; i0++) {
258 for (i1=0; i1<s1; i1++) {
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;
310 for (i=0; i<s0; i++) {
315 if (axis_offset==0) {
319 for (i0=0; i0<s0; i0++) {
320 for (i2=0; i2<s2; i2++) {
325 else if (axis_offset==1) {
329 for (i0=0; i0<s0; i0++) {
330 for (i1=0; i1<s1; i1++) {
337 if (axis_offset==0) {
342 for (i0=0; i0<s0; i0++) {
343 for (i2=0; i2<s2; i2++) {
344 for (i3=0; i3<s3; i3++) {
350 else if (axis_offset==1) {
355 for (i0=0; i0<s0; i0++) {
356 for (i1=0; i1<s1; i1++) {
357 for (i3=0; i3<s3; i3++) {
363 else if (axis_offset==2) {
368 for (i0=0; i0<s0; i0++) {
369 for (i1=0; i1<s1; i1++) {
370 for (i2=0; i2<s2; i2++) {
409 if (axis_offset==1) {
410 for (i0=0; i0<s0; i0++) {
411 for (i1=0; i1<s1; i1++) {
412 for (i2=0; i2<s2; i2++) {
413 for (i3=0; i3<s3; i3++) {
420 else if (axis_offset==2) {
421 for (i0=0; i0<s0; i0++) {
422 for (i1=0; i1<s1; i1++) {
423 for (i2=0; i2<s2; i2++) {
424 for (i3=0; i3<s3; i3++) {
431 else if (axis_offset==3) {
432 for (i0=0; i0<s0; i0++) {
433 for (i1=0; i1<s1; i1++) {
434 for (i2=0; i2<s2; i2++) {
435 for (i3=0; i3<s3; i3++) {
443 for (i0=0; i0<s0; i0++) {
444 for (i1=0; i1<s1; i1++) {
445 for (i2=0; i2<s2; i2++) {
446 for (i3=0; i3<s3; i3++) {
454 else if (inRank == 3) {
459 if (axis_offset==1) {
460 for (i0=0; i0<s0; i0++) {
461 for (i1=0; i1<s1; i1++) {
462 for (i2=0; i2<s2; i2++) {
468 else if (axis_offset==2) {
469 for (i0=0; i0<s0; i0++) {
470 for (i1=0; i1<s1; i1++) {
471 for (i2=0; i2<s2; i2++) {
479 for (i0=0; i0<s0; i0++) {
480 for (i1=0; i1<s1; i1++) {
481 for (i2=0; i2<s2; i2++) {
488 else if (inRank == 2) {
492 if (axis_offset==1) {
493 for (i0=0; i0<s0; i0++) {
494 for (i1=0; i1<s1; i1++) {
500 for (i0=0; i0<s0; i0++) {
501 for (i1=0; i1<s1; i1++) {
507 else if (inRank == 1) {
510 for (i0=0; i0<s0; i0++) {
514 else if (inRank == 0) {
515 ev[evOffset] = in[inOffset];
518 throw DataException(
"Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
556 for (i0=0; i0<s0; i0++) {
557 for (i1=0; i1<s1; i1++) {
558 for (i2=0; i2<s2; i2++) {
559 for (i3=0; i3<s3; i3++) {
565 }
else if (axis1==2) {
566 for (i0=0; i0<s0; i0++) {
567 for (i1=0; i1<s1; i1++) {
568 for (i2=0; i2<s2; i2++) {
569 for (i3=0; i3<s3; i3++) {
576 }
else if (axis1==3) {
577 for (i0=0; i0<s0; i0++) {
578 for (i1=0; i1<s1; i1++) {
579 for (i2=0; i2<s2; i2++) {
580 for (i3=0; i3<s3; i3++) {
587 }
else if (axis0==1) {
589 for (i0=0; i0<s0; i0++) {
590 for (i1=0; i1<s1; i1++) {
591 for (i2=0; i2<s2; i2++) {
592 for (i3=0; i3<s3; i3++) {
598 }
else if (axis1==3) {
599 for (i0=0; i0<s0; i0++) {
600 for (i1=0; i1<s1; i1++) {
601 for (i2=0; i2<s2; i2++) {
602 for (i3=0; i3<s3; i3++) {
609 }
else if (axis0==2) {
611 for (i0=0; i0<s0; i0++) {
612 for (i1=0; i1<s1; i1++) {
613 for (i2=0; i2<s2; i2++) {
614 for (i3=0; i3<s3; i3++) {
623 }
else if ( inRank == 3) {
630 for (i0=0; i0<s0; i0++) {
631 for (i1=0; i1<s1; i1++) {
632 for (i2=0; i2<s2; i2++) {
637 }
else if (axis1==2) {
638 for (i0=0; i0<s0; i0++) {
639 for (i1=0; i1<s1; i1++) {
640 for (i2=0; i2<s2; i2++) {
646 }
else if (axis0==1) {
648 for (i0=0; i0<s0; i0++) {
649 for (i1=0; i1<s1; i1++) {
650 for (i2=0; i2<s2; i2++) {
657 }
else if ( inRank == 2) {
663 for (i0=0; i0<s0; i0++) {
664 for (i1=0; i1<s1; i1++) {
671 throw DataException(
"Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
696 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
723 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
756 const double tol=1.e-13)
758 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
759 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
773 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
792 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
823 template <
class UnaryFunction>
828 UnaryFunction operation)
831 "Error - Couldn't perform unaryOp due to insufficient storage.");
834 data[offset+i]=operation(data[offset+i]);
839 template <
class BinaryFunction>
848 BinaryFunction operation)
851 "Error - Couldn't perform binaryOp due to shape mismatch,");
853 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
855 "Error - Couldn't perform binaryOp due to insufficient storage in right object.");
857 left[leftOffset+i]=operation(left[leftOffset+i],right[rightOffset+i]);
861 template <
class BinaryFunction>
868 BinaryFunction operation)
871 "Error - Couldn't perform binaryOp due to insufficient storage in left object.");
873 left[offset+i]=operation(left[offset+i],right);
877 template <
class BinaryFunction>
883 BinaryFunction operation,
884 double initial_value)
887 "Error - Couldn't perform reductionOp due to insufficient storage.");
888 double current_value=initial_value;
890 current_value=operation(current_value,left[offset+i]);
892 return current_value;
939 for (
size_t z=inOffset;z<inOffset+count;++z)
DataVector implements an arbitrarily long vector of data values. DataVector is the underlying data co...
Definition: DataVector.h:44
DataTypes::ShapeType determineResultShape(const DataTypes::ShapeType &left, const DataTypes::ShapeType &right)
Determine the shape of the result array for a matrix multiplication of the given views.
Definition: DataMaths.cpp:173
void matMult(const DataTypes::ValueType &left, const DataTypes::ShapeType &leftShape, DataTypes::ValueType::size_type leftOffset, const DataTypes::ValueType &right, const DataTypes::ShapeType &rightShape, DataTypes::ValueType::size_type rightOffset, DataTypes::ValueType &result, const DataTypes::ShapeType &resultShape)
Perform a matrix multiply of the given views.
Definition: DataMaths.cpp:42
escript::DataVector ValueType
Vector to store underlying data.
Definition: DataTypes.h:37
void matrixInverseError(int err)
throws an appropriate exception based on failure of matrix_inverse.
Definition: DataMaths.cpp:189
Definition: AbstractContinuousDomain.cpp:24
void transpose(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis_offset)
Transpose each data point of this Data object around the given axis.
Definition: DataMaths.h:394
void binaryOp(DataTypes::ValueType &left, const DataTypes::ShapeType &leftShape, DataTypes::ValueType::size_type leftOffset, const DataTypes::ValueType &right, const DataTypes::ShapeType &rightShape, DataTypes::ValueType::size_type rightOffset, BinaryFunction operation)
Perform the binary operation on the data points specified by the given offsets in the "left" and "rig...
Definition: DataMaths.h:842
size_type size() const
Return the number of elements in this DataVector.
Definition: DataVector.h:215
Definition: LapackInverseHelper.h:26
bool vectorHasNaN(const DataTypes::ValueType &in, DataTypes::ValueType::size_type inOffset, size_t count)
returns true if the vector contains NaN
Definition: DataMaths.h:937
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:38
void unaryOp(DataTypes::ValueType &data, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset, UnaryFunction operation)
Perform the unary operation on the data point specified by the given offset. Applies the specified op...
Definition: DataMaths.h:826
void swapaxes(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis0, int axis1)
swaps the components axis0 and axis1.
Definition: DataMaths.h:538
void eigenvalues(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
solves a local eigenvalue problem
Definition: DataMaths.h:689
void symmetric(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
Definition: DataMaths.h:197
void eigenvalues2(const double A00, const double A01, const double A11, double *ev0, double *ev1)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:94
double reductionOp(const DataTypes::ValueType &left, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset, BinaryFunction operation, double initial_value)
Perform the given data point reduction operation on the data point specified by the given offset into...
Definition: DataMaths.h:880
Describes binary operations performed on double*.
bool nancheck(double d)
acts as a wrapper to isnan.
Definition: LocalOps.h:41
DataTypes::ValueType::size_type getRelIndex(const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:183
#define EsysAssert(AssertTest, AssertMessage)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false...
Definition: EsysAssert.h:96
void eigenvalues_and_eigenvectors(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, DataTypes::ValueType &V, const DataTypes::ShapeType &VShape, DataTypes::ValueType::size_type VOffset, const double tol=1.e-13)
solves a local eigenvalue problem
Definition: DataMaths.h:750
int matrix_inverse(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &out, const DataTypes::ShapeType &outShape, DataTypes::ValueType::size_type outOffset, int count, LapackInverseHelper &helper)
computes the inverses of square (up to 3x3) matricies
Definition: DataMaths.cpp:210
void eigenvalues_and_eigenvectors1(const double A00, double *ev0, double *V00, const double tol)
solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:161
#define V(_K_, _I_)
Definition: ShapeFunctions.cpp:120
void eigenvalues1(const double A00, double *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: LocalOps.h:78
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:94
DataException exception class.
Definition: DataException.h:35
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:54
void eigenvalues_and_eigenvectors2(const double A00, const double A01, const double A11, double *ev0, double *ev1, double *V00, double *V10, double *V01, double *V11, const double tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: LocalOps.h:253
void nonsymmetric(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset)
computes a nonsymmetric matrix from your square matrix A: (A - transpose(A)) / 2
Definition: DataMaths.h:246
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:167
void eigenvalues3(const double A00, const double A01, const double A02, const double A11, const double A12, const double A22, double *ev0, double *ev1, double *ev2)
solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
Definition: LocalOps.h:118
void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02, const double A11, const double A12, const double A22, double *ev0, double *ev1, double *ev2, double *V00, double *V10, double *V20, double *V01, double *V11, double *V21, double *V02, double *V12, double *V22, const double tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: LocalOps.h:362
long size_type
Definition: DataVector.h:60
void trace(const DataTypes::ValueType &in, const DataTypes::ShapeType &inShape, DataTypes::ValueType::size_type inOffset, DataTypes::ValueType &ev, const DataTypes::ShapeType &evShape, DataTypes::ValueType::size_type evOffset, int axis_offset)
computes the trace of a matrix
Definition: DataMaths.h:295
bool checkOffset(const DataTypes::ValueType &data, const DataTypes::ShapeType &shape, DataTypes::ValueType::size_type offset)
Definition: DataMaths.h:816