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