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 
50 namespace escript {
51 
52 
53 
54 bool always_real(escript::ES_optype operation);
55 
56 
61 struct FMax
62 {
64  {
65  return std::max(x,y);
66  }
70 };
71 
76 struct FMin
77 {
79  {
80  return std::min(x,y);
81  }
85 };
86 
91 template<typename T>
92 struct 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 
104 inline
107 {
108  if (x == 0) {
109  return 0;
110  } else {
111  return x/fabs(x);
112  }
113 }
114 
119 inline
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
130 inline
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 
144 inline
146 {
147 #ifdef nan
148  return nan("");
149 #else
150  return sqrt(-1.);
151 #endif
152 
153 }
154 
155 
163 inline
165  *ev0=A00;
166 }
167 
168 inline
170  *ev0=A00;
171 }
172 
173 
184 template <class T>
185 inline
186 void 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 }
209 inline
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 }
252 inline
254 {
255  eigenvalues1(A00,ev0);
256  *V00=1.;
257  return;
258 }
270 inline
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 }
310 inline
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 
344 inline
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 }
399 inline
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 }
453 inline
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;
467  eigenvalues_and_eigenvectors2(A11,A12,A22,
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().
565 template <class LEFT, class RIGHT, class RES>
566 inline
567 void 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 
604 inline
606 {
607  return ::erf(x);
608 }
609 
610 inline
612 {
613  return makeNaN();
614 }
615 
617 {
618  return escript::fsign(x);
619 }
620 
622 {
623  return makeNaN();
624 }
625 
626 inline
628 {
629  return acos(x);
630 }
631 
632 inline
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 
650 inline DataTypes::real_t calc_gtzero(const DataTypes::real_t& x) {return x>0;}
652 
653 
654 inline DataTypes::real_t calc_gezero(const DataTypes::real_t& x) {return x>=0;}
656 
657 
658 inline DataTypes::real_t calc_ltzero(const DataTypes::real_t& x) {return x<0;}
660 
661 inline DataTypes::real_t calc_lezero(const DataTypes::real_t& x) {return x<=0;}
663 
664 template <typename T>
666 {
667  return fabs(i);
668 }
669 
670 template <>
672 {
673  return abs(i);
674 }
675 
676 
677 
678 
679 // deals with unary operations which return real, regardless of
680 // their input type
681 template <class T>
682 inline 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 
733 template <typename T1, typename T2>
734 inline T1 conjugate(const T2 i)
735 {
736  return conj(i);
737 }
738 
739 // This should never actually be called
740 template <>
742 {
743  return r;
744 }
745 
746 inline 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
759 template <class T1, typename T2>
760 inline 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 
920 bool supports_cplx(escript::ES_optype operation);
921 
922 
923 } // end of namespace
924 
925 #endif // __ESCRIPT_LOCALOPS_H__
926 
escript::calc_acos
DataTypes::real_t calc_acos(DataTypes::real_t x)
Definition: ArrayOps.h:627
escript::FMin
Return the minimum value of the two given values.
Definition: ArrayOps.h:77
escript::AbsMax::operator()
DataTypes::real_t operator()(T x, T y) const
Definition: ArrayOps.h:94
escript::SIGN
@ SIGN
Definition: ES_optype.h:53
escript::DataTypes::real_t
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:52
escript::CONJ
@ CONJ
Definition: ES_optype.h:79
escript::LZ
@ LZ
Definition: ES_optype.h:61
escript::COSH
@ COSH
Definition: ES_optype.h:45
escript::eigenvalues1
void eigenvalues1(const DataTypes::real_t A00, DataTypes::real_t *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: ArrayOps.h:164
escript::FMax
Return the maximum value of the two given values.
Definition: ArrayOps.h:62
escript::nancheck
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition: ArrayOps.h:120
escript::calc_gezero
DataTypes::real_t calc_gezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:654
escript::SQRT
@ SQRT
Definition: ES_optype.h:58
escript::IMAG
@ IMAG
Definition: ES_optype.h:78
escript::calc_erf
DataTypes::real_t calc_erf(DataTypes::real_t x)
Definition: ArrayOps.h:605
escript::EXP
@ EXP
Definition: ES_optype.h:57
escript::AbsMax::first_argument_type
T first_argument_type
Definition: ArrayOps.h:98
escript::calc_ltzero
DataTypes::real_t calc_ltzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:658
escript::SINH
@ SINH
Definition: ES_optype.h:44
escript::transpose
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
escript::opToString
const std::string & opToString(ES_optype op)
Definition: ES_optype.cpp:89
escript::ATANH
@ ATANH
Definition: ES_optype.h:50
escript::LOG
@ LOG
Definition: ES_optype.h:52
escript::fsign
DataTypes::real_t fsign(DataTypes::real_t x)
Definition: ArrayOps.h:106
escript::calc_lezero
DataTypes::real_t calc_lezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:661
escript::ACOS
@ ACOS
Definition: ES_optype.h:42
escript::fabs
escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
Definition: ArrayOps.h:643
escript::tensor_unary_array_operation
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
escript::AbsMax
Return the absolute maximum value of the two given values.
Definition: ArrayOps.h:93
escript::FMin::first_argument_type
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:82
escript::FMax::second_argument_type
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:68
escript::TANH
@ TANH
Definition: ES_optype.h:46
escript::eigenvalues_and_eigenvectors3
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
escript::PHS
@ PHS
Definition: ES_optype.h:84
escript::ERF
@ ERF
Definition: ES_optype.h:47
escript::NEG
@ NEG
Definition: ES_optype.h:55
escript::eigenvalues2
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
escript::GZ
@ GZ
Definition: ES_optype.h:60
ArrayOps.h
escript::DataException
Definition: DataException.h:28
escript::FMax::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:69
escript::ABS
@ ABS
Definition: ES_optype.h:54
escript::AbsMax::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:100
escript::FMax::first_argument_type
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:67
escript::matrix_matrix_product
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::REAL
@ REAL
Definition: ES_optype.h:77
escript::ACOSH
@ ACOSH
Definition: ES_optype.h:49
escript::abs_f
DataTypes::real_t abs_f(T i)
Definition: ArrayOps.h:665
escript::GEZ
@ GEZ
Definition: ES_optype.h:62
escript::LEZ
@ LEZ
Definition: ES_optype.h:63
escript::vectorInKernel3__nonZeroA00
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
escript::FMin::second_argument_type
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:83
escript::AbsMax::second_argument_type
T second_argument_type
Definition: ArrayOps.h:99
escript::FMin::operator()
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:78
escript::RECIP
@ RECIP
Definition: ES_optype.h:59
escript::normalizeVector3
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
escript::ES_optype
ES_optype
Definition: ES_optype.h:30
escript::conjugate
T1 conjugate(const T2 i)
Definition: ArrayOps.h:734
paso::util::scale
void scale(dim_t N, double *x, double a)
x = a*x
Definition: PasoUtil.h:94
M_PI
#define M_PI
Definition: ArrayOps.h:36
escript::eigenvalues_and_eigenvectors1
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
escript
Definition: AbstractContinuousDomain.cpp:23
escript::LOG10
@ LOG10
Definition: ES_optype.h:51
DataTypes.h
escript::always_real
bool always_real(escript::ES_optype operation)
Definition: ArrayOps.cpp:66
escript::eigenvalues_and_eigenvectors2
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
escript::supports_cplx
bool supports_cplx(escript::ES_optype operation)
Definition: ArrayOps.cpp:26
escript::SIN
@ SIN
Definition: ES_optype.h:38
escript::EZ
@ EZ
Definition: ES_optype.h:65
ES_optype.h
escript::makeNaN
DataTypes::real_t makeNaN()
returns a NaN.
Definition: ArrayOps.h:145
escript::FMax::operator()
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:63
escript::TAN
@ TAN
Definition: ES_optype.h:40
escript::DataTypes::cplx_t
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:55
escript::tensor_unary_array_operation_real
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
escript::calc_gtzero
DataTypes::real_t calc_gtzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:650
escript::vectorInKernel2
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
escript::eigenvalues3
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
escript::ASIN
@ ASIN
Definition: ES_optype.h:41
escript::COS
@ COS
Definition: ES_optype.h:39
escript::ASINH
@ ASINH
Definition: ES_optype.h:48
DataException.h
escript::FMin::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:84
escript::tensor_unary_promote
void tensor_unary_promote(const size_t size, const DataTypes::real_t *arg1, DataTypes::cplx_t *argRes)
Definition: ArrayOps.h:746
escript::ATAN
@ ATAN
Definition: ES_optype.h:43
escript::NEZ
@ NEZ
Definition: ES_optype.h:64
escript::calc_sign
DataTypes::real_t calc_sign(DataTypes::real_t x)
Definition: ArrayOps.h:616