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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6515 - (show annotations)
Mon Mar 6 01:54:53 2017 UTC (16 months, 1 week ago) by jfenwick
File MIME type: text/plain
File size: 25933 byte(s)
More helpful exception message

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
36 /**
37 \file LocalOps.h
38 \brief Describes binary operations performed on double*.
39
40 For operations on DataAbstract see BinaryOp.h.
41 For operations on DataVector see DataMaths.h.
42 */
43
44 #include "ES_optype.h"
45
46 namespace escript {
47
48
49
50 bool always_real(escript::ES_optype operation);
51
52
53 /**
54 \brief
55 Return the maximum value of the two given values.
56 */
57 struct FMax
58 {
59 inline DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
60 {
61 return std::max(x,y);
62 }
63 typedef DataTypes::real_t first_argument_type;
64 typedef DataTypes::real_t second_argument_type;
65 typedef DataTypes::real_t result_type;
66 };
67
68 /**
69 \brief
70 Return the minimum value of the two given values.
71 */
72 struct FMin
73 {
74 inline DataTypes::real_t operator()(DataTypes::real_t x, DataTypes::real_t y) const
75 {
76 return std::min(x,y);
77 }
78 typedef DataTypes::real_t first_argument_type;
79 typedef DataTypes::real_t second_argument_type;
80 typedef DataTypes::real_t result_type;
81 };
82
83 /**
84 \brief
85 Return the absolute maximum value of the two given values.
86 */
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 }
94 typedef T first_argument_type;
95 typedef T second_argument_type;
96 typedef DataTypes::real_t result_type;
97 };
98
99
100 inline
101 DataTypes::real_t
102 fsign(DataTypes::real_t x)
103 {
104 if (x == 0) {
105 return 0;
106 } else {
107 return x/fabs(x);
108 }
109 }
110
111 /**
112 \brief acts as a wrapper to isnan.
113 \warning if compiler does not support FP_NAN this function will always return false.
114 */
115 inline
116 bool nancheck(DataTypes::real_t d)
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 /**
125 \brief returns a NaN.
126 \warning Should probably only used where you know you can test for NaNs
127 */
128 inline
129 DataTypes::real_t makeNaN()
130 {
131 #ifdef nan
132 return nan("");
133 #else
134 return sqrt(-1.);
135 #endif
136
137 }
138
139
140 /**
141 \brief
142 solves a 1x1 eigenvalue A*V=ev*V problem
143
144 \param A00 Input - A_00
145 \param ev0 Output - eigenvalue
146 */
147 inline
148 void eigenvalues1(const DataTypes::real_t A00,DataTypes::real_t* ev0) {
149 *ev0=A00;
150 }
151
152 inline
153 void eigenvalues1(const DataTypes::cplx_t A00,DataTypes::cplx_t* ev0) {
154 *ev0=A00;
155 }
156
157
158 /**
159 \brief
160 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
161
162 \param A00 Input - A_00
163 \param A01 Input - A_01
164 \param A11 Input - A_11
165 \param ev0 Output - smallest eigenvalue
166 \param ev1 Output - largest eigenvalue
167 */
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 }
179 /**
180 \brief
181 solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
182
183 \param A00 Input - A_00
184 \param A01 Input - A_01
185 \param A02 Input - A_02
186 \param A11 Input - A_11
187 \param A12 Input - A_12
188 \param A22 Input - A_22
189 \param ev0 Output - smallest eigenvalue
190 \param ev1 Output - eigenvalue
191 \param ev2 Output - largest eigenvalue
192 */
193 inline
194 void eigenvalues3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02,
195 const DataTypes::real_t A11, const DataTypes::real_t A12,
196 const DataTypes::real_t A22,
197 DataTypes::real_t* ev0, DataTypes::real_t* ev1,DataTypes::real_t* ev2) {
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 }
227 /**
228 \brief
229 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
230
231 \param A00 Input - A_00
232 \param ev0 Output - eigenvalue
233 \param V00 Output - eigenvector
234 \param tol Input - tolerance to identify to eigenvalues
235 */
236 inline
237 void eigenvalues_and_eigenvectors1(const DataTypes::real_t A00,DataTypes::real_t* ev0,DataTypes::real_t* V00,const DataTypes::real_t tol)
238 {
239 eigenvalues1(A00,ev0);
240 *V00=1.;
241 return;
242 }
243 /**
244 \brief
245 returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
246
247 \param A00 Input - matrix component
248 \param A10 Input - matrix component
249 \param A01 Input - matrix component
250 \param A11 Input - matrix component
251 \param V0 Output - vector component
252 \param V1 Output - vector component
253 */
254 inline
255 void vectorInKernel2(const DataTypes::real_t A00,const DataTypes::real_t A10,const DataTypes::real_t A01,const DataTypes::real_t A11,
256 DataTypes::real_t* V0, DataTypes::real_t*V1)
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 }
276 /**
277 \brief
278 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
279 assuming that the kernel dimension is at least 1 and A00 is non zero.
280
281 \param A00 Input - matrix component
282 \param A10 Input - matrix component
283 \param A20 Input - matrix component
284 \param A01 Input - matrix component
285 \param A11 Input - matrix component
286 \param A21 Input - matrix component
287 \param A02 Input - matrix component
288 \param A12 Input - matrix component
289 \param A22 Input - matrix component
290 \param V0 Output - vector component
291 \param V1 Output - vector component
292 \param V2 Output - vector component
293 */
294 inline
295 void vectorInKernel3__nonZeroA00(const DataTypes::real_t A00,const DataTypes::real_t A10,const DataTypes::real_t A20,
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,
298 DataTypes::real_t* V0,DataTypes::real_t* V1,DataTypes::real_t* V2)
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
311 /**
312 \brief
313 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
314 ordered by increasing value and eigen vectors are normalizeVector3d such that
315 length is zero and first non-zero component is positive.
316
317 \param A00 Input - A_00
318 \param A01 Input - A_01
319 \param A11 Input - A_11
320 \param ev0 Output - smallest eigenvalue
321 \param ev1 Output - eigenvalue
322 \param V00 Output - eigenvector componenent coresponding to ev0
323 \param V10 Output - eigenvector componenent coresponding to ev0
324 \param V01 Output - eigenvector componenent coresponding to ev1
325 \param V11 Output - eigenvector componenent coresponding to ev1
326 \param tol Input - tolerance to identify to eigenvalues
327 */
328 inline
329 void eigenvalues_and_eigenvectors2(const DataTypes::real_t A00,const DataTypes::real_t A01,const DataTypes::real_t A11,
330 DataTypes::real_t* ev0, DataTypes::real_t* ev1,
331 DataTypes::real_t* V00, DataTypes::real_t* V10, DataTypes::real_t* V01, DataTypes::real_t* V11,
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 }
375 /**
376 \brief
377 nomalizes a 3-d vector such that length is one and first non-zero component is positive.
378
379 \param V0 - vector componenent
380 \param V1 - vector componenent
381 \param V2 - vector componenent
382 */
383 inline
384 void normalizeVector3(DataTypes::real_t* V0,DataTypes::real_t* V1,DataTypes::real_t* V2)
385 {
386 DataTypes::real_t s;
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 }
411 /**
412 \brief
413 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
414 ordered by increasing value and eigen vectors are normalizeVector3d such that
415 length is zero and first non-zero component is positive.
416
417 \param A00 Input - A_00
418 \param A01 Input - A_01
419 \param A02 Input - A_02
420 \param A11 Input - A_11
421 \param A12 Input - A_12
422 \param A22 Input - A_22
423 \param ev0 Output - smallest eigenvalue
424 \param ev1 Output - eigenvalue
425 \param ev2 Output -
426 \param V00 Output - eigenvector componenent coresponding to ev0
427 \param V10 Output - eigenvector componenent coresponding to ev0
428 \param V20 Output -
429 \param V01 Output - eigenvector componenent coresponding to ev1
430 \param V11 Output - eigenvector componenent coresponding to ev1
431 \param V21 Output -
432 \param V02 Output -
433 \param V12 Output -
434 \param V22 Output -
435 \param tol Input - tolerance to identify to eigenvalues
436 */
437 inline
438 void eigenvalues_and_eigenvectors3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02,
439 const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22,
440 DataTypes::real_t* ev0, DataTypes::real_t* ev1, DataTypes::real_t* ev2,
441 DataTypes::real_t* V00, DataTypes::real_t* V10, DataTypes::real_t* V20,
442 DataTypes::real_t* V01, DataTypes::real_t* V11, DataTypes::real_t* V21,
443 DataTypes::real_t* V02, DataTypes::real_t* V12, DataTypes::real_t* V22,
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
592 DataTypes::real_t calc_erf(DataTypes::real_t x)
593 {
594 return ::erf(x);
595 }
596
597 inline
598 DataTypes::cplx_t calc_erf(DataTypes::cplx_t x)
599 {
600 return makeNaN();
601 }
602
603 #endif
604
605 inline DataTypes::real_t calc_sign(DataTypes::real_t x)
606 {
607 return escript::fsign(x);
608 }
609
610 inline DataTypes::cplx_t calc_sign(DataTypes::cplx_t x)
611 {
612 return makeNaN();
613 }
614
615 inline
616 DataTypes::real_t calc_acos(DataTypes::real_t x)
617 {
618 return acos(x);
619 }
620
621 inline
622 DataTypes::cplx_t calc_acos(DataTypes::cplx_t x)
623 {
624 #ifdef ESYS_USE_BOOST_ACOS
625 return boost::math::acos(x);
626 #else
627 return acos(x);
628 #endif
629 }
630
631
632 inline escript::DataTypes::real_t fabs(const escript::DataTypes::cplx_t c)
633 {
634 return abs(c);
635 }
636
637
638
639 inline DataTypes::real_t calc_gtzero(const DataTypes::real_t& x) {return x>0;}
640 inline DataTypes::cplx_t calc_gtzero(const DataTypes::cplx_t& x) {return makeNaN();}
641
642
643 inline DataTypes::real_t calc_gezero(const DataTypes::real_t& x) {return x>=0;}
644 inline DataTypes::cplx_t calc_gezero(const DataTypes::cplx_t& x) {return makeNaN();}
645
646
647 inline DataTypes::real_t calc_ltzero(const DataTypes::real_t& x) {return x<0;}
648 inline DataTypes::cplx_t calc_ltzero(const DataTypes::cplx_t& x) {return makeNaN();}
649
650 inline DataTypes::real_t calc_lezero(const DataTypes::real_t& x) {return x<=0;}
651 inline DataTypes::cplx_t calc_lezero(const DataTypes::cplx_t& x) {return makeNaN();}
652
653 template <typename IN>
654 inline DataTypes::real_t abs_f(IN i)
655 {
656 return fabs(i);
657 }
658
659 template <>
660 inline DataTypes::real_t abs_f(DataTypes::cplx_t i)
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 case PHS:
705 for (size_t i = 0; i < size; ++i) {
706 argRes[i] = std::arg(arg1[i]);
707 }
708 break;
709 default:
710 std::ostringstream oss;
711 oss << "Unsupported unary operation=";
712 oss << opToString(operation);
713 oss << '/';
714 oss << operation;
715 oss << " (Was expecting an operation with real results)";
716 throw DataException(oss.str());
717 }
718 }
719
720
721
722 template <typename OUT, typename IN>
723 inline OUT conjugate(const IN i)
724 {
725 return conj(i);
726 }
727
728 // This should never actually be called
729 template <>
730 inline DataTypes::real_t conjugate(const DataTypes::real_t r)
731 {
732 return r;
733 }
734
735 // No openmp because it's called by Lazy
736 // In most cases, IN and OUT will be the same
737 // but not ruling out putting Re() and Im()
738 // through this
739 template <class IN, typename OUT>
740 inline void tensor_unary_array_operation(const size_t size,
741 const IN *arg1,
742 OUT * argRes,
743 escript::ES_optype operation,
744 DataTypes::real_t tol=0)
745 {
746 switch (operation)
747 {
748 case NEG:
749 for (size_t i = 0; i < size; ++i) {
750 argRes[i] = -arg1[i];
751 }
752 break;
753 case SIN:
754 for (size_t i = 0; i < size; ++i) {
755 argRes[i] = sin(arg1[i]);
756 }
757 break;
758 case COS:
759 for (size_t i = 0; i < size; ++i) {
760 argRes[i] = cos(arg1[i]);
761 }
762 break;
763 case TAN:
764 for (size_t i = 0; i < size; ++i) {
765 argRes[i] = tan(arg1[i]);
766 }
767 break;
768 case ASIN:
769 for (size_t i = 0; i < size; ++i) {
770 argRes[i] = asin(arg1[i]);
771 }
772 break;
773 case ACOS:
774 for (size_t i = 0; i < size; ++i) {
775 argRes[i]=calc_acos(arg1[i]);
776 }
777 break;
778 case ATAN:
779 for (size_t i = 0; i < size; ++i) {
780 argRes[i] = atan(arg1[i]);
781 }
782 break;
783 case ABS:
784 for (size_t i = 0; i < size; ++i) {
785 argRes[i] = std::abs(arg1[i]);
786 }
787 break;
788 case SINH:
789 for (size_t i = 0; i < size; ++i) {
790 argRes[i] = sinh(arg1[i]);
791 }
792 break;
793 case COSH:
794 for (size_t i = 0; i < size; ++i) {
795 argRes[i] = cosh(arg1[i]);
796 }
797 break;
798 case TANH:
799 for (size_t i = 0; i < size; ++i) {
800 argRes[i] = tanh(arg1[i]);
801 }
802 break;
803 case ERF:
804 for (size_t i = 0; i < size; ++i) {
805 argRes[i] = calc_erf(arg1[i]);
806 }
807 break;
808 case ASINH:
809 for (size_t i = 0; i < size; ++i) {
810 argRes[i] = asinh(arg1[i]);
811 }
812 break;
813 case ACOSH:
814 for (size_t i = 0; i < size; ++i) {
815 argRes[i] = acosh(arg1[i]);
816 }
817 break;
818 case ATANH:
819 for (size_t i = 0; i < size; ++i) {
820 argRes[i] = atanh(arg1[i]);
821 }
822 break;
823 case LOG10:
824 for (size_t i = 0; i < size; ++i) {
825 argRes[i] = log10(arg1[i]);
826 }
827 break;
828 case LOG:
829 for (size_t i = 0; i < size; ++i) {
830 argRes[i] = log(arg1[i]);
831 }
832 break;
833 case SIGN:
834 for (size_t i = 0; i < size; ++i) {
835 argRes[i] = calc_sign(arg1[i]);
836 }
837 break;
838 case EXP:
839 for (size_t i = 0; i < size; ++i) {
840 argRes[i] = exp(arg1[i]);
841 }
842 break;
843 case SQRT:
844 for (size_t i = 0; i < size; ++i) {
845 argRes[i] = sqrt(arg1[i]);
846 }
847 break;
848 case GZ:
849 for (size_t i = 0; i < size; ++i) {
850 argRes[i] = calc_gtzero(arg1[i]);
851 }
852 break;
853 case GEZ:
854 for (size_t i = 0; i < size; ++i) {
855 argRes[i] = calc_gezero(arg1[i]);
856 }
857 break;
858 case LZ:
859 for (size_t i = 0; i < size; ++i) {
860 argRes[i] = calc_ltzero(arg1[i]);
861 }
862 break;
863 case LEZ:
864 for (size_t i = 0; i < size; ++i) {
865 argRes[i] = calc_lezero(arg1[i]);
866 }
867 break;
868 case CONJ:
869 for (size_t i = 0; i < size; ++i) {
870 argRes[i] = conjugate<OUT,IN>(arg1[i]);
871 }
872 break;
873 case RECIP:
874 for (size_t i = 0; i < size; ++i) {
875 argRes[i] = 1.0/arg1[i];
876 }
877 break;
878 case EZ:
879 for (size_t i = 0; i < size; ++i) {
880 argRes[i] = fabs(arg1[i])<=tol;
881 }
882 break;
883 case NEZ:
884 for (size_t i = 0; i < size; ++i) {
885 argRes[i] = fabs(arg1[i])>tol;
886 }
887 break;
888
889 default:
890 std::ostringstream oss;
891 oss << "Unsupported unary operation=";
892 oss << opToString(operation);
893 oss << '/';
894 oss << operation;
895 throw DataException(oss.str());
896 }
897 return;
898 }
899
900 bool supports_cplx(escript::ES_optype operation);
901
902
903 } // end of namespace
904
905 #endif // __ESCRIPT_LOCALOPS_H__
906

  ViewVC Help
Powered by ViewVC 1.1.26