escript  Revision_
ArrayOps.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 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 
128 inline
130 {
131 #ifdef nan
132  return nan("");
133 #else
134  return sqrt(-1.);
135 #endif
136 
137 }
138 
139 
147 inline
149  *ev0=A00;
150 }
151 
152 inline
154  *ev0=A00;
155 }
156 
157 
168 template <class T>
169 inline
170 void eigenvalues2(const T A00,const T A01,const T A11,
171  T* ev0, T* ev1) {
172  const T trA=(A00+A11)/2.;
173  const T A_00=A00-trA;
174  const T A_11=A11-trA;
175  const T s=sqrt(A01*A01-A_00*A_11);
176  *ev0=trA-s;
177  *ev1=trA+s;
178 }
193 inline
195  const DataTypes::real_t A11, const DataTypes::real_t A12,
196  const DataTypes::real_t A22,
198 
199  const DataTypes::real_t trA=(A00+A11+A22)/3.;
200  const DataTypes::real_t A_00=A00-trA;
201  const DataTypes::real_t A_11=A11-trA;
202  const DataTypes::real_t A_22=A22-trA;
203  const DataTypes::real_t A01_2=A01*A01;
204  const DataTypes::real_t A02_2=A02*A02;
205  const DataTypes::real_t A12_2=A12*A12;
206  const DataTypes::real_t p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
207  if (p<=0.) {
208  *ev2=trA;
209  *ev1=trA;
210  *ev0=trA;
211 
212  } else {
213  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);
214  const DataTypes::real_t sq_p=sqrt(p/3.);
215  DataTypes::real_t z=-q/(2*pow(sq_p,3));
216  if (z<-1.) {
217  z=-1.;
218  } else if (z>1.) {
219  z=1.;
220  }
221  const DataTypes::real_t alpha_3=acos(z)/3.;
222  *ev2=trA+2.*sq_p*cos(alpha_3);
223  *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
224  *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
225  }
226 }
236 inline
238 {
239  eigenvalues1(A00,ev0);
240  *V00=1.;
241  return;
242 }
254 inline
257 {
258  DataTypes::real_t absA00=fabs(A00);
259  DataTypes::real_t absA10=fabs(A10);
260  DataTypes::real_t absA01=fabs(A01);
261  DataTypes::real_t absA11=fabs(A11);
262  DataTypes::real_t m=absA11>absA10 ? absA11 : absA10;
263  if (absA00>m || absA01>m) {
264  *V0=-A01;
265  *V1=A00;
266  } else {
267  if (m<=0) {
268  *V0=1.;
269  *V1=0.;
270  } else {
271  *V0=A11;
272  *V1=-A10;
273  }
274  }
275 }
294 inline
296  const DataTypes::real_t A01,const DataTypes::real_t A11,const DataTypes::real_t A21,
297  const DataTypes::real_t A02,const DataTypes::real_t A12,const DataTypes::real_t A22,
299 {
300  DataTypes::real_t TEMP0,TEMP1;
301  const DataTypes::real_t I00=1./A00;
302  const DataTypes::real_t IA10=I00*A10;
303  const DataTypes::real_t IA20=I00*A20;
304  vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
305  A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
306  *V0=-(A10*TEMP0+A20*TEMP1);
307  *V1=A00*TEMP0;
308  *V2=A00*TEMP1;
309 }
310 
328 inline
332  const DataTypes::real_t tol)
333 {
334  DataTypes::real_t TEMP0,TEMP1;
335  eigenvalues2(A00,A01,A11,ev0,ev1);
336  const DataTypes::real_t absev0=fabs(*ev0);
337  const DataTypes::real_t absev1=fabs(*ev1);
338  DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
339  if (fabs((*ev0)-(*ev1))<tol*max_ev) {
340  *V00=1.;
341  *V10=0.;
342  *V01=0.;
343  *V11=1.;
344  } else {
345  vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
346  const DataTypes::real_t scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
347  if (TEMP0<0.) {
348  *V00=-TEMP0*scale;
349  *V10=-TEMP1*scale;
350  if (TEMP1<0.) {
351  *V01= *V10;
352  *V11=-(*V00);
353  } else {
354  *V01=-(*V10);
355  *V11= (*V00);
356  }
357  } else if (TEMP0>0.) {
358  *V00=TEMP0*scale;
359  *V10=TEMP1*scale;
360  if (TEMP1<0.) {
361  *V01=-(*V10);
362  *V11= (*V00);
363  } else {
364  *V01= (*V10);
365  *V11=-(*V00);
366  }
367  } else {
368  *V00=0.;
369  *V10=1;
370  *V11=0.;
371  *V01=1.;
372  }
373  }
374 }
383 inline
385 {
387  if (*V0>0) {
388  s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
389  *V0*=s;
390  *V1*=s;
391  *V2*=s;
392  } else if (*V0<0) {
393  s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
394  *V0*=s;
395  *V1*=s;
396  *V2*=s;
397  } else {
398  if (*V1>0) {
399  s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
400  *V1*=s;
401  *V2*=s;
402  } else if (*V1<0) {
403  s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
404  *V1*=s;
405  *V2*=s;
406  } else {
407  *V2=1.;
408  }
409  }
410 }
437 inline
439  const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22,
444  const DataTypes::real_t tol)
445 {
446  const DataTypes::real_t absA01=fabs(A01);
447  const DataTypes::real_t absA02=fabs(A02);
448  const DataTypes::real_t m=absA01>absA02 ? absA01 : absA02;
449  if (m<=0) {
450  DataTypes::real_t TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
451  eigenvalues_and_eigenvectors2(A11,A12,A22,
452  &TEMP_EV0,&TEMP_EV1,
453  &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
454  if (A00<=TEMP_EV0) {
455  *V00=1.;
456  *V10=0.;
457  *V20=0.;
458  *V01=0.;
459  *V11=TEMP_V00;
460  *V21=TEMP_V10;
461  *V02=0.;
462  *V12=TEMP_V01;
463  *V22=TEMP_V11;
464  *ev0=A00;
465  *ev1=TEMP_EV0;
466  *ev2=TEMP_EV1;
467  } else if (A00>TEMP_EV1) {
468  *V02=1.;
469  *V12=0.;
470  *V22=0.;
471  *V00=0.;
472  *V10=TEMP_V00;
473  *V20=TEMP_V10;
474  *V01=0.;
475  *V11=TEMP_V01;
476  *V21=TEMP_V11;
477  *ev0=TEMP_EV0;
478  *ev1=TEMP_EV1;
479  *ev2=A00;
480  } else {
481  *V01=1.;
482  *V11=0.;
483  *V21=0.;
484  *V00=0.;
485  *V10=TEMP_V00;
486  *V20=TEMP_V10;
487  *V02=0.;
488  *V12=TEMP_V01;
489  *V22=TEMP_V11;
490  *ev0=TEMP_EV0;
491  *ev1=A00;
492  *ev2=TEMP_EV1;
493  }
494  } else {
495  eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
496  const DataTypes::real_t absev0=fabs(*ev0);
497  const DataTypes::real_t absev1=fabs(*ev1);
498  const DataTypes::real_t absev2=fabs(*ev2);
499  DataTypes::real_t max_ev=absev0>absev1 ? absev0 : absev1;
500  max_ev=max_ev>absev2 ? max_ev : absev2;
501  const DataTypes::real_t d_01=fabs((*ev0)-(*ev1));
502  const DataTypes::real_t d_12=fabs((*ev1)-(*ev2));
503  const DataTypes::real_t max_d=d_01>d_12 ? d_01 : d_12;
504  if (max_d<=tol*max_ev) {
505  *V00=1.;
506  *V10=0;
507  *V20=0;
508  *V01=0;
509  *V11=1.;
510  *V21=0;
511  *V02=0;
512  *V12=0;
513  *V22=1.;
514  } else {
515  const DataTypes::real_t S00=A00-(*ev0);
516  const DataTypes::real_t absS00=fabs(S00);
517  if (absS00>m) {
518  vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
519  } else if (absA02<m) {
520  vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
521  } else {
522  vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
523  }
524  normalizeVector3(V00,V10,V20);;
525  const DataTypes::real_t T00=A00-(*ev2);
526  const DataTypes::real_t absT00=fabs(T00);
527  if (absT00>m) {
528  vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
529  } else if (absA02<m) {
530  vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
531  } else {
532  vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
533  }
534  const DataTypes::real_t dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
535  *V02-=dot*(*V00);
536  *V12-=dot*(*V10);
537  *V22-=dot*(*V20);
538  normalizeVector3(V02,V12,V22);
539  *V01=(*V10)*(*V22)-(*V12)*(*V20);
540  *V11=(*V20)*(*V02)-(*V00)*(*V22);
541  *V21=(*V00)*(*V12)-(*V02)*(*V10);
542  normalizeVector3(V01,V11,V21);
543  }
544  }
545 }
546 
547 // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
548 // SM is the product of the last axis_offset entries in arg_0.getShape().
549 template <class LEFT, class RIGHT, class RES>
550 inline
551 void matrix_matrix_product(const int SL, const int SM, const int SR, const LEFT* A, const RIGHT* B, RES* C, int transpose)
552 {
553  if (transpose == 0) {
554  for (int i=0; i<SL; i++) {
555  for (int j=0; j<SR; j++) {
556  RES sum = 0.0;
557  for (int l=0; l<SM; l++) {
558  sum += A[i+SL*l] * B[l+SM*j];
559  }
560  C[i+SL*j] = sum;
561  }
562  }
563  }
564  else if (transpose == 1) {
565  for (int i=0; i<SL; i++) {
566  for (int j=0; j<SR; j++) {
567  RES sum = 0.0;
568  for (int l=0; l<SM; l++) {
569  sum += A[i*SM+l] * B[l+SM*j];
570  }
571  C[i+SL*j] = sum;
572  }
573  }
574  }
575  else if (transpose == 2) {
576  for (int i=0; i<SL; i++) {
577  for (int j=0; j<SR; j++) {
578  RES sum = 0.0;
579  for (int l=0; l<SM; l++) {
580  sum += A[i+SL*l] * B[l*SR+j];
581  }
582  C[i+SL*j] = sum;
583  }
584  }
585  }
586 }
587 
588 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
589 #else
590 
591 inline
593 {
594  return ::erf(x);
595 }
596 
597 inline
599 {
600  return makeNaN();
601 }
602 
603 #endif
604 
606 {
607  return escript::fsign(x);
608 }
609 
611 {
612  return makeNaN();
613 }
614 
615 inline
617 {
618  return acos(x);
619 }
620 
621 inline
623 {
624 #ifdef ESYS_USE_BOOST_ACOS
625  return boost::math::acos(x);
626 #else
627  return acos(x);
628 #endif
629 }
630 
631 
633 {
634  return abs(c);
635 }
636 
637 
638 
639 inline DataTypes::real_t calc_gtzero(const DataTypes::real_t& x) {return x>0;}
641 
642 
643 inline DataTypes::real_t calc_gezero(const DataTypes::real_t& x) {return x>=0;}
645 
646 
647 inline DataTypes::real_t calc_ltzero(const DataTypes::real_t& x) {return x<0;}
649 
650 inline DataTypes::real_t calc_lezero(const DataTypes::real_t& x) {return x<=0;}
652 
653 template <typename IN>
655 {
656  return fabs(i);
657 }
658 
659 template <>
661 {
662  return abs(i);
663 }
664 
665 
666 
667 
668 // deals with unary operations which return real, regardless of
669 // their input type
670 template <class IN>
671 inline void tensor_unary_array_operation_real(const size_t size,
672  const IN *arg1,
673  DataTypes::real_t * argRes,
674  escript::ES_optype operation,
675  DataTypes::real_t tol=0)
676 {
677  switch (operation)
678  {
679  case REAL:
680  for (int i = 0; i < size; ++i) {
681  argRes[i] = std::real(arg1[i]);
682  }
683  break;
684  case IMAG:
685  for (int i = 0; i < size; ++i) {
686  argRes[i] = std::imag(arg1[i]);
687  }
688  break;
689  case EZ:
690  for (size_t i = 0; i < size; ++i) {
691  argRes[i] = (fabs(arg1[i])<=tol);
692  }
693  break;
694  case NEZ:
695  for (size_t i = 0; i < size; ++i) {
696  argRes[i] = (fabs(arg1[i])>tol);
697  }
698  break;
699  case ABS:
700  for (size_t i = 0; i < size; ++i) {
701  argRes[i] = abs_f(arg1[i]);
702  }
703  break;
704  default:
705  throw DataException("Unsupported unary operation");
706  }
707 }
708 
709 
710 
711 template <typename OUT, typename IN>
712 inline OUT conjugate(const IN i)
713 {
714  return conj(i);
715 }
716 
717 // This should never actually be called
718 template <>
720 {
721  return r;
722 }
723 
724 // No openmp because it's called by Lazy
725 // In most cases, IN and OUT will be the same
726 // but not ruling out putting Re() and Im()
727 // through this
728 template <class IN, typename OUT>
729 inline void tensor_unary_array_operation(const size_t size,
730  const IN *arg1,
731  OUT * argRes,
732  escript::ES_optype operation,
733  DataTypes::real_t tol=0)
734 {
735  switch (operation)
736  {
737  case NEG:
738  for (size_t i = 0; i < size; ++i) {
739  argRes[i] = -arg1[i];
740  }
741  break;
742  case SIN:
743  for (size_t i = 0; i < size; ++i) {
744  argRes[i] = sin(arg1[i]);
745  }
746  break;
747  case COS:
748  for (size_t i = 0; i < size; ++i) {
749  argRes[i] = cos(arg1[i]);
750  }
751  break;
752  case TAN:
753  for (size_t i = 0; i < size; ++i) {
754  argRes[i] = tan(arg1[i]);
755  }
756  break;
757  case ASIN:
758  for (size_t i = 0; i < size; ++i) {
759  argRes[i] = asin(arg1[i]);
760  }
761  break;
762  case ACOS:
763  for (size_t i = 0; i < size; ++i) {
764  argRes[i]=calc_acos(arg1[i]);
765  }
766  break;
767  case ATAN:
768  for (size_t i = 0; i < size; ++i) {
769  argRes[i] = atan(arg1[i]);
770  }
771  break;
772  case ABS:
773  for (size_t i = 0; i < size; ++i) {
774  argRes[i] = std::abs(arg1[i]);
775  }
776  break;
777  case SINH:
778  for (size_t i = 0; i < size; ++i) {
779  argRes[i] = sinh(arg1[i]);
780  }
781  break;
782  case COSH:
783  for (size_t i = 0; i < size; ++i) {
784  argRes[i] = cosh(arg1[i]);
785  }
786  break;
787  case TANH:
788  for (size_t i = 0; i < size; ++i) {
789  argRes[i] = tanh(arg1[i]);
790  }
791  break;
792  case ERF:
793  for (size_t i = 0; i < size; ++i) {
794  argRes[i] = calc_erf(arg1[i]);
795  }
796  break;
797  case ASINH:
798  for (size_t i = 0; i < size; ++i) {
799  argRes[i] = asinh(arg1[i]);
800  }
801  break;
802  case ACOSH:
803  for (size_t i = 0; i < size; ++i) {
804  argRes[i] = acosh(arg1[i]);
805  }
806  break;
807  case ATANH:
808  for (size_t i = 0; i < size; ++i) {
809  argRes[i] = atanh(arg1[i]);
810  }
811  break;
812  case LOG10:
813  for (size_t i = 0; i < size; ++i) {
814  argRes[i] = log10(arg1[i]);
815  }
816  break;
817  case LOG:
818  for (size_t i = 0; i < size; ++i) {
819  argRes[i] = log(arg1[i]);
820  }
821  break;
822  case SIGN:
823  for (size_t i = 0; i < size; ++i) {
824  argRes[i] = calc_sign(arg1[i]);
825  }
826  break;
827  case EXP:
828  for (size_t i = 0; i < size; ++i) {
829  argRes[i] = exp(arg1[i]);
830  }
831  break;
832  case SQRT:
833  for (size_t i = 0; i < size; ++i) {
834  argRes[i] = sqrt(arg1[i]);
835  }
836  break;
837  case GZ:
838  for (size_t i = 0; i < size; ++i) {
839  argRes[i] = calc_gtzero(arg1[i]);
840  }
841  break;
842  case GEZ:
843  for (size_t i = 0; i < size; ++i) {
844  argRes[i] = calc_gezero(arg1[i]);
845  }
846  break;
847  case LZ:
848  for (size_t i = 0; i < size; ++i) {
849  argRes[i] = calc_ltzero(arg1[i]);
850  }
851  break;
852  case LEZ:
853  for (size_t i = 0; i < size; ++i) {
854  argRes[i] = calc_lezero(arg1[i]);
855  }
856  break;
857  case CONJ:
858  for (size_t i = 0; i < size; ++i) {
859  argRes[i] = conjugate<OUT,IN>(arg1[i]);
860  }
861  break;
862  case RECIP:
863  for (size_t i = 0; i < size; ++i) {
864  argRes[i] = 1.0/arg1[i];
865  }
866  break;
867  case EZ:
868  for (size_t i = 0; i < size; ++i) {
869  argRes[i] = fabs(arg1[i])<=tol;
870  }
871  break;
872  case NEZ:
873  for (size_t i = 0; i < size; ++i) {
874  argRes[i] = fabs(arg1[i])>tol;
875  }
876  break;
877 
878  default:
879  std::string s="Unsupported unary operation ";
880  s+=operation;
881  throw DataException(s);
882  }
883  return;
884 }
885 
886 bool supports_cplx(escript::ES_optype operation);
887 
888 
889 } // end of namespace
890 
891 #endif // __ESCRIPT_LOCALOPS_H__
892 
Definition: ES_optype.h:45
DataTypes::real_t calc_sign(DataTypes::real_t x)
Definition: ArrayOps.h:605
DataTypes::real_t calc_erf(DataTypes::real_t x)
Definition: ArrayOps.h:592
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:671
Definition: ES_optype.h:38
Definition: ES_optype.h:42
escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
Definition: ArrayOps.h:632
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:148
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:438
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:643
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:647
DataTypes::real_t abs_f(IN i)
Definition: ArrayOps.h:654
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:650
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:170
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:384
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:551
Definition: ES_optype.h:55
Definition: ES_optype.h:46
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:329
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:295
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:237
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:255
DataTypes::real_t calc_gtzero(const DataTypes::real_t &x)
Definition: ArrayOps.h:639
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:129
T second_argument_type
Definition: ArrayOps.h:95
OUT conjugate(const IN i)
Definition: ArrayOps.h:712
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:729
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:64
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:616
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:194