escript Revision_
ArrayOps.h
Go to the documentation of this file.
1
2/*****************************************************************************
3*
4* Copyright (c) 2003-2020 by The University of Queensland
5* http://www.uq.edu.au
6*
7* Primary Business: Queensland, Australia
8* Licensed under the Apache License, version 2.0
9* http://www.apache.org/licenses/LICENSE-2.0
10*
11* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12* Development 2012-2013 by School of Earth Sciences
13* Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14* Development from 2019 by School of Earth and Environmental Sciences
15**
16*****************************************************************************/
17
18
19#ifndef __ESCRIPT_LOCALOPS_H__
20#define __ESCRIPT_LOCALOPS_H__
21
22#include "DataTypes.h"
23#include "DataException.h"
24#include <iostream>
25#include <cmath>
26#include <complex>
27#ifdef _WIN32
28# include <algorithm>
29#endif // _WIN32
30
31#ifdef ESYS_USE_BOOST_ACOS
32#include <boost/math/complex/acos.hpp> // std::acos for complex on OSX (elcapitan) is wrong
33#endif
34
35#ifndef M_PI
36# define M_PI 3.14159265358979323846 /* pi */
37#endif
38
39
48#include "ES_optype.h"
49
50namespace escript {
51
52
53
54bool always_real(escript::ES_optype operation);
55
56
71
86
91template<typename T>
92struct AbsMax
93{
94 inline DataTypes::real_t operator()(T x, T y) const
95 {
96 return std::max(std::abs(x),std::abs(y));
97 }
101};
102
103
104inline
107{
108 if (x == 0) {
109 return 0;
110 } else {
111 return x/fabs(x);
112 }
113}
114
119inline
121{
122 // Q: so why not just test d!=d?
123 // A: Coz it doesn't always work [I've checked].
124 // One theory is that the optimizer skips the test.
125 return std::isnan(d); // isNan should be a function in C++ land
126}
127
128// I'm not sure if there is agreement about what a complex nan would be
129// so, in this case I've just checked both components
130inline
132{
133 // Q: so why not just test d!=d?
134 // A: Coz it doesn't always work [I've checked].
135 // One theory is that the optimizer skips the test.
136 return std::isnan( real(d) ) || std::isnan( imag(d) ); // isNan should be a function in C++ land
137}
138
139
144inline
146{
147#ifdef nan
148 return nan("");
149#else
150 return sqrt(-1.);
151#endif
152
153}
154
155
163inline
165 *ev0=A00;
166}
167
168inline
170 *ev0=A00;
171}
172
173
184template <class T>
185inline
186void eigenvalues2(const T A00,const T A01,const T A11,
187 T* ev0, T* ev1) {
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);
192 *ev0=trA-s;
193 *ev1=trA+s;
194}
209inline
211 const DataTypes::real_t A11, const DataTypes::real_t A12,
212 const DataTypes::real_t A22,
214
215 const DataTypes::real_t trA=(A00+A11+A22)/3.;
216 const DataTypes::real_t A_00=A00-trA;
217 const DataTypes::real_t A_11=A11-trA;
218 const DataTypes::real_t A_22=A22-trA;
219 const DataTypes::real_t A01_2=A01*A01;
220 const DataTypes::real_t A02_2=A02*A02;
221 const DataTypes::real_t A12_2=A12*A12;
222 const DataTypes::real_t p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
223 if (p<=0.) {
224 *ev2=trA;
225 *ev1=trA;
226 *ev0=trA;
227
228 } else {
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);
230 const DataTypes::real_t sq_p=sqrt(p/3.);
231 DataTypes::real_t z=-q/(2*pow(sq_p,3));
232 if (z<-1.) {
233 z=-1.;
234 } else if (z>1.) {
235 z=1.;
236 }
237 const DataTypes::real_t alpha_3=acos(z)/3.;
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.);
241 }
242}
252inline
254{
255 eigenvalues1(A00,ev0);
256 *V00=1.;
257 return;
258}
270inline
273{
274 DataTypes::real_t absA00=fabs(A00);
275 DataTypes::real_t absA10=fabs(A10);
276 DataTypes::real_t absA01=fabs(A01);
277 DataTypes::real_t absA11=fabs(A11);
278 DataTypes::real_t m=absA11>absA10 ? absA11 : absA10;
279 if (absA00>m || absA01>m) {
280 *V0=-A01;
281 *V1=A00;
282 } else {
283 if (m<=0) {
284 *V0=1.;
285 *V1=0.;
286 } else {
287 *V0=A11;
288 *V1=-A10;
289 }
290 }
291}
310inline
312 const DataTypes::real_t A01,const DataTypes::real_t A11,const DataTypes::real_t A21,
313 const DataTypes::real_t A02,const DataTypes::real_t A12,const DataTypes::real_t A22,
315{
316 DataTypes::real_t TEMP0,TEMP1;
317 const DataTypes::real_t I00=1./A00;
318 const DataTypes::real_t IA10=I00*A10;
319 const DataTypes::real_t IA20=I00*A20;
320 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
321 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
322 *V0=-(A10*TEMP0+A20*TEMP1);
323 *V1=A00*TEMP0;
324 *V2=A00*TEMP1;
325}
326
344inline
348 const DataTypes::real_t tol)
349{
350 DataTypes::real_t TEMP0,TEMP1;
351 eigenvalues2(A00,A01,A11,ev0,ev1);
352 const DataTypes::real_t absev0=fabs(*ev0);
353 const DataTypes::real_t absev1=fabs(*ev1);
354 DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
355 if (fabs((*ev0)-(*ev1))<tol*max_ev) {
356 *V00=1.;
357 *V10=0.;
358 *V01=0.;
359 *V11=1.;
360 } else {
361 vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
362 const DataTypes::real_t scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
363 if (TEMP0<0.) {
364 *V00=-TEMP0*scale;
365 *V10=-TEMP1*scale;
366 if (TEMP1<0.) {
367 *V01= *V10;
368 *V11=-(*V00);
369 } else {
370 *V01=-(*V10);
371 *V11= (*V00);
372 }
373 } else if (TEMP0>0.) {
374 *V00=TEMP0*scale;
375 *V10=TEMP1*scale;
376 if (TEMP1<0.) {
377 *V01=-(*V10);
378 *V11= (*V00);
379 } else {
380 *V01= (*V10);
381 *V11=-(*V00);
382 }
383 } else {
384 *V00=0.;
385 *V10=1;
386 *V11=0.;
387 *V01=1.;
388 }
389 }
390}
399inline
401{
403 if (*V0>0) {
404 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
405 *V0*=s;
406 *V1*=s;
407 *V2*=s;
408 } else if (*V0<0) {
409 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
410 *V0*=s;
411 *V1*=s;
412 *V2*=s;
413 } else {
414 if (*V1>0) {
415 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
416 *V1*=s;
417 *V2*=s;
418 } else if (*V1<0) {
419 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
420 *V1*=s;
421 *V2*=s;
422 } else {
423 *V2=1.;
424 }
425 }
426}
453inline
455 const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22,
460 const DataTypes::real_t tol)
461{
462 const DataTypes::real_t absA01=fabs(A01);
463 const DataTypes::real_t absA02=fabs(A02);
464 const DataTypes::real_t m=absA01>absA02 ? absA01 : absA02;
465 if (m<=0) {
466 DataTypes::real_t TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
468 &TEMP_EV0,&TEMP_EV1,
469 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
470 if (A00<=TEMP_EV0) {
471 *V00=1.;
472 *V10=0.;
473 *V20=0.;
474 *V01=0.;
475 *V11=TEMP_V00;
476 *V21=TEMP_V10;
477 *V02=0.;
478 *V12=TEMP_V01;
479 *V22=TEMP_V11;
480 *ev0=A00;
481 *ev1=TEMP_EV0;
482 *ev2=TEMP_EV1;
483 } else if (A00>TEMP_EV1) {
484 *V02=1.;
485 *V12=0.;
486 *V22=0.;
487 *V00=0.;
488 *V10=TEMP_V00;
489 *V20=TEMP_V10;
490 *V01=0.;
491 *V11=TEMP_V01;
492 *V21=TEMP_V11;
493 *ev0=TEMP_EV0;
494 *ev1=TEMP_EV1;
495 *ev2=A00;
496 } else {
497 *V01=1.;
498 *V11=0.;
499 *V21=0.;
500 *V00=0.;
501 *V10=TEMP_V00;
502 *V20=TEMP_V10;
503 *V02=0.;
504 *V12=TEMP_V01;
505 *V22=TEMP_V11;
506 *ev0=TEMP_EV0;
507 *ev1=A00;
508 *ev2=TEMP_EV1;
509 }
510 } else {
511 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
512 const DataTypes::real_t absev0=fabs(*ev0);
513 const DataTypes::real_t absev1=fabs(*ev1);
514 const DataTypes::real_t absev2=fabs(*ev2);
515 DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
516 max_ev=max_ev>absev2 ? max_ev : absev2;
517 const DataTypes::real_t d_01=fabs((*ev0)-(*ev1));
518 const DataTypes::real_t d_12=fabs((*ev1)-(*ev2));
519 const DataTypes::real_t max_d=d_01>d_12 ? d_01 : d_12;
520 if (max_d<=tol*max_ev) {
521 *V00=1.;
522 *V10=0;
523 *V20=0;
524 *V01=0;
525 *V11=1.;
526 *V21=0;
527 *V02=0;
528 *V12=0;
529 *V22=1.;
530 } else {
531 const DataTypes::real_t S00=A00-(*ev0);
532 const DataTypes::real_t absS00=fabs(S00);
533 if (absS00>m) {
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);
537 } else {
538 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
539 }
540 normalizeVector3(V00,V10,V20);;
541 const DataTypes::real_t T00=A00-(*ev2);
542 const DataTypes::real_t absT00=fabs(T00);
543 if (absT00>m) {
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);
547 } else {
548 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
549 }
550 const DataTypes::real_t dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
551 *V02-=dot*(*V00);
552 *V12-=dot*(*V10);
553 *V22-=dot*(*V20);
554 normalizeVector3(V02,V12,V22);
555 *V01=(*V10)*(*V22)-(*V12)*(*V20);
556 *V11=(*V20)*(*V02)-(*V00)*(*V22);
557 *V21=(*V00)*(*V12)-(*V02)*(*V10);
558 normalizeVector3(V01,V11,V21);
559 }
560 }
561}
562
563// General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
564// SM is the product of the last axis_offset entries in arg_0.getShape().
565template <class LEFT, class RIGHT, class RES>
566inline
567void matrix_matrix_product(const int SL, const int SM, const int SR, const LEFT* A, const RIGHT* B, RES* C, int transpose)
568{
569 if (transpose == 0) {
570 for (int i=0; i<SL; i++) {
571 for (int j=0; j<SR; j++) {
572 RES sum = 0.0;
573 for (int l=0; l<SM; l++) {
574 sum += A[i+SL*l] * B[l+SM*j];
575 }
576 C[i+SL*j] = sum;
577 }
578 }
579 }
580 else if (transpose == 1) {
581 for (int i=0; i<SL; i++) {
582 for (int j=0; j<SR; j++) {
583 RES sum = 0.0;
584 for (int l=0; l<SM; l++) {
585 sum += A[i*SM+l] * B[l+SM*j];
586 }
587 C[i+SL*j] = sum;
588 }
589 }
590 }
591 else if (transpose == 2) {
592 for (int i=0; i<SL; i++) {
593 for (int j=0; j<SR; j++) {
594 RES sum = 0.0;
595 for (int l=0; l<SM; l++) {
596 sum += A[i+SL*l] * B[l*SR+j];
597 }
598 C[i+SL*j] = sum;
599 }
600 }
601 }
602}
603
604inline
606{
607 return ::erf(x);
608}
609
610inline
615
620
622{
623 return makeNaN();
624}
625
626inline
628{
629 return acos(x);
630}
631
632inline
634{
635#ifdef ESYS_USE_BOOST_ACOS
636 return boost::math::acos(x);
637#else
638 return acos(x);
639#endif
640}
641
642
644{
645 return abs(c);
646}
647
648
649
650inline DataTypes::real_t calc_gtzero(const DataTypes::real_t& x) {return x>0;}
652
653
654inline DataTypes::real_t calc_gezero(const DataTypes::real_t& x) {return x>=0;}
656
657
658inline DataTypes::real_t calc_ltzero(const DataTypes::real_t& x) {return x<0;}
660
661inline DataTypes::real_t calc_lezero(const DataTypes::real_t& x) {return x<=0;}
663
664template <typename T>
666{
667 return fabs(i);
668}
669
670template <>
672{
673 return abs(i);
674}
675
676
677
678
679// deals with unary operations which return real, regardless of
680// their input type
681template <class T>
682inline void tensor_unary_array_operation_real(const size_t size,
683 const T *arg1,
684 DataTypes::real_t * argRes,
685 escript::ES_optype operation,
686 DataTypes::real_t tol=0)
687{
688 switch (operation)
689 {
690 case REAL:
691 for (int i = 0; i < size; ++i) {
692 argRes[i] = std::real(arg1[i]);
693 }
694 break;
695 case IMAG:
696 for (int i = 0; i < size; ++i) {
697 argRes[i] = std::imag(arg1[i]);
698 }
699 break;
700 case EZ:
701 for (size_t i = 0; i < size; ++i) {
702 argRes[i] = (fabs(arg1[i])<=tol);
703 }
704 break;
705 case NEZ:
706 for (size_t i = 0; i < size; ++i) {
707 argRes[i] = (fabs(arg1[i])>tol);
708 }
709 break;
710 case ABS:
711 for (size_t i = 0; i < size; ++i) {
712 argRes[i] = abs_f(arg1[i]);
713 }
714 break;
715 case PHS:
716 for (size_t i = 0; i < size; ++i) {
717 argRes[i] = std::arg(arg1[i]);
718 }
719 break;
720 default:
721 std::ostringstream oss;
722 oss << "Unsupported unary operation=";
723 oss << opToString(operation);
724 oss << '/';
725 oss << operation;
726 oss << " (Was expecting an operation with real results)";
727 throw DataException(oss.str());
728 }
729}
730
731
732
733template <typename T1, typename T2>
734inline T1 conjugate(const T2 i)
735{
736 return conj(i);
737}
738
739// This should never actually be called
740template <>
742{
743 return r;
744}
745
746inline void tensor_unary_promote(const size_t size,
747 const DataTypes::real_t *arg1,
748 DataTypes::cplx_t * argRes)
749{
750 for (size_t i = 0; i < size; ++i) {
751 argRes[i] = arg1[i];
752 }
753}
754
755// No openmp because it's called by Lazy
756// In most cases, T1 and T2 will be the same
757// but not ruling out putting Re() and Im()
758// through this
759template <class T1, typename T2>
760inline void tensor_unary_array_operation(const size_t size,
761 const T1 *arg1,
762 T2 * argRes,
763 escript::ES_optype operation,
764 DataTypes::real_t tol=0)
765{
766 switch (operation)
767 {
768 case NEG:
769 for (size_t i = 0; i < size; ++i) {
770 argRes[i] = -arg1[i];
771 }
772 break;
773 case SIN:
774 for (size_t i = 0; i < size; ++i) {
775 argRes[i] = sin(arg1[i]);
776 }
777 break;
778 case COS:
779 for (size_t i = 0; i < size; ++i) {
780 argRes[i] = cos(arg1[i]);
781 }
782 break;
783 case TAN:
784 for (size_t i = 0; i < size; ++i) {
785 argRes[i] = tan(arg1[i]);
786 }
787 break;
788 case ASIN:
789 for (size_t i = 0; i < size; ++i) {
790 argRes[i] = asin(arg1[i]);
791 }
792 break;
793 case ACOS:
794 for (size_t i = 0; i < size; ++i) {
795 argRes[i]=calc_acos(arg1[i]);
796 }
797 break;
798 case ATAN:
799 for (size_t i = 0; i < size; ++i) {
800 argRes[i] = atan(arg1[i]);
801 }
802 break;
803 case ABS:
804 for (size_t i = 0; i < size; ++i) {
805 argRes[i] = std::abs(arg1[i]);
806 }
807 break;
808 case SINH:
809 for (size_t i = 0; i < size; ++i) {
810 argRes[i] = sinh(arg1[i]);
811 }
812 break;
813 case COSH:
814 for (size_t i = 0; i < size; ++i) {
815 argRes[i] = cosh(arg1[i]);
816 }
817 break;
818 case TANH:
819 for (size_t i = 0; i < size; ++i) {
820 argRes[i] = tanh(arg1[i]);
821 }
822 break;
823 case ERF:
824 for (size_t i = 0; i < size; ++i) {
825 argRes[i] = calc_erf(arg1[i]);
826 }
827 break;
828 case ASINH:
829 for (size_t i = 0; i < size; ++i) {
830 argRes[i] = asinh(arg1[i]);
831 }
832 break;
833 case ACOSH:
834 for (size_t i = 0; i < size; ++i) {
835 argRes[i] = acosh(arg1[i]);
836 }
837 break;
838 case ATANH:
839 for (size_t i = 0; i < size; ++i) {
840 argRes[i] = atanh(arg1[i]);
841 }
842 break;
843 case LOG10:
844 for (size_t i = 0; i < size; ++i) {
845 argRes[i] = log10(arg1[i]);
846 }
847 break;
848 case LOG:
849 for (size_t i = 0; i < size; ++i) {
850 argRes[i] = log(arg1[i]);
851 }
852 break;
853 case SIGN:
854 for (size_t i = 0; i < size; ++i) {
855 argRes[i] = calc_sign(arg1[i]);
856 }
857 break;
858 case EXP:
859 for (size_t i = 0; i < size; ++i) {
860 argRes[i] = exp(arg1[i]);
861 }
862 break;
863 case SQRT:
864 for (size_t i = 0; i < size; ++i) {
865 argRes[i] = sqrt(arg1[i]);
866 }
867 break;
868 case GZ:
869 for (size_t i = 0; i < size; ++i) {
870 argRes[i] = calc_gtzero(arg1[i]);
871 }
872 break;
873 case GEZ:
874 for (size_t i = 0; i < size; ++i) {
875 argRes[i] = calc_gezero(arg1[i]);
876 }
877 break;
878 case LZ:
879 for (size_t i = 0; i < size; ++i) {
880 argRes[i] = calc_ltzero(arg1[i]);
881 }
882 break;
883 case LEZ:
884 for (size_t i = 0; i < size; ++i) {
885 argRes[i] = calc_lezero(arg1[i]);
886 }
887 break;
888 case CONJ:
889 for (size_t i = 0; i < size; ++i) {
890 argRes[i] = conjugate<T2,T1>(arg1[i]);
891 }
892 break;
893 case RECIP:
894 for (size_t i = 0; i < size; ++i) {
895 argRes[i] = 1.0/arg1[i];
896 }
897 break;
898 case EZ:
899 for (size_t i = 0; i < size; ++i) {
900 argRes[i] = fabs(arg1[i])<=tol;
901 }
902 break;
903 case NEZ:
904 for (size_t i = 0; i < size; ++i) {
905 argRes[i] = fabs(arg1[i])>tol;
906 }
907 break;
908
909 default:
910 std::ostringstream oss;
911 oss << "Unsupported unary operation=";
912 oss << opToString(operation);
913 oss << '/';
914 oss << operation;
915 throw DataException(oss.str());
916 }
917 return;
918}
919
920bool supports_cplx(escript::ES_optype operation);
921
922
923} // end of namespace
924
925#endif // __ESCRIPT_LOCALOPS_H__
926
#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 result_type
Definition ArrayOps.h:100
DataTypes::real_t operator()(T x, T y) const
Definition ArrayOps.h:94
T first_argument_type
Definition ArrayOps.h:98
T second_argument_type
Definition ArrayOps.h:99
Return the maximum value of the two given values.
Definition ArrayOps.h:62
DataTypes::real_t first_argument_type
Definition ArrayOps.h:67
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
Return the minimum value of the two given values.
Definition ArrayOps.h:77
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition ArrayOps.h:78
DataTypes::real_t result_type
Definition ArrayOps.h:84
DataTypes::real_t first_argument_type
Definition ArrayOps.h:82
DataTypes::real_t second_argument_type
Definition ArrayOps.h:83