18 #if !defined escript_LocalOps_H 19 #define escript_LocalOps_H 22 # define M_PI 3.14159265358979323846 94 void eigenvalues2(
const double A00,
const double A01,
const double A11,
95 double* ev0,
double* ev1) {
96 const register double trA=(A00+A11)/2.;
97 const register double A_00=A00-trA;
98 const register double A_11=A11-trA;
99 const register double s=sqrt(A01*A01-A_00*A_11);
118 void eigenvalues3(
const double A00,
const double A01,
const double A02,
119 const double A11,
const double A12,
121 double* ev0,
double* ev1,
double* ev2) {
123 const register double trA=(A00+A11+A22)/3.;
124 const register double A_00=A00-trA;
125 const register double A_11=A11-trA;
126 const register double A_22=A22-trA;
127 const register double A01_2=A01*A01;
128 const register double A02_2=A02*A02;
129 const register double A12_2=A12*A12;
130 const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
137 const register double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
138 const register double sq_p=sqrt(p/3.);
139 register double z=-q/(2*pow(sq_p,3));
145 const register double alpha_3=acos(z)/3.;
146 *ev2=trA+2.*sq_p*cos(alpha_3);
147 *ev1=trA-2.*sq_p*cos(alpha_3+
M_PI/3.);
148 *ev0=trA-2.*sq_p*cos(alpha_3-
M_PI/3.);
179 void vectorInKernel2(
const double A00,
const double A10,
const double A01,
const double A11,
180 double* V0,
double*V1)
182 register double absA00=fabs(A00);
183 register double absA10=fabs(A10);
184 register double absA01=fabs(A01);
185 register double absA11=fabs(A11);
186 register double m=absA11>absA10 ? absA11 : absA10;
187 if (absA00>m || absA01>m) {
220 const double A01,
const double A11,
const double A21,
221 const double A02,
const double A12,
const double A22,
222 double* V0,
double* V1,
double* V2)
225 register const double I00=1./A00;
226 register const double IA10=I00*A10;
227 register const double IA20=I00*A20;
229 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
230 *V0=-(A10*TEMP0+A20*TEMP1);
254 double* ev0,
double* ev1,
255 double* V00,
double* V10,
double* V01,
double* V11,
260 const register double absev0=fabs(*ev0);
261 const register double absev1=fabs(*ev1);
262 register double max_ev=absev0>absev1 ? absev0 : absev1;
263 if (fabs((*ev0)-(*ev1))<tol*max_ev) {
270 const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
281 }
else if (TEMP0>0.) {
312 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
317 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
323 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
327 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
363 const double A11,
const double A12,
const double A22,
364 double* ev0,
double* ev1,
double* ev2,
365 double* V00,
double* V10,
double* V20,
366 double* V01,
double* V11,
double* V21,
367 double* V02,
double* V12,
double* V22,
370 register const double absA01=fabs(A01);
371 register const double absA02=fabs(A02);
372 register const double m=absA01>absA02 ? absA01 : absA02;
374 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
377 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
391 }
else if (A00>TEMP_EV1) {
420 const double absev0=fabs(*ev0);
421 const double absev1=fabs(*ev1);
422 const double absev2=fabs(*ev2);
423 double max_ev=absev0>absev1 ? absev0 : absev1;
424 max_ev=max_ev>absev2 ? max_ev : absev2;
425 const double d_01=fabs((*ev0)-(*ev1));
426 const double d_12=fabs((*ev1)-(*ev2));
427 const double max_d=d_01>d_12 ? d_01 : d_12;
428 if (max_d<=tol*max_ev) {
439 const double S00=A00-(*ev0);
440 const double absS00=fabs(S00);
442 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
443 }
else if (absA02<m) {
444 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
446 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
449 const double T00=A00-(*ev2);
450 const double absT00=fabs(T00);
452 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
453 }
else if (absA02<m) {
454 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
456 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
458 const double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
463 *V01=(*V10)*(*V22)-(*V12)*(*V20);
464 *V11=(*V20)*(*V02)-(*V00)*(*V22);
465 *V21=(*V00)*(*V12)-(*V02)*(*V10);
476 if (transpose == 0) {
477 for (
int i=0; i<SL; i++) {
478 for (
int j=0; j<SR; j++) {
480 for (
int l=0; l<SM; l++) {
481 sum += A[i+SL*l] * B[l+SM*j];
487 else if (transpose == 1) {
488 for (
int i=0; i<SL; i++) {
489 for (
int j=0; j<SR; j++) {
491 for (
int l=0; l<SM; l++) {
492 sum += A[i*SM+l] * B[l+SM*j];
498 else if (transpose == 2) {
499 for (
int i=0; i<SL; i++) {
500 for (
int j=0; j<SR; j++) {
502 for (
int l=0; l<SM; l++) {
503 sum += A[i+SL*l] * B[l*SR+j];
511 template <
typename UnaryFunction>
515 UnaryFunction operation)
517 for (
int i = 0; i < size; ++i) {
518 argRes[i] = operation(arg1[i]);
523 template <
typename BinaryFunction>
528 BinaryFunction operation)
530 for (
int i = 0; i < size; ++i) {
531 argRes[i] = operation(arg1[i], arg2[i]);
536 template <
typename BinaryFunction>
541 BinaryFunction operation)
543 for (
int i = 0; i < size; ++i) {
544 argRes[i] = operation(arg1, arg2[i]);
549 template <
typename BinaryFunction>
554 BinaryFunction operation)
556 for (
int i = 0; i < size; ++i) {
557 argRes[i] = operation(arg1[i], arg2);
Definition: AbstractContinuousDomain.cpp:24
void tensor_unary_operation(const int size, const double *arg1, double *argRes, UnaryFunction operation)
Definition: LocalOps.h:512
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 vectorInKernel2(const double A00, const double A10, const double A01, const double A11, double *V0, double *V1)
returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension i...
Definition: LocalOps.h:179
double makeNaN()
returns a NaN.
Definition: LocalOps.h:59
#define M_PI
Definition: LocalOps.h:22
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
bool nancheck(double d)
acts as a wrapper to isnan.
Definition: LocalOps.h:41
void scale(dim_t N, double *x, double a)
x = a*x
Definition: PasoUtil.h:93
void matrix_matrix_product(const int SL, const int SM, const int SR, const double *A, const double *B, double *C, int transpose)
Definition: LocalOps.h:474
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
void eigenvalues1(const double A00, double *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: LocalOps.h:78
void normalizeVector3(double *V0, double *V1, double *V2)
nomalizes a 3-d vector such that length is one and first non-zero component is positive.
Definition: LocalOps.h:308
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 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
void vectorInKernel3__nonZeroA00(const double A00, const double A10, const double A20, const double A01, const double A11, const double A21, const double A02, const double A12, const double A22, double *V0, double *V1, double *V2)
returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] assuming that ...
Definition: LocalOps.h:219
void tensor_binary_operation(const int size, const double *arg1, const double *arg2, double *argRes, BinaryFunction operation)
Definition: LocalOps.h:524