Go to the documentation of this file.
19 #ifndef __ESCRIPT_LOCALOPS_H__
20 #define __ESCRIPT_LOCALOPS_H__
28 #ifdef ESYS_USE_BOOST_ACOS
29 #include <boost/math/complex/acos.hpp>
33 # define M_PI 3.14159265358979323846
93 return std::max(std::abs(x),std::abs(y));
122 return std::isnan(d);
133 return std::isnan( real(d) ) || std::isnan( imag(d) );
185 const T trA=(A00+A11)/2.;
186 const T A_00=A00-trA;
187 const T A_11=A11-trA;
188 const T s=sqrt(A01*A01-A_00*A_11);
226 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);
235 *ev2=trA+2.*sq_p*cos(alpha_3);
236 *ev1=trA-2.*sq_p*cos(alpha_3+
M_PI/3.);
237 *ev0=trA-2.*sq_p*cos(alpha_3-
M_PI/3.);
276 if (absA00>m || absA01>m) {
318 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
319 *V0=-(A10*TEMP0+A20*TEMP1);
352 if (
fabs((*ev0)-(*ev1))<tol*max_ev) {
370 }
else if (TEMP0>0.) {
401 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
406 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
412 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
416 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
466 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
480 }
else if (A00>TEMP_EV1) {
513 max_ev=max_ev>absev2 ? max_ev : absev2;
517 if (max_d<=tol*max_ev) {
531 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
532 }
else if (absA02<m) {
533 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
535 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
541 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
542 }
else if (absA02<m) {
543 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
545 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
552 *V01=(*V10)*(*V22)-(*V12)*(*V20);
553 *V11=(*V20)*(*V02)-(*V00)*(*V22);
554 *V21=(*V00)*(*V12)-(*V02)*(*V10);
562 template <
class LEFT,
class RIGHT,
class RES>
567 for (
int i=0; i<SL; i++) {
568 for (
int j=0; j<SR; j++) {
570 for (
int l=0; l<SM; l++) {
571 sum += A[i+SL*l] * B[l+SM*j];
578 for (
int i=0; i<SL; i++) {
579 for (
int j=0; j<SR; j++) {
581 for (
int l=0; l<SM; l++) {
582 sum += A[i*SM+l] * B[l+SM*j];
589 for (
int i=0; i<SL; i++) {
590 for (
int j=0; j<SR; j++) {
592 for (
int l=0; l<SM; l++) {
593 sum += A[i+SL*l] * B[l*SR+j];
601 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
637 #ifdef ESYS_USE_BOOST_ACOS
638 return boost::math::acos(x);
666 template <
typename IN>
693 for (
int i = 0; i < size; ++i) {
694 argRes[i] = std::real(arg1[i]);
698 for (
int i = 0; i < size; ++i) {
699 argRes[i] = std::imag(arg1[i]);
703 for (
size_t i = 0; i < size; ++i) {
704 argRes[i] = (
fabs(arg1[i])<=tol);
708 for (
size_t i = 0; i < size; ++i) {
709 argRes[i] = (
fabs(arg1[i])>tol);
713 for (
size_t i = 0; i < size; ++i) {
714 argRes[i] =
abs_f(arg1[i]);
718 for (
size_t i = 0; i < size; ++i) {
719 argRes[i] = std::arg(arg1[i]);
723 std::ostringstream oss;
724 oss <<
"Unsupported unary operation=";
728 oss <<
" (Was expecting an operation with real results)";
735 template <
typename OUT,
typename IN>
752 for (
size_t i = 0; i < size; ++i) {
761 template <
class IN,
typename OUT>
771 for (
size_t i = 0; i < size; ++i) {
772 argRes[i] = -arg1[i];
776 for (
size_t i = 0; i < size; ++i) {
777 argRes[i] = sin(arg1[i]);
781 for (
size_t i = 0; i < size; ++i) {
782 argRes[i] = cos(arg1[i]);
786 for (
size_t i = 0; i < size; ++i) {
787 argRes[i] = tan(arg1[i]);
791 for (
size_t i = 0; i < size; ++i) {
792 argRes[i] = asin(arg1[i]);
796 for (
size_t i = 0; i < size; ++i) {
801 for (
size_t i = 0; i < size; ++i) {
802 argRes[i] = atan(arg1[i]);
806 for (
size_t i = 0; i < size; ++i) {
807 argRes[i] = std::abs(arg1[i]);
811 for (
size_t i = 0; i < size; ++i) {
812 argRes[i] = sinh(arg1[i]);
816 for (
size_t i = 0; i < size; ++i) {
817 argRes[i] = cosh(arg1[i]);
821 for (
size_t i = 0; i < size; ++i) {
822 argRes[i] = tanh(arg1[i]);
826 for (
size_t i = 0; i < size; ++i) {
831 for (
size_t i = 0; i < size; ++i) {
832 argRes[i] = asinh(arg1[i]);
836 for (
size_t i = 0; i < size; ++i) {
837 argRes[i] = acosh(arg1[i]);
841 for (
size_t i = 0; i < size; ++i) {
842 argRes[i] = atanh(arg1[i]);
846 for (
size_t i = 0; i < size; ++i) {
847 argRes[i] = log10(arg1[i]);
851 for (
size_t i = 0; i < size; ++i) {
852 argRes[i] = log(arg1[i]);
856 for (
size_t i = 0; i < size; ++i) {
861 for (
size_t i = 0; i < size; ++i) {
862 argRes[i] = exp(arg1[i]);
866 for (
size_t i = 0; i < size; ++i) {
867 argRes[i] = sqrt(arg1[i]);
871 for (
size_t i = 0; i < size; ++i) {
876 for (
size_t i = 0; i < size; ++i) {
881 for (
size_t i = 0; i < size; ++i) {
886 for (
size_t i = 0; i < size; ++i) {
891 for (
size_t i = 0; i < size; ++i) {
892 argRes[i] = conjugate<OUT,IN>(arg1[i]);
896 for (
size_t i = 0; i < size; ++i) {
897 argRes[i] = 1.0/arg1[i];
901 for (
size_t i = 0; i < size; ++i) {
902 argRes[i] =
fabs(arg1[i])<=tol;
906 for (
size_t i = 0; i < size; ++i) {
907 argRes[i] =
fabs(arg1[i])>tol;
912 std::ostringstream oss;
913 oss <<
"Unsupported unary operation=";
927 #endif // __ESCRIPT_LOCALOPS_H__
void tensor_unary_array_operation_real(const size_t size, const IN *arg1, DataTypes::real_t *argRes, escript::ES_optype operation, DataTypes::real_t tol=0)
Definition: ArrayOps.h:683
DataTypes::real_t calc_acos(DataTypes::real_t x)
Definition: ArrayOps.h:628
Return the minimum value of the two given values.
Definition: ArrayOps.h:72
DataTypes::real_t operator()(T x, T y) const
Definition: ArrayOps.h:90
Definition: ES_optype.h:76
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:76
OUT conjugate(const IN i)
Definition: ArrayOps.h:735
Definition: ES_optype.h:102
Definition: ES_optype.h:84
Definition: ES_optype.h:68
void eigenvalues1(const DataTypes::real_t A00, DataTypes::real_t *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: ArrayOps.h:160
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition: ArrayOps.h:116
DataTypes::real_t calc_gezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:655
Definition: ES_optype.h:81
Definition: ES_optype.h:101
DataTypes::real_t calc_erf(DataTypes::real_t x)
Definition: ArrayOps.h:604
Definition: ES_optype.h:80
T first_argument_type
Definition: ArrayOps.h:94
DataTypes::real_t calc_ltzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:659
Definition: ES_optype.h:67
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:342
const std::string & opToString(ES_optype op)
Definition: ES_optype.cpp:88
Definition: ES_optype.h:73
Definition: ES_optype.h:75
DataTypes::real_t abs_f(IN i)
Definition: ArrayOps.h:666
DataTypes::real_t fsign(DataTypes::real_t x)
Definition: ArrayOps.h:102
DataTypes::real_t calc_lezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:662
Definition: ES_optype.h:65
escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
Definition: ArrayOps.h:644
Return the absolute maximum value of the two given values.
Definition: ArrayOps.h:88
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:78
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:64
Definition: ES_optype.h:69
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:450
Definition: ES_optype.h:107
Definition: ES_optype.h:70
Definition: ES_optype.h:78
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:182
Definition: ES_optype.h:83
Definition: DataException.h:39
DataTypes::real_t result_type
Definition: ArrayOps.h:65
Definition: ES_optype.h:77
DataTypes::real_t result_type
Definition: ArrayOps.h:96
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:63
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:563
Definition: ES_optype.h:100
Definition: ES_optype.h:72
Definition: ES_optype.h:85
Definition: ES_optype.h:86
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:307
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:79
T second_argument_type
Definition: ArrayOps.h:95
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:74
Definition: ES_optype.h:82
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:396
ES_optype
Definition: ES_optype.h:39
void scale(dim_t N, double *x, double a)
x = a*x
Definition: PasoUtil.h:119
#define M_PI
Definition: ArrayOps.h:32
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:249
Definition: AbstractContinuousDomain.cpp:23
Definition: ES_optype.h:74
bool always_real(escript::ES_optype operation)
Definition: ArrayOps.cpp:65
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:341
bool supports_cplx(escript::ES_optype operation)
Definition: ArrayOps.cpp:25
Definition: ES_optype.h:61
Definition: ES_optype.h:88
DataTypes::real_t makeNaN()
returns a NaN.
Definition: ArrayOps.h:141
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:59
Definition: ES_optype.h:63
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:79
DataTypes::real_t calc_gtzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:651
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:267
void tensor_unary_array_operation(const size_t size, const IN *arg1, OUT *argRes, escript::ES_optype operation, DataTypes::real_t tol=0)
Definition: ArrayOps.h:761
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:206
Definition: ES_optype.h:64
Definition: ES_optype.h:62
Definition: ES_optype.h:71
DataTypes::real_t result_type
Definition: ArrayOps.h:80
void tensor_unary_promote(const size_t size, const DataTypes::real_t *arg1, DataTypes::cplx_t *argRes)
Definition: ArrayOps.h:747
Definition: ES_optype.h:66
Definition: ES_optype.h:87
DataTypes::real_t calc_sign(DataTypes::real_t x)
Definition: ArrayOps.h:617