19#ifndef __ESCRIPT_LOCALOPS_H__
20#define __ESCRIPT_LOCALOPS_H__
31#ifdef ESYS_USE_BOOST_ACOS
32#include <boost/math/complex/acos.hpp>
36# define M_PI 3.14159265358979323846
96 return std::max(std::abs(x),std::abs(y));
125 return std::isnan(d);
136 return std::isnan( real(d) ) || std::isnan( imag(d) );
188 const T trA=(A00+A11)/2.;
189 const T A_00=A00-trA;
190 const T A_11=A11-trA;
191 const T s=sqrt(A01*A01-A_00*A_11);
229 const DataTypes::real_t q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
238 *ev2=trA+2.*sq_p*cos(alpha_3);
239 *ev1=trA-2.*sq_p*cos(alpha_3+
M_PI/3.);
240 *ev0=trA-2.*sq_p*cos(alpha_3-
M_PI/3.);
279 if (absA00>m || absA01>m) {
321 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
322 *V0=-(A10*TEMP0+A20*TEMP1);
355 if (
fabs((*ev0)-(*ev1))<tol*max_ev) {
373 }
else if (TEMP0>0.) {
404 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
409 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
415 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
419 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
469 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
483 }
else if (A00>TEMP_EV1) {
516 max_ev=max_ev>absev2 ? max_ev : absev2;
520 if (max_d<=tol*max_ev) {
534 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
535 }
else if (absA02<m) {
536 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
538 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
544 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
545 }
else if (absA02<m) {
546 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
548 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
555 *V01=(*V10)*(*V22)-(*V12)*(*V20);
556 *V11=(*V20)*(*V02)-(*V00)*(*V22);
557 *V21=(*V00)*(*V12)-(*V02)*(*V10);
565template <
class LEFT,
class RIGHT,
class RES>
570 for (
int i=0; i<SL; i++) {
571 for (
int j=0; j<SR; j++) {
573 for (
int l=0; l<SM; l++) {
574 sum += A[i+SL*l] * B[l+SM*j];
581 for (
int i=0; i<SL; i++) {
582 for (
int j=0; j<SR; j++) {
584 for (
int l=0; l<SM; l++) {
585 sum += A[i*SM+l] * B[l+SM*j];
592 for (
int i=0; i<SL; i++) {
593 for (
int j=0; j<SR; j++) {
595 for (
int l=0; l<SM; l++) {
596 sum += A[i+SL*l] * B[l*SR+j];
635#ifdef ESYS_USE_BOOST_ACOS
636 return boost::math::acos(x);
691 for (
int i = 0; i < size; ++i) {
692 argRes[i] = std::real(arg1[i]);
696 for (
int i = 0; i < size; ++i) {
697 argRes[i] = std::imag(arg1[i]);
701 for (
size_t i = 0; i < size; ++i) {
702 argRes[i] = (
fabs(arg1[i])<=tol);
706 for (
size_t i = 0; i < size; ++i) {
707 argRes[i] = (
fabs(arg1[i])>tol);
711 for (
size_t i = 0; i < size; ++i) {
712 argRes[i] =
abs_f(arg1[i]);
716 for (
size_t i = 0; i < size; ++i) {
717 argRes[i] = std::arg(arg1[i]);
721 std::ostringstream oss;
722 oss <<
"Unsupported unary operation=";
726 oss <<
" (Was expecting an operation with real results)";
733template <
typename T1,
typename T2>
750 for (
size_t i = 0; i < size; ++i) {
759template <
class T1,
typename T2>
769 for (
size_t i = 0; i < size; ++i) {
770 argRes[i] = -arg1[i];
774 for (
size_t i = 0; i < size; ++i) {
775 argRes[i] = sin(arg1[i]);
779 for (
size_t i = 0; i < size; ++i) {
780 argRes[i] = cos(arg1[i]);
784 for (
size_t i = 0; i < size; ++i) {
785 argRes[i] = tan(arg1[i]);
789 for (
size_t i = 0; i < size; ++i) {
790 argRes[i] = asin(arg1[i]);
794 for (
size_t i = 0; i < size; ++i) {
799 for (
size_t i = 0; i < size; ++i) {
800 argRes[i] = atan(arg1[i]);
804 for (
size_t i = 0; i < size; ++i) {
805 argRes[i] = std::abs(arg1[i]);
809 for (
size_t i = 0; i < size; ++i) {
810 argRes[i] = sinh(arg1[i]);
814 for (
size_t i = 0; i < size; ++i) {
815 argRes[i] = cosh(arg1[i]);
819 for (
size_t i = 0; i < size; ++i) {
820 argRes[i] = tanh(arg1[i]);
824 for (
size_t i = 0; i < size; ++i) {
829 for (
size_t i = 0; i < size; ++i) {
830 argRes[i] = asinh(arg1[i]);
834 for (
size_t i = 0; i < size; ++i) {
835 argRes[i] = acosh(arg1[i]);
839 for (
size_t i = 0; i < size; ++i) {
840 argRes[i] = atanh(arg1[i]);
844 for (
size_t i = 0; i < size; ++i) {
845 argRes[i] = log10(arg1[i]);
849 for (
size_t i = 0; i < size; ++i) {
850 argRes[i] = log(arg1[i]);
854 for (
size_t i = 0; i < size; ++i) {
859 for (
size_t i = 0; i < size; ++i) {
860 argRes[i] = exp(arg1[i]);
864 for (
size_t i = 0; i < size; ++i) {
865 argRes[i] = sqrt(arg1[i]);
869 for (
size_t i = 0; i < size; ++i) {
874 for (
size_t i = 0; i < size; ++i) {
879 for (
size_t i = 0; i < size; ++i) {
884 for (
size_t i = 0; i < size; ++i) {
889 for (
size_t i = 0; i < size; ++i) {
890 argRes[i] = conjugate<T2,T1>(arg1[i]);
894 for (
size_t i = 0; i < size; ++i) {
895 argRes[i] = 1.0/arg1[i];
899 for (
size_t i = 0; i < size; ++i) {
900 argRes[i] =
fabs(arg1[i])<=tol;
904 for (
size_t i = 0; i < size; ++i) {
905 argRes[i] =
fabs(arg1[i])>tol;
910 std::ostringstream oss;
911 oss <<
"Unsupported unary operation=";
#define M_PI
Definition ArrayOps.h:36
Definition DataException.h:28
std::complex< real_t > cplx_t
complex data type
Definition DataTypes.h:55
double real_t
type of all real-valued scalars in escript
Definition DataTypes.h:52
Definition AbstractContinuousDomain.cpp:23
void eigenvalues1(const DataTypes::real_t A00, DataTypes::real_t *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition ArrayOps.h:164
DataTypes::real_t calc_gezero(const DataTypes::real_t &x)
Definition ArrayOps.h:654
void vectorInKernel2(const DataTypes::real_t A00, const DataTypes::real_t A10, const DataTypes::real_t A01, const DataTypes::real_t A11, DataTypes::real_t *V0, DataTypes::real_t *V1)
returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension i...
Definition ArrayOps.h:271
DataTypes::real_t calc_ltzero(const DataTypes::real_t &x)
Definition ArrayOps.h:658
void eigenvalues3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2)
solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
Definition ArrayOps.h:210
void eigenvalues_and_eigenvectors3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V20, DataTypes::real_t *V01, DataTypes::real_t *V11, DataTypes::real_t *V21, DataTypes::real_t *V02, DataTypes::real_t *V12, DataTypes::real_t *V22, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition ArrayOps.h:454
DataTypes::real_t calc_acos(DataTypes::real_t x)
Definition ArrayOps.h:627
DataTypes::real_t calc_gtzero(const DataTypes::real_t &x)
Definition ArrayOps.h:650
void tensor_unary_promote(const size_t size, const DataTypes::real_t *arg1, DataTypes::cplx_t *argRes)
Definition ArrayOps.h:746
T1 conjugate(const T2 i)
Definition ArrayOps.h:734
DataTypes::real_t calc_lezero(const DataTypes::real_t &x)
Definition ArrayOps.h:661
void eigenvalues_and_eigenvectors2(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A11, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V01, DataTypes::real_t *V11, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition ArrayOps.h:345
void transpose(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis_offset)
Transpose each data point of this Data object around the given axis.
Definition DataVectorOps.h:343
DataTypes::real_t makeNaN()
returns a NaN.
Definition ArrayOps.h:145
void normalizeVector3(DataTypes::real_t *V0, DataTypes::real_t *V1, DataTypes::real_t *V2)
nomalizes a 3-d vector such that length is one and first non-zero component is positive.
Definition ArrayOps.h:400
void eigenvalues2(const T A00, const T A01, const T A11, T *ev0, T *ev1)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
Definition ArrayOps.h:186
void tensor_unary_array_operation_real(const size_t size, const T *arg1, DataTypes::real_t *argRes, escript::ES_optype operation, DataTypes::real_t tol=0)
Definition ArrayOps.h:682
void eigenvalues_and_eigenvectors1(const DataTypes::real_t A00, DataTypes::real_t *ev0, DataTypes::real_t *V00, const DataTypes::real_t tol)
solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
Definition ArrayOps.h:253
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition ArrayOps.h:120
bool supports_cplx(escript::ES_optype operation)
Definition ArrayOps.cpp:26
DataTypes::real_t calc_erf(DataTypes::real_t x)
Definition ArrayOps.h:605
const std::string & opToString(ES_optype op)
Definition ES_optype.cpp:89
void vectorInKernel3__nonZeroA00(const DataTypes::real_t A00, const DataTypes::real_t A10, const DataTypes::real_t A20, const DataTypes::real_t A01, const DataTypes::real_t A11, const DataTypes::real_t A21, const DataTypes::real_t A02, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *V0, DataTypes::real_t *V1, DataTypes::real_t *V2)
returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] assuming that ...
Definition ArrayOps.h:311
ES_optype
Definition ES_optype.h:30
@ ABS
Definition ES_optype.h:54
@ RECIP
Definition ES_optype.h:59
@ NEZ
Definition ES_optype.h:64
@ LZ
Definition ES_optype.h:61
@ LEZ
Definition ES_optype.h:63
@ LOG10
Definition ES_optype.h:51
@ SINH
Definition ES_optype.h:44
@ ATANH
Definition ES_optype.h:50
@ REAL
Definition ES_optype.h:77
@ GZ
Definition ES_optype.h:60
@ EXP
Definition ES_optype.h:57
@ TANH
Definition ES_optype.h:46
@ COSH
Definition ES_optype.h:45
@ IMAG
Definition ES_optype.h:78
@ ERF
Definition ES_optype.h:47
@ SIN
Definition ES_optype.h:38
@ ASINH
Definition ES_optype.h:48
@ ACOSH
Definition ES_optype.h:49
@ LOG
Definition ES_optype.h:52
@ EZ
Definition ES_optype.h:65
@ SIGN
Definition ES_optype.h:53
@ TAN
Definition ES_optype.h:40
@ PHS
Definition ES_optype.h:84
@ ASIN
Definition ES_optype.h:41
@ GEZ
Definition ES_optype.h:62
@ SQRT
Definition ES_optype.h:58
@ CONJ
Definition ES_optype.h:79
@ ACOS
Definition ES_optype.h:42
@ ATAN
Definition ES_optype.h:43
@ COS
Definition ES_optype.h:39
@ NEG
Definition ES_optype.h:55
DataTypes::real_t calc_sign(DataTypes::real_t x)
Definition ArrayOps.h:616
DataTypes::real_t abs_f(T i)
Definition ArrayOps.h:665
DataTypes::real_t fsign(DataTypes::real_t x)
Definition ArrayOps.h:106
bool always_real(escript::ES_optype operation)
Definition ArrayOps.cpp:66
void matrix_matrix_product(const int SL, const int SM, const int SR, const LEFT *A, const RIGHT *B, RES *C, int transpose)
Definition ArrayOps.h:567
escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
Definition ArrayOps.h:643
void tensor_unary_array_operation(const size_t size, const T1 *arg1, T2 *argRes, escript::ES_optype operation, DataTypes::real_t tol=0)
Definition ArrayOps.h:760
Return the absolute maximum value of the two given values.
Definition ArrayOps.h:93
DataTypes::real_t operator()(T x, T y) const
Definition ArrayOps.h:94
T second_argument_type
Definition ArrayOps.h:99
T first_argument_type
Definition ArrayOps.h:98
DataTypes::real_t result_type
Definition ArrayOps.h:100
Return the maximum value of the two given values.
Definition ArrayOps.h:62
DataTypes::real_t second_argument_type
Definition ArrayOps.h:68
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition ArrayOps.h:63
DataTypes::real_t result_type
Definition ArrayOps.h:69
DataTypes::real_t first_argument_type
Definition ArrayOps.h:67
Return the minimum value of the two given values.
Definition ArrayOps.h:77
DataTypes::real_t result_type
Definition ArrayOps.h:84
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition ArrayOps.h:78
DataTypes::real_t first_argument_type
Definition ArrayOps.h:82
DataTypes::real_t second_argument_type
Definition ArrayOps.h:83