/[escript]/trunk/escriptcore/src/LocalOps.h
ViewVC logotype

Contents of /trunk/escriptcore/src/LocalOps.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6000 - (show annotations)
Tue Mar 1 00:24:43 2016 UTC (3 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 30618 byte(s)
a few more include rearrangements.

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 Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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 "UnaryFuncs.h"
24
25 #include <cmath>
26 #include <complex>
27
28 #ifndef M_PI
29 # define M_PI 3.14159265358979323846 /* pi */
30 #endif
31
32
33 /**
34 \file LocalOps.h
35 \brief Describes binary operations performed on double*.
36
37 For operations on DataAbstract see BinaryOp.h.
38 For operations on DataVector see DataMaths.h.
39 */
40
41 namespace escript {
42
43
44 typedef enum
45 {
46 NEGF,
47 SINF,
48 COSF,
49 TANF,
50 ASINF,
51 ACOSF,
52 ATANF,
53 SINHF,
54 COSHF,
55 TANHF,
56 ERFF,
57 ASINHF,
58 ACOSHF,
59 ATANHF,
60 LOG10F,
61 LOGF,
62 SIGNF,
63 ABSF,
64 EXPF,
65 SQRTF,
66 POWF,
67 PLUSF,
68 MINUSF,
69 MULTIPLIESF,
70 DIVIDESF,
71 LESSF,
72 GREATERF,
73 GREATER_EQUALF,
74 LESS_EQUALF,
75 EQZEROF,
76 NEQZEROF,
77 GTZEROF,
78 GEZEROF,
79 LTZEROF,
80 LEZEROF,
81 CONJF,
82 REALF,
83 IMAGF,
84 INVF
85 } ESFunction;
86
87 bool always_real(ESFunction operation);
88
89 /**
90 \brief acts as a wrapper to isnan.
91 \warning if compiler does not support FP_NAN this function will always return false.
92 */
93 inline
94 bool nancheck(double d)
95 {
96 // Q: so why not just test d!=d?
97 // A: Coz it doesn't always work [I've checked].
98 // One theory is that the optimizer skips the test.
99 return std::isnan(d); // isNan should be a function in C++ land
100 }
101
102 /**
103 \brief returns a NaN.
104 \warning Should probably only used where you know you can test for NaNs
105 */
106 inline
107 double makeNaN()
108 {
109 #ifdef nan
110 return nan("");
111 #else
112 return sqrt(-1.);
113 #endif
114
115 }
116
117
118 /**
119 \brief
120 solves a 1x1 eigenvalue A*V=ev*V problem
121
122 \param A00 Input - A_00
123 \param ev0 Output - eigenvalue
124 */
125 inline
126 void eigenvalues1(const double A00,double* ev0) {
127
128 *ev0=A00;
129
130 }
131 /**
132 \brief
133 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
134
135 \param A00 Input - A_00
136 \param A01 Input - A_01
137 \param A11 Input - A_11
138 \param ev0 Output - smallest eigenvalue
139 \param ev1 Output - largest eigenvalue
140 */
141 inline
142 void eigenvalues2(const double A00,const double A01,const double A11,
143 double* ev0, double* ev1) {
144 const double trA=(A00+A11)/2.;
145 const double A_00=A00-trA;
146 const double A_11=A11-trA;
147 const double s=sqrt(A01*A01-A_00*A_11);
148 *ev0=trA-s;
149 *ev1=trA+s;
150 }
151 /**
152 \brief
153 solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
154
155 \param A00 Input - A_00
156 \param A01 Input - A_01
157 \param A02 Input - A_02
158 \param A11 Input - A_11
159 \param A12 Input - A_12
160 \param A22 Input - A_22
161 \param ev0 Output - smallest eigenvalue
162 \param ev1 Output - eigenvalue
163 \param ev2 Output - largest eigenvalue
164 */
165 inline
166 void eigenvalues3(const double A00, const double A01, const double A02,
167 const double A11, const double A12,
168 const double A22,
169 double* ev0, double* ev1,double* ev2) {
170
171 const double trA=(A00+A11+A22)/3.;
172 const double A_00=A00-trA;
173 const double A_11=A11-trA;
174 const double A_22=A22-trA;
175 const double A01_2=A01*A01;
176 const double A02_2=A02*A02;
177 const double A12_2=A12*A12;
178 const double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
179 if (p<=0.) {
180 *ev2=trA;
181 *ev1=trA;
182 *ev0=trA;
183
184 } else {
185 const double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
186 const double sq_p=sqrt(p/3.);
187 double z=-q/(2*pow(sq_p,3));
188 if (z<-1.) {
189 z=-1.;
190 } else if (z>1.) {
191 z=1.;
192 }
193 const double alpha_3=acos(z)/3.;
194 *ev2=trA+2.*sq_p*cos(alpha_3);
195 *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
196 *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
197 }
198 }
199 /**
200 \brief
201 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
202
203 \param A00 Input - A_00
204 \param ev0 Output - eigenvalue
205 \param V00 Output - eigenvector
206 \param tol Input - tolerance to identify to eigenvalues
207 */
208 inline
209 void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
210 {
211 eigenvalues1(A00,ev0);
212 *V00=1.;
213 return;
214 }
215 /**
216 \brief
217 returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
218
219 \param A00 Input - matrix component
220 \param A10 Input - matrix component
221 \param A01 Input - matrix component
222 \param A11 Input - matrix component
223 \param V0 Output - vector component
224 \param V1 Output - vector component
225 */
226 inline
227 void vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
228 double* V0, double*V1)
229 {
230 double absA00=fabs(A00);
231 double absA10=fabs(A10);
232 double absA01=fabs(A01);
233 double absA11=fabs(A11);
234 double m=absA11>absA10 ? absA11 : absA10;
235 if (absA00>m || absA01>m) {
236 *V0=-A01;
237 *V1=A00;
238 } else {
239 if (m<=0) {
240 *V0=1.;
241 *V1=0.;
242 } else {
243 *V0=A11;
244 *V1=-A10;
245 }
246 }
247 }
248 /**
249 \brief
250 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
251 assuming that the kernel dimension is at least 1 and A00 is non zero.
252
253 \param A00 Input - matrix component
254 \param A10 Input - matrix component
255 \param A20 Input - matrix component
256 \param A01 Input - matrix component
257 \param A11 Input - matrix component
258 \param A21 Input - matrix component
259 \param A02 Input - matrix component
260 \param A12 Input - matrix component
261 \param A22 Input - matrix component
262 \param V0 Output - vector component
263 \param V1 Output - vector component
264 \param V2 Output - vector component
265 */
266 inline
267 void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
268 const double A01,const double A11,const double A21,
269 const double A02,const double A12,const double A22,
270 double* V0,double* V1,double* V2)
271 {
272 double TEMP0,TEMP1;
273 const double I00=1./A00;
274 const double IA10=I00*A10;
275 const double IA20=I00*A20;
276 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
277 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
278 *V0=-(A10*TEMP0+A20*TEMP1);
279 *V1=A00*TEMP0;
280 *V2=A00*TEMP1;
281 }
282
283 /**
284 \brief
285 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
286 ordered by increasing value and eigen vectors are normalizeVector3d such that
287 length is zero and first non-zero component is positive.
288
289 \param A00 Input - A_00
290 \param A01 Input - A_01
291 \param A11 Input - A_11
292 \param ev0 Output - smallest eigenvalue
293 \param ev1 Output - eigenvalue
294 \param V00 Output - eigenvector componenent coresponding to ev0
295 \param V10 Output - eigenvector componenent coresponding to ev0
296 \param V01 Output - eigenvector componenent coresponding to ev1
297 \param V11 Output - eigenvector componenent coresponding to ev1
298 \param tol Input - tolerance to identify to eigenvalues
299 */
300 inline
301 void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
302 double* ev0, double* ev1,
303 double* V00, double* V10, double* V01, double* V11,
304 const double tol)
305 {
306 double TEMP0,TEMP1;
307 eigenvalues2(A00,A01,A11,ev0,ev1);
308 const double absev0=fabs(*ev0);
309 const double absev1=fabs(*ev1);
310 double max_ev=absev0>absev1 ? absev0 : absev1;
311 if (fabs((*ev0)-(*ev1))<tol*max_ev) {
312 *V00=1.;
313 *V10=0.;
314 *V01=0.;
315 *V11=1.;
316 } else {
317 vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
318 const double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
319 if (TEMP0<0.) {
320 *V00=-TEMP0*scale;
321 *V10=-TEMP1*scale;
322 if (TEMP1<0.) {
323 *V01= *V10;
324 *V11=-(*V00);
325 } else {
326 *V01=-(*V10);
327 *V11= (*V00);
328 }
329 } else if (TEMP0>0.) {
330 *V00=TEMP0*scale;
331 *V10=TEMP1*scale;
332 if (TEMP1<0.) {
333 *V01=-(*V10);
334 *V11= (*V00);
335 } else {
336 *V01= (*V10);
337 *V11=-(*V00);
338 }
339 } else {
340 *V00=0.;
341 *V10=1;
342 *V11=0.;
343 *V01=1.;
344 }
345 }
346 }
347 /**
348 \brief
349 nomalizes a 3-d vector such that length is one and first non-zero component is positive.
350
351 \param V0 - vector componenent
352 \param V1 - vector componenent
353 \param V2 - vector componenent
354 */
355 inline
356 void normalizeVector3(double* V0,double* V1,double* V2)
357 {
358 double s;
359 if (*V0>0) {
360 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
361 *V0*=s;
362 *V1*=s;
363 *V2*=s;
364 } else if (*V0<0) {
365 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
366 *V0*=s;
367 *V1*=s;
368 *V2*=s;
369 } else {
370 if (*V1>0) {
371 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
372 *V1*=s;
373 *V2*=s;
374 } else if (*V1<0) {
375 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
376 *V1*=s;
377 *V2*=s;
378 } else {
379 *V2=1.;
380 }
381 }
382 }
383 /**
384 \brief
385 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
386 ordered by increasing value and eigen vectors are normalizeVector3d such that
387 length is zero and first non-zero component is positive.
388
389 \param A00 Input - A_00
390 \param A01 Input - A_01
391 \param A02 Input - A_02
392 \param A11 Input - A_11
393 \param A12 Input - A_12
394 \param A22 Input - A_22
395 \param ev0 Output - smallest eigenvalue
396 \param ev1 Output - eigenvalue
397 \param ev2 Output -
398 \param V00 Output - eigenvector componenent coresponding to ev0
399 \param V10 Output - eigenvector componenent coresponding to ev0
400 \param V20 Output -
401 \param V01 Output - eigenvector componenent coresponding to ev1
402 \param V11 Output - eigenvector componenent coresponding to ev1
403 \param V21 Output -
404 \param V02 Output -
405 \param V12 Output -
406 \param V22 Output -
407 \param tol Input - tolerance to identify to eigenvalues
408 */
409 inline
410 void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
411 const double A11, const double A12, const double A22,
412 double* ev0, double* ev1, double* ev2,
413 double* V00, double* V10, double* V20,
414 double* V01, double* V11, double* V21,
415 double* V02, double* V12, double* V22,
416 const double tol)
417 {
418 const double absA01=fabs(A01);
419 const double absA02=fabs(A02);
420 const double m=absA01>absA02 ? absA01 : absA02;
421 if (m<=0) {
422 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
423 eigenvalues_and_eigenvectors2(A11,A12,A22,
424 &TEMP_EV0,&TEMP_EV1,
425 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
426 if (A00<=TEMP_EV0) {
427 *V00=1.;
428 *V10=0.;
429 *V20=0.;
430 *V01=0.;
431 *V11=TEMP_V00;
432 *V21=TEMP_V10;
433 *V02=0.;
434 *V12=TEMP_V01;
435 *V22=TEMP_V11;
436 *ev0=A00;
437 *ev1=TEMP_EV0;
438 *ev2=TEMP_EV1;
439 } else if (A00>TEMP_EV1) {
440 *V02=1.;
441 *V12=0.;
442 *V22=0.;
443 *V00=0.;
444 *V10=TEMP_V00;
445 *V20=TEMP_V10;
446 *V01=0.;
447 *V11=TEMP_V01;
448 *V21=TEMP_V11;
449 *ev0=TEMP_EV0;
450 *ev1=TEMP_EV1;
451 *ev2=A00;
452 } else {
453 *V01=1.;
454 *V11=0.;
455 *V21=0.;
456 *V00=0.;
457 *V10=TEMP_V00;
458 *V20=TEMP_V10;
459 *V02=0.;
460 *V12=TEMP_V01;
461 *V22=TEMP_V11;
462 *ev0=TEMP_EV0;
463 *ev1=A00;
464 *ev2=TEMP_EV1;
465 }
466 } else {
467 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
468 const double absev0=fabs(*ev0);
469 const double absev1=fabs(*ev1);
470 const double absev2=fabs(*ev2);
471 double max_ev=absev0>absev1 ? absev0 : absev1;
472 max_ev=max_ev>absev2 ? max_ev : absev2;
473 const double d_01=fabs((*ev0)-(*ev1));
474 const double d_12=fabs((*ev1)-(*ev2));
475 const double max_d=d_01>d_12 ? d_01 : d_12;
476 if (max_d<=tol*max_ev) {
477 *V00=1.;
478 *V10=0;
479 *V20=0;
480 *V01=0;
481 *V11=1.;
482 *V21=0;
483 *V02=0;
484 *V12=0;
485 *V22=1.;
486 } else {
487 const double S00=A00-(*ev0);
488 const double absS00=fabs(S00);
489 if (absS00>m) {
490 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
491 } else if (absA02<m) {
492 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
493 } else {
494 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
495 }
496 normalizeVector3(V00,V10,V20);;
497 const double T00=A00-(*ev2);
498 const double absT00=fabs(T00);
499 if (absT00>m) {
500 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
501 } else if (absA02<m) {
502 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
503 } else {
504 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
505 }
506 const double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
507 *V02-=dot*(*V00);
508 *V12-=dot*(*V10);
509 *V22-=dot*(*V20);
510 normalizeVector3(V02,V12,V22);
511 *V01=(*V10)*(*V22)-(*V12)*(*V20);
512 *V11=(*V20)*(*V02)-(*V00)*(*V22);
513 *V21=(*V00)*(*V12)-(*V02)*(*V10);
514 normalizeVector3(V01,V11,V21);
515 }
516 }
517 }
518
519 // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
520 // SM is the product of the last axis_offset entries in arg_0.getShape().
521 inline
522 void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
523 {
524 if (transpose == 0) {
525 for (int i=0; i<SL; i++) {
526 for (int j=0; j<SR; j++) {
527 double sum = 0.0;
528 for (int l=0; l<SM; l++) {
529 sum += A[i+SL*l] * B[l+SM*j];
530 }
531 C[i+SL*j] = sum;
532 }
533 }
534 }
535 else if (transpose == 1) {
536 for (int i=0; i<SL; i++) {
537 for (int j=0; j<SR; j++) {
538 double sum = 0.0;
539 for (int l=0; l<SM; l++) {
540 sum += A[i*SM+l] * B[l+SM*j];
541 }
542 C[i+SL*j] = sum;
543 }
544 }
545 }
546 else if (transpose == 2) {
547 for (int i=0; i<SL; i++) {
548 for (int j=0; j<SR; j++) {
549 double sum = 0.0;
550 for (int l=0; l<SM; l++) {
551 sum += A[i+SL*l] * B[l*SR+j];
552 }
553 C[i+SL*j] = sum;
554 }
555 }
556 }
557 }
558
559 template <typename UnaryFunction>
560 inline void tensor_unary_operation(const int size,
561 const double *arg1,
562 double * argRes,
563 UnaryFunction operation)
564 {
565 for (int i = 0; i < size; ++i) {
566 argRes[i] = operation(arg1[i]);
567 }
568 return;
569 }
570
571 // ----------------------
572
573
574 // -------------------------------------
575
576 template <typename BinaryFunction, typename T, typename U, typename V>
577 inline void tensor_binary_operation(const int size,
578 const T *arg1,
579 const U *arg2,
580 V * argRes,
581 BinaryFunction operation)
582 {
583 for (int i = 0; i < size; ++i) {
584 argRes[i] = operation(arg1[i], arg2[i]);
585 }
586 return;
587 }
588
589 template <typename BinaryFunction, typename T, typename U, typename V>
590 inline void tensor_binary_operation(const int size,
591 T arg1,
592 const U *arg2,
593 V *argRes,
594 BinaryFunction operation)
595 {
596 for (int i = 0; i < size; ++i) {
597 argRes[i] = operation(arg1, arg2[i]);
598 }
599 return;
600 }
601
602 template <typename BinaryFunction, typename T, typename U, typename V>
603 inline void tensor_binary_operation(const int size,
604 const T *arg1,
605 U arg2,
606 V *argRes,
607 BinaryFunction operation)
608 {
609 for (int i = 0; i < size; ++i) {
610 argRes[i] = operation(arg1[i], arg2);
611 }
612 return;
613 }
614
615 // following the form of negate from <functional>
616 template <typename T>
617 struct sin_func
618 {
619 T operator() (const T& x) const {return sin(x);}
620 typedef T argument_type;
621 typedef T result_type;
622 };
623
624 template <typename T>
625 struct cos_func
626 {
627 T operator() (const T& x) const {return cos(x);}
628 typedef T argument_type;
629 typedef T result_type;
630 };
631
632 template <typename T>
633 struct tan_func
634 {
635 T operator() (const T& x) const {return tan(x);}
636 typedef T argument_type;
637 typedef T result_type;
638 };
639
640 template <typename T>
641 struct asin_func
642 {
643 T operator() (const T& x) const {return asin(x);}
644 typedef T argument_type;
645 typedef T result_type;
646 };
647
648 template <typename T>
649 struct acos_func
650 {
651 T operator() (const T& x) const {return acos(x);}
652 typedef T argument_type;
653 typedef T result_type;
654 };
655
656 template <typename T>
657 struct atan_func
658 {
659 T operator() (const T& x) const {return atan(x);}
660 typedef T argument_type;
661 typedef T result_type;
662 };
663
664 template <typename T>
665 struct sinh_func
666 {
667 T operator() (const T& x) const {return sinh(x);}
668 typedef T argument_type;
669 typedef T result_type;
670 };
671
672 template <typename T>
673 struct cosh_func
674 {
675 T operator() (const T& x) const {return cosh(x);}
676 typedef T argument_type;
677 typedef T result_type;
678 };
679
680
681 template <typename T>
682 struct tanh_func
683 {
684 T operator() (const T& x) const {return tanh(x);}
685 typedef T argument_type;
686 typedef T result_type;
687 };
688
689 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
690 #else
691 template <typename T>
692 struct erf_func
693 {
694 T operator() (const T& x) const {return ::erf(x);}
695 typedef T argument_type;
696 typedef T result_type;
697 };
698
699 template <>
700 struct erf_func<escript::DataTypes::cplx_t> // dummy instantiation
701 {
702 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
703 typedef DataTypes::cplx_t argument_type;
704 typedef DataTypes::cplx_t result_type;
705 };
706
707 #endif
708
709 template <typename T>
710 struct asinh_func
711 {
712 T operator() (const T& x) const
713 {
714 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
715 return escript::asinh_substitute(x);
716 #else
717 return asinh(x);
718 #endif
719 }
720 typedef T argument_type;
721 typedef T result_type;
722 };
723
724 template <typename T>
725 struct acosh_func
726 {
727 T operator() (const T& x) const
728 {
729 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
730 return escript::acosh_substitute(x);
731 #else
732 return acosh(x);
733 #endif
734 }
735 typedef T argument_type;
736 typedef T result_type;
737 };
738
739 template <typename T>
740 struct atanh_func
741 {
742 T operator() (const T& x) const
743 {
744 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
745 return escript::atanh_substitute(x);
746 #else
747 return atanh(x);
748 #endif
749 }
750 typedef T argument_type;
751 typedef T result_type;
752 };
753
754 template <typename T>
755 struct log10_func
756 {
757 T operator() (const T& x) const {return log10(x);}
758 typedef T argument_type;
759 typedef T result_type;
760 };
761
762 template <typename T>
763 struct log_func
764 {
765 T operator() (const T& x) const {return log(x);}
766 typedef T argument_type;
767 typedef T result_type;
768 };
769
770 template <typename T>
771 struct sign_func
772 {
773 T operator() (const T& x) const {return escript::fsign(x);}
774 typedef T argument_type;
775 typedef T result_type;
776 };
777
778 template <>
779 struct sign_func<DataTypes::cplx_t> // dummy instantiation
780 {
781 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
782 typedef DataTypes::cplx_t argument_type;
783 typedef DataTypes::cplx_t result_type;
784 };
785
786 inline escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
787 {
788 return abs(c);
789 }
790
791
792 template <typename T>
793 struct abs_func
794 {
795 T operator() (const T& x) const {return fabs(x);}
796 typedef T argument_type;
797 typedef T result_type;
798 };
799
800 template <typename T>
801 struct exp_func
802 {
803 T operator() (const T& x) const {return exp(x);}
804 typedef T argument_type;
805 typedef T result_type;
806 };
807
808 template <typename T>
809 struct sqrt_func
810 {
811 T operator() (const T& x) const {return sqrt(x);}
812 typedef T argument_type;
813 typedef T result_type;
814 };
815
816 // following the form of plus from <functional>
817 template <typename T, typename U, typename V>
818 struct pow_func
819 {
820 V operator() (const T& x, const U& y) const {return pow(static_cast<V>(x),static_cast<V>(y));}
821 typedef T first_argument_type;
822 typedef U second_argument_type;
823 typedef V result_type;
824 };
825
826 // following the form of plus from <functional>
827 template <typename T, typename U, typename V>
828 struct plus_func
829 {
830 V operator() (const T& x, const U& y) const {return x+y;}
831 typedef T first_argument_type;
832 typedef U second_argument_type;
833 typedef V result_type;
834 };
835
836 template <typename T, typename U, typename V>
837 struct minus_func
838 {
839 V operator() (const T& x, const U& y) const {return x-y;}
840 typedef T first_argument_type;
841 typedef U second_argument_type;
842 typedef V result_type;
843 };
844
845 template <typename T, typename U, typename V>
846 struct multiplies_func
847 {
848 V operator() (const T& x, const U& y) const {return x*y;}
849 typedef T first_argument_type;
850 typedef U second_argument_type;
851 typedef V result_type;
852 };
853
854 template <typename T, typename U, typename V>
855 struct divides_func
856 {
857 V operator() (const T& x, const U& y) const {return x/y;}
858 typedef T first_argument_type;
859 typedef U second_argument_type;
860 typedef V result_type;
861 };
862
863
864 // using this instead of ::less because that returns bool and we need a result type of T
865 template <typename T>
866 struct less_func
867 {
868 T operator() (const T& x, const T& y) const {return x<y;}
869 typedef T first_argument_type;
870 typedef T second_argument_type;
871 typedef T result_type;
872 };
873
874 // using this instead of ::less because that returns bool and we need a result type of T
875 template <typename T>
876 struct greater_func
877 {
878 T operator() (const T& x, const T& y) const {return x>y;}
879 typedef T first_argument_type;
880 typedef T second_argument_type;
881 typedef T result_type;
882 };
883
884 template <typename T>
885 struct greater_equal_func
886 {
887 T operator() (const T& x, const T& y) const {return x>=y;}
888 typedef T first_argument_type;
889 typedef T second_argument_type;
890 typedef T result_type;
891 };
892
893 template <typename T>
894 struct less_equal_func
895 {
896 T operator() (const T& x, const T& y) const {return x<=y;}
897 typedef T first_argument_type;
898 typedef T second_argument_type;
899 typedef T result_type;
900 };
901
902 template <typename T>
903 struct gtzero_func
904 {
905 T operator() (const T& x) const {return x>0;}
906 typedef T first_argument_type;
907 typedef T result_type;
908 };
909
910 template <>
911 struct gtzero_func<DataTypes::cplx_t> // to keep the templater happy
912 {
913 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
914 typedef DataTypes::cplx_t first_argument_type;
915 typedef DataTypes::cplx_t result_type;
916 };
917
918
919
920 template <typename T>
921 struct gezero_func
922 {
923 T operator() (const T& x) const {return x>=0;}
924 typedef T first_argument_type;
925 typedef T result_type;
926 };
927
928 template <>
929 struct gezero_func<DataTypes::cplx_t> // to keep the templater happy
930 {
931 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
932 typedef DataTypes::cplx_t first_argument_type;
933 typedef DataTypes::cplx_t result_type;
934 };
935
936
937 template <typename T>
938 struct ltzero_func
939 {
940 T operator() (const T& x) const {return x<0;}
941 typedef T first_argument_type;
942 typedef T result_type;
943 };
944
945 template <>
946 struct ltzero_func<DataTypes::cplx_t> // to keep the templater happy
947 {
948 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
949 typedef DataTypes::cplx_t first_argument_type;
950 typedef DataTypes::cplx_t result_type;
951 };
952
953
954
955 template <typename T>
956 struct lezero_func
957 {
958 T operator() (const T& x) const {return x<=0;}
959 typedef T first_argument_type;
960 typedef T result_type;
961 };
962
963 template <>
964 struct lezero_func<DataTypes::cplx_t> // to keep the templater happy
965 {
966 DataTypes::cplx_t operator() (const DataTypes::cplx_t& x) const {return makeNaN();}
967 typedef DataTypes::cplx_t first_argument_type;
968 typedef DataTypes::cplx_t result_type;
969 };
970
971
972 template <class IN, typename OUT, class UnaryFunction>
973 inline void tensor_unary_operation_helper(const size_t size,
974 const IN *arg1,
975 OUT * argRes,
976 UnaryFunction operation)
977 {
978
979 for (int i = 0; i < size; ++i) {
980 argRes[i] = operation(arg1[i]);
981 }
982 }
983
984 template <typename IN>
985 inline DataTypes::real_t abs_f(IN i)
986 {
987 return fabs(i);
988 }
989
990 template <>
991 inline DataTypes::real_t abs_f(DataTypes::cplx_t i)
992 {
993 return abs(i);
994 }
995
996
997
998
999 // deals with unary operations which return real, regardless of
1000 // their input type
1001 template <class IN>
1002 inline void tensor_unary_array_operation_real(const size_t size,
1003 const IN *arg1,
1004 DataTypes::real_t * argRes,
1005 escript::ESFunction operation,
1006 DataTypes::real_t tol=0)
1007 {
1008 switch (operation)
1009 {
1010 case REALF:
1011 for (int i = 0; i < size; ++i) {
1012 argRes[i] = std::real(arg1[i]);
1013 }
1014 break;
1015 case IMAGF:
1016 for (int i = 0; i < size; ++i) {
1017 argRes[i] = std::imag(arg1[i]);
1018 }
1019 break;
1020 case EQZEROF:
1021 for (size_t i = 0; i < size; ++i) {
1022 argRes[i] = (fabs(arg1[i])<=tol);
1023 }
1024 break;
1025 case NEQZEROF:
1026 for (size_t i = 0; i < size; ++i) {
1027 argRes[i] = (fabs(arg1[i])>tol);
1028 }
1029 break;
1030 case ABSF:
1031 for (size_t i = 0; i < size; ++i) {
1032 argRes[i] = abs_f(arg1[i]);
1033 }
1034 break;
1035 default:
1036 throw DataException("Unsupported unary operation");
1037 }
1038 }
1039
1040
1041
1042 template <typename OUT, typename IN>
1043 inline OUT conjugate(const IN i)
1044 {
1045 return conj(i);
1046 }
1047
1048 // This should never actually be called
1049 template <>
1050 inline DataTypes::real_t conjugate(const DataTypes::real_t r)
1051 {
1052 return r;
1053 }
1054
1055 // In most cases, IN and OUT will be the same
1056 // but not ruling out putting Re() and Im()
1057 // through this
1058 template <class IN, typename OUT>
1059 inline void tensor_unary_array_operation(const size_t size,
1060 const IN *arg1,
1061 OUT * argRes,
1062 escript::ESFunction operation,
1063 DataTypes::real_t tol=0)
1064 {
1065 switch (operation)
1066 {
1067 case NEGF:
1068 for (size_t i = 0; i < size; ++i) {
1069 argRes[i] = -arg1[i];
1070 }
1071 break;
1072 case SINF: tensor_unary_operation_helper(size, arg1, argRes, sin_func<IN>()); break;
1073 case COSF: tensor_unary_operation_helper(size, arg1, argRes, cos_func<IN>()); break;
1074 case TANF: tensor_unary_operation_helper(size, arg1, argRes, tan_func<IN>()); break;
1075 case ASINF: tensor_unary_operation_helper(size, arg1, argRes, asin_func<IN>()); break;
1076 case ACOSF: tensor_unary_operation_helper(size, arg1, argRes, acos_func<IN>()); break;
1077 case ATANF: tensor_unary_operation_helper(size, arg1, argRes, atan_func<IN>()); break;
1078 case SINHF: tensor_unary_operation_helper(size, arg1, argRes, sinh_func<IN>()); break;
1079 case COSHF: tensor_unary_operation_helper(size, arg1, argRes, cosh_func<IN>()); break;
1080 case TANHF: tensor_unary_operation_helper(size, arg1, argRes, tanh_func<IN>()); break;
1081 case ERFF: tensor_unary_operation_helper(size, arg1, argRes, erf_func<IN>()); break;
1082 case ASINHF: tensor_unary_operation_helper(size, arg1, argRes, asinh_func<IN>()); break;
1083 case ACOSHF: tensor_unary_operation_helper(size, arg1, argRes, acosh_func<IN>()); break;
1084 case ATANHF: tensor_unary_operation_helper(size, arg1, argRes, atanh_func<IN>()); break;
1085 case LOG10F: tensor_unary_operation_helper(size, arg1, argRes, log10_func<IN>()); break;
1086 case LOGF: tensor_unary_operation_helper(size, arg1, argRes, log_func<IN>()); break;
1087 case SIGNF: tensor_unary_operation_helper(size, arg1, argRes, sign_func<IN>()); break;
1088 case EXPF: tensor_unary_operation_helper(size, arg1, argRes, exp_func<IN>()); break;
1089 case SQRTF: tensor_unary_operation_helper(size, arg1, argRes, sqrt_func<IN>()); break;
1090
1091 case GTZEROF: tensor_unary_operation_helper(size, arg1, argRes, gtzero_func<IN>()); break;
1092 case GEZEROF: tensor_unary_operation_helper(size, arg1, argRes, gezero_func<IN>()); break;
1093 case LTZEROF: tensor_unary_operation_helper(size, arg1, argRes, ltzero_func<IN>()); break;
1094 case LEZEROF: tensor_unary_operation_helper(size, arg1, argRes, lezero_func<IN>()); break;
1095 case CONJF:
1096 for (size_t i = 0; i < size; ++i) {
1097 argRes[i] = conjugate<OUT,IN>(arg1[i]);
1098 }
1099 break;
1100 case INVF:
1101 for (size_t i = 0; i < size; ++i) {
1102 argRes[i] = 1.0/arg1[i];
1103 }
1104 break;
1105
1106 default:
1107 throw DataException("Unsupported unary operation");
1108 }
1109 return;
1110 }
1111
1112 bool supports_cplx(escript::ESFunction operation);
1113
1114
1115 } // end of namespace
1116
1117 #endif // __ESCRIPT_LOCALOPS_H__
1118

  ViewVC Help
Powered by ViewVC 1.1.26