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 
28 #ifdef ESYS_USE_BOOST_ACOS
29 #include <boost/math/complex/acos.hpp> // std::acos for complex on OSX (elcapitan) is wrong
30 #endif
31 
32 #ifndef M_PI
33 # define M_PI 3.14159265358979323846 /* pi */
34 #endif
35 
36 
45 #include "ES_optype.h"
46 
47 namespace escript {
48 
49 
50 
51 bool always_real(escript::ES_optype operation);
52 
53 
58 struct FMax
59 {
61  {
62  return std::max(x,y);
63  }
67 };
68 
73 struct FMin
74 {
76  {
77  return std::min(x,y);
78  }
82 };
83 
88 template<typename T>
89 struct AbsMax
90 {
91  inline DataTypes::real_t operator()(T x, T y) const
92  {
93  return std::max(std::abs(x),std::abs(y));
94  }
98 };
99 
100 
101 inline
104 {
105  if (x == 0) {
106  return 0;
107  } else {
108  return x/fabs(x);
109  }
110 }
111 
116 inline
118 {
119  // Q: so why not just test d!=d?
120  // A: Coz it doesn't always work [I've checked].
121  // One theory is that the optimizer skips the test.
122  return std::isnan(d); // isNan should be a function in C++ land
123 }
124 
125 // I'm not sure if there is agreement about what a complex nan would be
126 // so, in this case I've just checked both components
127 inline
129 {
130  // Q: so why not just test d!=d?
131  // A: Coz it doesn't always work [I've checked].
132  // One theory is that the optimizer skips the test.
133  return std::isnan( real(d) ) || std::isnan( imag(d) ); // isNan should be a function in C++ land
134 }
135 
136 
141 inline
143 {
144 #ifdef nan
145  return nan("");
146 #else
147  return sqrt(-1.);
148 #endif
149 
150 }
151 
152 
160 inline
162  *ev0=A00;
163 }
164 
165 inline
167  *ev0=A00;
168 }
169 
170 
181 template <class T>
182 inline
183 void eigenvalues2(const T A00,const T A01,const T A11,
184  T* ev0, T* ev1) {
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);
189  *ev0=trA-s;
190  *ev1=trA+s;
191 }
206 inline
208  const DataTypes::real_t A11, const DataTypes::real_t A12,
209  const DataTypes::real_t A22,
211 
212  const DataTypes::real_t trA=(A00+A11+A22)/3.;
213  const DataTypes::real_t A_00=A00-trA;
214  const DataTypes::real_t A_11=A11-trA;
215  const DataTypes::real_t A_22=A22-trA;
216  const DataTypes::real_t A01_2=A01*A01;
217  const DataTypes::real_t A02_2=A02*A02;
218  const DataTypes::real_t A12_2=A12*A12;
219  const DataTypes::real_t p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
220  if (p<=0.) {
221  *ev2=trA;
222  *ev1=trA;
223  *ev0=trA;
224 
225  } else {
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);
227  const DataTypes::real_t sq_p=sqrt(p/3.);
228  DataTypes::real_t z=-q/(2*pow(sq_p,3));
229  if (z<-1.) {
230  z=-1.;
231  } else if (z>1.) {
232  z=1.;
233  }
234  const DataTypes::real_t alpha_3=acos(z)/3.;
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.);
238  }
239 }
249 inline
251 {
252  eigenvalues1(A00,ev0);
253  *V00=1.;
254  return;
255 }
267 inline
270 {
271  DataTypes::real_t absA00=fabs(A00);
272  DataTypes::real_t absA10=fabs(A10);
273  DataTypes::real_t absA01=fabs(A01);
274  DataTypes::real_t absA11=fabs(A11);
275  DataTypes::real_t m=absA11>absA10 ? absA11 : absA10;
276  if (absA00>m || absA01>m) {
277  *V0=-A01;
278  *V1=A00;
279  } else {
280  if (m<=0) {
281  *V0=1.;
282  *V1=0.;
283  } else {
284  *V0=A11;
285  *V1=-A10;
286  }
287  }
288 }
307 inline
309  const DataTypes::real_t A01,const DataTypes::real_t A11,const DataTypes::real_t A21,
310  const DataTypes::real_t A02,const DataTypes::real_t A12,const DataTypes::real_t A22,
312 {
313  DataTypes::real_t TEMP0,TEMP1;
314  const DataTypes::real_t I00=1./A00;
315  const DataTypes::real_t IA10=I00*A10;
316  const DataTypes::real_t IA20=I00*A20;
317  vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
318  A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
319  *V0=-(A10*TEMP0+A20*TEMP1);
320  *V1=A00*TEMP0;
321  *V2=A00*TEMP1;
322 }
323 
341 inline
345  const DataTypes::real_t tol)
346 {
347  DataTypes::real_t TEMP0,TEMP1;
348  eigenvalues2(A00,A01,A11,ev0,ev1);
349  const DataTypes::real_t absev0=fabs(*ev0);
350  const DataTypes::real_t absev1=fabs(*ev1);
351  DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
352  if (fabs((*ev0)-(*ev1))<tol*max_ev) {
353  *V00=1.;
354  *V10=0.;
355  *V01=0.;
356  *V11=1.;
357  } else {
358  vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
359  const DataTypes::real_t scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
360  if (TEMP0<0.) {
361  *V00=-TEMP0*scale;
362  *V10=-TEMP1*scale;
363  if (TEMP1<0.) {
364  *V01= *V10;
365  *V11=-(*V00);
366  } else {
367  *V01=-(*V10);
368  *V11= (*V00);
369  }
370  } else if (TEMP0>0.) {
371  *V00=TEMP0*scale;
372  *V10=TEMP1*scale;
373  if (TEMP1<0.) {
374  *V01=-(*V10);
375  *V11= (*V00);
376  } else {
377  *V01= (*V10);
378  *V11=-(*V00);
379  }
380  } else {
381  *V00=0.;
382  *V10=1;
383  *V11=0.;
384  *V01=1.;
385  }
386  }
387 }
396 inline
398 {
400  if (*V0>0) {
401  s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
402  *V0*=s;
403  *V1*=s;
404  *V2*=s;
405  } else if (*V0<0) {
406  s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
407  *V0*=s;
408  *V1*=s;
409  *V2*=s;
410  } else {
411  if (*V1>0) {
412  s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
413  *V1*=s;
414  *V2*=s;
415  } else if (*V1<0) {
416  s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
417  *V1*=s;
418  *V2*=s;
419  } else {
420  *V2=1.;
421  }
422  }
423 }
450 inline
452  const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22,
457  const DataTypes::real_t tol)
458 {
459  const DataTypes::real_t absA01=fabs(A01);
460  const DataTypes::real_t absA02=fabs(A02);
461  const DataTypes::real_t m=absA01>absA02 ? absA01 : absA02;
462  if (m<=0) {
463  DataTypes::real_t TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
464  eigenvalues_and_eigenvectors2(A11,A12,A22,
465  &TEMP_EV0,&TEMP_EV1,
466  &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
467  if (A00<=TEMP_EV0) {
468  *V00=1.;
469  *V10=0.;
470  *V20=0.;
471  *V01=0.;
472  *V11=TEMP_V00;
473  *V21=TEMP_V10;
474  *V02=0.;
475  *V12=TEMP_V01;
476  *V22=TEMP_V11;
477  *ev0=A00;
478  *ev1=TEMP_EV0;
479  *ev2=TEMP_EV1;
480  } else if (A00>TEMP_EV1) {
481  *V02=1.;
482  *V12=0.;
483  *V22=0.;
484  *V00=0.;
485  *V10=TEMP_V00;
486  *V20=TEMP_V10;
487  *V01=0.;
488  *V11=TEMP_V01;
489  *V21=TEMP_V11;
490  *ev0=TEMP_EV0;
491  *ev1=TEMP_EV1;
492  *ev2=A00;
493  } else {
494  *V01=1.;
495  *V11=0.;
496  *V21=0.;
497  *V00=0.;
498  *V10=TEMP_V00;
499  *V20=TEMP_V10;
500  *V02=0.;
501  *V12=TEMP_V01;
502  *V22=TEMP_V11;
503  *ev0=TEMP_EV0;
504  *ev1=A00;
505  *ev2=TEMP_EV1;
506  }
507  } else {
508  eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
509  const DataTypes::real_t absev0=fabs(*ev0);
510  const DataTypes::real_t absev1=fabs(*ev1);
511  const DataTypes::real_t absev2=fabs(*ev2);
512  DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
513  max_ev=max_ev>absev2 ? max_ev : absev2;
514  const DataTypes::real_t d_01=fabs((*ev0)-(*ev1));
515  const DataTypes::real_t d_12=fabs((*ev1)-(*ev2));
516  const DataTypes::real_t max_d=d_01>d_12 ? d_01 : d_12;
517  if (max_d<=tol*max_ev) {
518  *V00=1.;
519  *V10=0;
520  *V20=0;
521  *V01=0;
522  *V11=1.;
523  *V21=0;
524  *V02=0;
525  *V12=0;
526  *V22=1.;
527  } else {
528  const DataTypes::real_t S00=A00-(*ev0);
529  const DataTypes::real_t absS00=fabs(S00);
530  if (absS00>m) {
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);
534  } else {
535  vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
536  }
537  normalizeVector3(V00,V10,V20);;
538  const DataTypes::real_t T00=A00-(*ev2);
539  const DataTypes::real_t absT00=fabs(T00);
540  if (absT00>m) {
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);
544  } else {
545  vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
546  }
547  const DataTypes::real_t dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
548  *V02-=dot*(*V00);
549  *V12-=dot*(*V10);
550  *V22-=dot*(*V20);
551  normalizeVector3(V02,V12,V22);
552  *V01=(*V10)*(*V22)-(*V12)*(*V20);
553  *V11=(*V20)*(*V02)-(*V00)*(*V22);
554  *V21=(*V00)*(*V12)-(*V02)*(*V10);
555  normalizeVector3(V01,V11,V21);
556  }
557  }
558 }
559 
560 // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
561 // SM is the product of the last axis_offset entries in arg_0.getShape().
562 template <class LEFT, class RIGHT, class RES>
563 inline
564 void matrix_matrix_product(const int SL, const int SM, const int SR, const LEFT* A, const RIGHT* B, RES* C, int transpose)
565 {
566  if (transpose == 0) {
567  for (int i=0; i<SL; i++) {
568  for (int j=0; j<SR; j++) {
569  RES sum = 0.0;
570  for (int l=0; l<SM; l++) {
571  sum += A[i+SL*l] * B[l+SM*j];
572  }
573  C[i+SL*j] = sum;
574  }
575  }
576  }
577  else if (transpose == 1) {
578  for (int i=0; i<SL; i++) {
579  for (int j=0; j<SR; j++) {
580  RES sum = 0.0;
581  for (int l=0; l<SM; l++) {
582  sum += A[i*SM+l] * B[l+SM*j];
583  }
584  C[i+SL*j] = sum;
585  }
586  }
587  }
588  else if (transpose == 2) {
589  for (int i=0; i<SL; i++) {
590  for (int j=0; j<SR; j++) {
591  RES sum = 0.0;
592  for (int l=0; l<SM; l++) {
593  sum += A[i+SL*l] * B[l*SR+j];
594  }
595  C[i+SL*j] = sum;
596  }
597  }
598  }
599 }
600 
601 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
602 #else
603 
604 inline
606 {
607  return ::erf(x);
608 }
609 
610 inline
612 {
613  return makeNaN();
614 }
615 
616 #endif
617 
619 {
620  return escript::fsign(x);
621 }
622 
624 {
625  return makeNaN();
626 }
627 
628 inline
630 {
631  return acos(x);
632 }
633 
634 inline
636 {
637 #ifdef ESYS_USE_BOOST_ACOS
638  return boost::math::acos(x);
639 #else
640  return acos(x);
641 #endif
642 }
643 
644 
646 {
647  return abs(c);
648 }
649 
650 
651 
652 inline DataTypes::real_t calc_gtzero(const DataTypes::real_t& x) {return x>0;}
654 
655 
656 inline DataTypes::real_t calc_gezero(const DataTypes::real_t& x) {return x>=0;}
658 
659 
660 inline DataTypes::real_t calc_ltzero(const DataTypes::real_t& x) {return x<0;}
662 
663 inline DataTypes::real_t calc_lezero(const DataTypes::real_t& x) {return x<=0;}
665 
666 template <typename IN>
668 {
669  return fabs(i);
670 }
671 
672 template <>
674 {
675  return abs(i);
676 }
677 
678 
679 
680 
681 // deals with unary operations which return real, regardless of
682 // their input type
683 template <class IN>
684 inline void tensor_unary_array_operation_real(const size_t size,
685  const IN *arg1,
686  DataTypes::real_t * argRes,
687  escript::ES_optype operation,
688  DataTypes::real_t tol=0)
689 {
690  switch (operation)
691  {
692  case REAL:
693  for (int i = 0; i < size; ++i) {
694  argRes[i] = std::real(arg1[i]);
695  }
696  break;
697  case IMAG:
698  for (int i = 0; i < size; ++i) {
699  argRes[i] = std::imag(arg1[i]);
700  }
701  break;
702  case EZ:
703  for (size_t i = 0; i < size; ++i) {
704  argRes[i] = (fabs(arg1[i])<=tol);
705  }
706  break;
707  case NEZ:
708  for (size_t i = 0; i < size; ++i) {
709  argRes[i] = (fabs(arg1[i])>tol);
710  }
711  break;
712  case ABS:
713  for (size_t i = 0; i < size; ++i) {
714  argRes[i] = abs_f(arg1[i]);
715  }
716  break;
717  case PHS:
718  for (size_t i = 0; i < size; ++i) {
719  argRes[i] = std::arg(arg1[i]);
720  }
721  break;
722  default:
723  std::ostringstream oss;
724  oss << "Unsupported unary operation=";
725  oss << opToString(operation);
726  oss << '/';
727  oss << operation;
728  oss << " (Was expecting an operation with real results)";
729  throw DataException(oss.str());
730  }
731 }
732 
733 
734 
735 template <typename OUT, typename IN>
736 inline OUT conjugate(const IN i)
737 {
738  return conj(i);
739 }
740 
741 // This should never actually be called
742 template <>
744 {
745  return r;
746 }
747 
748 inline void tensor_unary_promote(const size_t size,
749  const DataTypes::real_t *arg1,
750  DataTypes::cplx_t * argRes)
751 {
752  for (size_t i = 0; i < size; ++i) {
753  argRes[i] = arg1[i];
754  }
755 }
756 
757 // No openmp because it's called by Lazy
758 // In most cases, IN and OUT will be the same
759 // but not ruling out putting Re() and Im()
760 // through this
761 template <class IN, typename OUT>
762 inline void tensor_unary_array_operation(const size_t size,
763  const IN *arg1,
764  OUT * argRes,
765  escript::ES_optype operation,
766  DataTypes::real_t tol=0)
767 {
768  switch (operation)
769  {
770  case NEG:
771  for (size_t i = 0; i < size; ++i) {
772  argRes[i] = -arg1[i];
773  }
774  break;
775  case SIN:
776  for (size_t i = 0; i < size; ++i) {
777  argRes[i] = sin(arg1[i]);
778  }
779  break;
780  case COS:
781  for (size_t i = 0; i < size; ++i) {
782  argRes[i] = cos(arg1[i]);
783  }
784  break;
785  case TAN:
786  for (size_t i = 0; i < size; ++i) {
787  argRes[i] = tan(arg1[i]);
788  }
789  break;
790  case ASIN:
791  for (size_t i = 0; i < size; ++i) {
792  argRes[i] = asin(arg1[i]);
793  }
794  break;
795  case ACOS:
796  for (size_t i = 0; i < size; ++i) {
797  argRes[i]=calc_acos(arg1[i]);
798  }
799  break;
800  case ATAN:
801  for (size_t i = 0; i < size; ++i) {
802  argRes[i] = atan(arg1[i]);
803  }
804  break;
805  case ABS:
806  for (size_t i = 0; i < size; ++i) {
807  argRes[i] = std::abs(arg1[i]);
808  }
809  break;
810  case SINH:
811  for (size_t i = 0; i < size; ++i) {
812  argRes[i] = sinh(arg1[i]);
813  }
814  break;
815  case COSH:
816  for (size_t i = 0; i < size; ++i) {
817  argRes[i] = cosh(arg1[i]);
818  }
819  break;
820  case TANH:
821  for (size_t i = 0; i < size; ++i) {
822  argRes[i] = tanh(arg1[i]);
823  }
824  break;
825  case ERF:
826  for (size_t i = 0; i < size; ++i) {
827  argRes[i] = calc_erf(arg1[i]);
828  }
829  break;
830  case ASINH:
831  for (size_t i = 0; i < size; ++i) {
832  argRes[i] = asinh(arg1[i]);
833  }
834  break;
835  case ACOSH:
836  for (size_t i = 0; i < size; ++i) {
837  argRes[i] = acosh(arg1[i]);
838  }
839  break;
840  case ATANH:
841  for (size_t i = 0; i < size; ++i) {
842  argRes[i] = atanh(arg1[i]);
843  }
844  break;
845  case LOG10:
846  for (size_t i = 0; i < size; ++i) {
847  argRes[i] = log10(arg1[i]);
848  }
849  break;
850  case LOG:
851  for (size_t i = 0; i < size; ++i) {
852  argRes[i] = log(arg1[i]);
853  }
854  break;
855  case SIGN:
856  for (size_t i = 0; i < size; ++i) {
857  argRes[i] = calc_sign(arg1[i]);
858  }
859  break;
860  case EXP:
861  for (size_t i = 0; i < size; ++i) {
862  argRes[i] = exp(arg1[i]);
863  }
864  break;
865  case SQRT:
866  for (size_t i = 0; i < size; ++i) {
867  argRes[i] = sqrt(arg1[i]);
868  }
869  break;
870  case GZ:
871  for (size_t i = 0; i < size; ++i) {
872  argRes[i] = calc_gtzero(arg1[i]);
873  }
874  break;
875  case GEZ:
876  for (size_t i = 0; i < size; ++i) {
877  argRes[i] = calc_gezero(arg1[i]);
878  }
879  break;
880  case LZ:
881  for (size_t i = 0; i < size; ++i) {
882  argRes[i] = calc_ltzero(arg1[i]);
883  }
884  break;
885  case LEZ:
886  for (size_t i = 0; i < size; ++i) {
887  argRes[i] = calc_lezero(arg1[i]);
888  }
889  break;
890  case CONJ:
891  for (size_t i = 0; i < size; ++i) {
892  argRes[i] = conjugate<OUT,IN>(arg1[i]);
893  }
894  break;
895  case RECIP:
896  for (size_t i = 0; i < size; ++i) {
897  argRes[i] = 1.0/arg1[i];
898  }
899  break;
900  case EZ:
901  for (size_t i = 0; i < size; ++i) {
902  argRes[i] = fabs(arg1[i])<=tol;
903  }
904  break;
905  case NEZ:
906  for (size_t i = 0; i < size; ++i) {
907  argRes[i] = fabs(arg1[i])>tol;
908  }
909  break;
910 
911  default:
912  std::ostringstream oss;
913  oss << "Unsupported unary operation=";
914  oss << opToString(operation);
915  oss << '/';
916  oss << operation;
917  throw DataException(oss.str());
918  }
919  return;
920 }
921 
922 bool supports_cplx(escript::ES_optype operation);
923 
924 
925 } // end of namespace
926 
927 #endif // __ESCRIPT_LOCALOPS_H__
928 
escript::tensor_unary_array_operation_real
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:684
escript::calc_acos
DataTypes::real_t calc_acos(DataTypes::real_t x)
Definition: ArrayOps.h:629
escript::FMin
Return the minimum value of the two given values.
Definition: ArrayOps.h:73
escript::AbsMax::operator()
DataTypes::real_t operator()(T x, T y) const
Definition: ArrayOps.h:91
escript::SIGN
@ SIGN
Definition: ES_optype.h:79
escript::DataTypes::real_t
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:79
escript::conjugate
OUT conjugate(const IN i)
Definition: ArrayOps.h:736
escript::CONJ
@ CONJ
Definition: ES_optype.h:105
escript::LZ
@ LZ
Definition: ES_optype.h:87
escript::COSH
@ COSH
Definition: ES_optype.h:71
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:161
escript::FMax
Return the maximum value of the two given values.
Definition: ArrayOps.h:58
escript::nancheck
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition: ArrayOps.h:117
escript::calc_gezero
DataTypes::real_t calc_gezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:656
escript::SQRT
@ SQRT
Definition: ES_optype.h:84
escript::IMAG
@ IMAG
Definition: ES_optype.h:104
escript::calc_erf
DataTypes::real_t calc_erf(DataTypes::real_t x)
Definition: ArrayOps.h:605
escript::EXP
@ EXP
Definition: ES_optype.h:83
escript::AbsMax::first_argument_type
T first_argument_type
Definition: ArrayOps.h:95
escript::calc_ltzero
DataTypes::real_t calc_ltzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:660
escript::SINH
@ SINH
Definition: ES_optype.h:70
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:76
escript::LOG
@ LOG
Definition: ES_optype.h:78
escript::abs_f
DataTypes::real_t abs_f(IN i)
Definition: ArrayOps.h:667
escript::fsign
DataTypes::real_t fsign(DataTypes::real_t x)
Definition: ArrayOps.h:103
escript::calc_lezero
DataTypes::real_t calc_lezero(const DataTypes::real_t &x)
Definition: ArrayOps.h:663
escript::ACOS
@ ACOS
Definition: ES_optype.h:68
escript::fabs
escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
Definition: ArrayOps.h:645
escript::AbsMax
Return the absolute maximum value of the two given values.
Definition: ArrayOps.h:89
escript::FMin::first_argument_type
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:79
escript::FMax::second_argument_type
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:65
escript::TANH
@ TANH
Definition: ES_optype.h:72
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:451
escript::PHS
@ PHS
Definition: ES_optype.h:110
escript::ERF
@ ERF
Definition: ES_optype.h:73
escript::NEG
@ NEG
Definition: ES_optype.h:81
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:183
escript::GZ
@ GZ
Definition: ES_optype.h:86
ArrayOps.h
escript::DataException
Definition: DataException.h:41
escript::FMax::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:66
escript::ABS
@ ABS
Definition: ES_optype.h:80
escript::AbsMax::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:97
escript::FMax::first_argument_type
DataTypes::real_t first_argument_type
Definition: ArrayOps.h:64
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:564
escript::REAL
@ REAL
Definition: ES_optype.h:103
escript::ACOSH
@ ACOSH
Definition: ES_optype.h:75
escript::GEZ
@ GEZ
Definition: ES_optype.h:88
escript::LEZ
@ LEZ
Definition: ES_optype.h:89
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:308
escript::FMin::second_argument_type
DataTypes::real_t second_argument_type
Definition: ArrayOps.h:80
escript::AbsMax::second_argument_type
T second_argument_type
Definition: ArrayOps.h:96
escript::FMin::operator()
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:75
escript::RECIP
@ RECIP
Definition: ES_optype.h:85
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:397
escript::ES_optype
ES_optype
Definition: ES_optype.h:41
paso::util::scale
void scale(dim_t N, double *x, double a)
x = a*x
Definition: PasoUtil.h:122
M_PI
#define M_PI
Definition: ArrayOps.h:33
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:250
escript
Definition: AbstractContinuousDomain.cpp:23
escript::LOG10
@ LOG10
Definition: ES_optype.h:77
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:342
escript::supports_cplx
bool supports_cplx(escript::ES_optype operation)
Definition: ArrayOps.cpp:26
escript::SIN
@ SIN
Definition: ES_optype.h:64
escript::EZ
@ EZ
Definition: ES_optype.h:91
ES_optype.h
escript::makeNaN
DataTypes::real_t makeNaN()
returns a NaN.
Definition: ArrayOps.h:142
escript::FMax::operator()
DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
Definition: ArrayOps.h:60
escript::TAN
@ TAN
Definition: ES_optype.h:66
escript::DataTypes::cplx_t
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:82
escript::calc_gtzero
DataTypes::real_t calc_gtzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:652
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:268
escript::tensor_unary_array_operation
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:762
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:207
escript::ASIN
@ ASIN
Definition: ES_optype.h:67
escript::COS
@ COS
Definition: ES_optype.h:65
escript::ASINH
@ ASINH
Definition: ES_optype.h:74
DataException.h
escript::FMin::result_type
DataTypes::real_t result_type
Definition: ArrayOps.h:81
escript::tensor_unary_promote
void tensor_unary_promote(const size_t size, const DataTypes::real_t *arg1, DataTypes::cplx_t *argRes)
Definition: ArrayOps.h:748
escript::ATAN
@ ATAN
Definition: ES_optype.h:69
escript::NEZ
@ NEZ
Definition: ES_optype.h:90
escript::calc_sign
DataTypes::real_t calc_sign(DataTypes::real_t x)
Definition: ArrayOps.h:618