1 /* 2 3 PETSc mathematics include file. Defines certain basic mathematical 4 constants and functions for working with single and double precision 5 floating point numbers as well as complex and integers. 6 7 This file is included by petscsys.h and should not be used directly. 8 9 */ 10 11 #if !defined(__PETSCMATH_H) 12 #define __PETSCMATH_H 13 #include <math.h> 14 PETSC_EXTERN_CXX_BEGIN 15 16 extern MPI_Datatype MPIU_2SCALAR; 17 extern MPI_Datatype MPIU_2INT; 18 /* 19 20 Defines operations that are different for complex and real numbers; 21 note that one cannot really mix the use of complex and real in the same 22 PETSc program. All PETSc objects in one program are built around the object 23 PetscScalar which is either always a double or a complex. 24 25 */ 26 27 #define PetscExpPassiveScalar(a) PetscExpScalar() 28 29 #if defined(PETSC_USE_COMPLEX) 30 #if defined(PETSC_CLANGUAGE_CXX) 31 /* 32 C++ support of complex numbers: Original support 33 */ 34 #include <complex> 35 36 #if defined(PETSC_USE_SCALAR_SINGLE) 37 /* 38 For d double and c single complex defines the following operations 39 d == c 40 c == d 41 d != c 42 c != d 43 d / c 44 c /d 45 d * c 46 c * d 47 d - c 48 c - d 49 d + c 50 c + d 51 */ 52 namespace std 53 { 54 template<typename _Tp> 55 inline bool 56 operator==(const double& __x, const complex<_Tp>& __y) 57 { return __x == __y.real() && _Tp() == __y.imag(); } 58 template<typename _Tp> 59 inline bool 60 operator==(const complex<_Tp>& __x, const double& __y) 61 { return __x.real() == __y && __x.imag() == _Tp(); } 62 template<typename _Tp> 63 inline bool 64 operator!=(const complex<_Tp>& __x, const double& __y) 65 { return __x.real() != __y || __x.imag() != _Tp(); } 66 template<typename _Tp> 67 inline bool 68 operator!=(const double& __x, const complex<_Tp>& __y) 69 { return __x != __y.real() || _Tp() != __y.imag(); } 70 template<typename _Tp> 71 inline complex<_Tp> 72 operator/(const complex<_Tp>& __x, const double& __y) 73 { 74 complex<_Tp> __r = __x; 75 __r /= ((float)__y); 76 return __r; 77 } 78 template<typename _Tp> 79 inline complex<_Tp> 80 operator/(const double& __x, const complex<_Tp>& __y) 81 { 82 complex<_Tp> __r = (float)__x; 83 __r /= __y; 84 return __r; 85 } 86 template<typename _Tp> 87 inline complex<_Tp> 88 operator*(const complex<_Tp>& __x, const double& __y) 89 { 90 complex<_Tp> __r = __x; 91 __r *= ((float)__y); 92 return __r; 93 } 94 template<typename _Tp> 95 inline complex<_Tp> 96 operator*(const double& __x, const complex<_Tp>& __y) 97 { 98 complex<_Tp> __r = (float)__x; 99 __r *= __y; 100 return __r; 101 } 102 template<typename _Tp> 103 inline complex<_Tp> 104 operator-(const complex<_Tp>& __x, const double& __y) 105 { 106 complex<_Tp> __r = __x; 107 __r -= ((float)__y); 108 return __r; 109 } 110 template<typename _Tp> 111 inline complex<_Tp> 112 operator-(const double& __x, const complex<_Tp>& __y) 113 { 114 complex<_Tp> __r = (float)__x; 115 __r -= __y; 116 return __r; 117 } 118 template<typename _Tp> 119 inline complex<_Tp> 120 operator+(const complex<_Tp>& __x, const double& __y) 121 { 122 complex<_Tp> __r = __x; 123 __r += ((float)__y); 124 return __r; 125 } 126 template<typename _Tp> 127 inline complex<_Tp> 128 operator+(const double& __x, const complex<_Tp>& __y) 129 { 130 complex<_Tp> __r = (float)__x; 131 __r += __y; 132 return __r; 133 } 134 } 135 #endif 136 137 138 139 #define PetscRealPart(a) (a).real() 140 #define PetscImaginaryPart(a) (a).imag() 141 #define PetscAbsScalar(a) std::abs(a) 142 #define PetscConj(a) std::conj(a) 143 #define PetscSqrtScalar(a) std::sqrt(a) 144 #define PetscPowScalar(a,b) std::pow(a,b) 145 #define PetscExpScalar(a) std::exp(a) 146 #define PetscLogScalar(a) std::log(a) 147 #define PetscSinScalar(a) std::sin(a) 148 #define PetscCosScalar(a) std::cos(a) 149 150 #if defined(PETSC_USE_SCALAR_SINGLE) 151 typedef std::complex<float> PetscScalar; 152 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 153 typedef std::complex<long double> PetscScalar; 154 #elif defined(PETSC_USE_SCALAR_INT) 155 typedef std::complex<int> PetscScalar; 156 #else 157 typedef std::complex<double> PetscScalar; 158 #endif 159 #else 160 #include <complex.h> 161 162 /* 163 C support of complex numbers: Warning it needs a 164 C90 compliant compiler to work... 165 */ 166 167 #if defined(PETSC_USE_SCALAR_SINGLE) 168 typedef float complex PetscScalar; 169 170 #define PetscRealPart(a) crealf(a) 171 #define PetscImaginaryPart(a) cimagf(a) 172 #define PetscAbsScalar(a) cabsf(a) 173 #define PetscConj(a) conjf(a) 174 #define PetscSqrtScalar(a) csqrtf(a) 175 #define PetscPowScalar(a,b) cpowf(a,b) 176 #define PetscExpScalar(a) cexpf(a) 177 #define PetscLogScalar(a) clogf(a) 178 #define PetscSinScalar(a) csinf(a) 179 #define PetscCosScalar(a) ccosf(a) 180 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 181 typedef long double complex PetscScalar; 182 183 #define PetscRealPart(a) creall(a) 184 #define PetscImaginaryPart(a) cimagl(a) 185 #define PetscAbsScalar(a) cabsl(a) 186 #define PetscConj(a) conjl(a) 187 #define PetscSqrtScalar(a) csqrtl(a) 188 #define PetscPowScalar(a,b) cpowl(a,b) 189 #define PetscExpScalar(a) cexpl(a) 190 #define PetscLogScalar(a) clogl(a) 191 #define PetscSinScalar(a) csinl(a) 192 #define PetscCosScalar(a) ccosl(a) 193 194 #else 195 typedef double complex PetscScalar; 196 197 #define PetscRealPart(a) creal(a) 198 #define PetscImaginaryPart(a) cimag(a) 199 #define PetscAbsScalar(a) cabs(a) 200 #define PetscConj(a) conj(a) 201 #define PetscSqrtScalar(a) csqrt(a) 202 #define PetscPowScalar(a,b) cpow(a,b) 203 #define PetscExpScalar(a) cexp(a) 204 #define PetscLogScalar(a) clog(a) 205 #define PetscSinScalar(a) csin(a) 206 #define PetscCosScalar(a) ccos(a) 207 #endif 208 #endif 209 210 #if !defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX) 211 extern MPI_Datatype MPI_C_DOUBLE_COMPLEX; 212 extern MPI_Datatype MPI_C_COMPLEX; 213 #endif 214 215 #if defined(PETSC_USE_SCALAR_SINGLE) 216 #define MPIU_SCALAR MPI_C_COMPLEX 217 #else 218 #define MPIU_SCALAR MPI_C_DOUBLE_COMPLEX 219 #endif 220 #if defined(PETSC_USE_SCALAR_MAT_SINGLE) 221 #define MPIU_MATSCALAR ??Notdone 222 #else 223 #define MPIU_MATSCALAR MPI_C_DOUBLE_COMPLEX 224 #endif 225 226 227 /* Compiling for real numbers only */ 228 #else 229 # if defined(PETSC_USE_SCALAR_SINGLE) 230 # define MPIU_SCALAR MPI_FLOAT 231 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 232 # define MPIU_SCALAR MPI_LONG_DOUBLE 233 # elif defined(PETSC_USE_SCALAR_INT) 234 # define MPIU_SCALAR MPI_INT 235 # else 236 # define MPIU_SCALAR MPI_DOUBLE 237 # endif 238 # if defined(PETSC_USE_SCALAR_MAT_SINGLE) || defined(PETSC_USE_SCALAR_SINGLE) 239 # define MPIU_MATSCALAR MPI_FLOAT 240 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 241 # define MPIU_MATSCALAR MPI_LONG_DOUBLE 242 # elif defined(PETSC_USE_SCALAR_INT) 243 # define MPIU_MATSCALAR MPI_INT 244 # else 245 # define MPIU_MATSCALAR MPI_DOUBLE 246 # endif 247 # define PetscRealPart(a) (a) 248 # define PetscImaginaryPart(a) (0.) 249 # define PetscAbsScalar(a) (((a)<0.0) ? -(a) : (a)) 250 # define PetscConj(a) (a) 251 # define PetscSqrtScalar(a) sqrt(a) 252 # define PetscPowScalar(a,b) pow(a,b) 253 # define PetscExpScalar(a) exp(a) 254 # define PetscLogScalar(a) log(a) 255 # define PetscSinScalar(a) sin(a) 256 # define PetscCosScalar(a) cos(a) 257 258 # if defined(PETSC_USE_SCALAR_SINGLE) 259 typedef float PetscScalar; 260 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 261 typedef long double PetscScalar; 262 # elif defined(PETSC_USE_SCALAR_INT) 263 typedef int PetscScalar; 264 # else 265 typedef double PetscScalar; 266 # endif 267 #endif 268 269 #if defined(PETSC_USE_SCALAR_SINGLE) 270 # define MPIU_REAL MPI_FLOAT 271 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 272 # define MPIU_REAL MPI_LONG_DOUBLE 273 #elif defined(PETSC_USE_SCALAR_INT) 274 # define MPIU_REAL MPI_INT 275 #else 276 # define MPIU_REAL MPI_DOUBLE 277 #endif 278 279 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1) 280 #define PetscAbs(a) (((a) >= 0) ? (a) : -(a)) 281 /* 282 Allows compiling PETSc so that matrix values are stored in 283 single precision but all other objects still use double 284 precision. This does not work for complex numbers in that case 285 it remains double 286 287 EXPERIMENTAL! NOT YET COMPLETELY WORKING 288 */ 289 290 #if defined(PETSC_USE_SCALAR_MAT_SINGLE) 291 typedef float MatScalar; 292 #else 293 typedef PetscScalar MatScalar; 294 #endif 295 296 #if defined(PETSC_USE_SCALAR_SINGLE) 297 typedef float PetscReal; 298 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 299 typedef long double PetscReal; 300 #elif defined(PETSC_USE_SCALAR_INT) 301 typedef int PetscReal; 302 #else 303 typedef double PetscReal; 304 #endif 305 306 #if defined(PETSC_USE_COMPLEX) 307 typedef PetscReal MatReal; 308 #elif defined(PETSC_USE_SCALAR_MAT_SINGLE) || defined(PETSC_USE_SCALAR_SINGLE) 309 typedef float MatReal; 310 #else 311 typedef PetscReal MatReal; 312 #endif 313 314 315 /* --------------------------------------------------------------------------*/ 316 317 /* 318 Certain objects may be created using either single or double precision. 319 This is currently not used. 320 */ 321 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE } PetscScalarPrecision; 322 323 /* PETSC_i is the imaginary number, i */ 324 extern PetscScalar PETSC_i; 325 326 /*MC 327 PetscMin - Returns minimum of two numbers 328 329 Synopsis: 330 type PetscMin(type v1,type v2) 331 332 Not Collective 333 334 Input Parameter: 335 + v1 - first value to find minimum of 336 - v2 - second value to find minimum of 337 338 339 Notes: type can be integer or floating point value 340 341 Level: beginner 342 343 344 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 345 346 M*/ 347 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 348 349 /*MC 350 PetscMax - Returns maxium of two numbers 351 352 Synopsis: 353 type max PetscMax(type v1,type v2) 354 355 Not Collective 356 357 Input Parameter: 358 + v1 - first value to find maximum of 359 - v2 - second value to find maximum of 360 361 Notes: type can be integer or floating point value 362 363 Level: beginner 364 365 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 366 367 M*/ 368 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 369 370 /*MC 371 PetscAbsInt - Returns the absolute value of an integer 372 373 Synopsis: 374 int abs PetscAbsInt(int v1) 375 376 Not Collective 377 378 Input Parameter: 379 . v1 - the integer 380 381 Level: beginner 382 383 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 384 385 M*/ 386 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 387 388 /*MC 389 PetscAbsReal - Returns the absolute value of an real number 390 391 Synopsis: 392 Real abs PetscAbsReal(PetscReal v1) 393 394 Not Collective 395 396 Input Parameter: 397 . v1 - the double 398 399 400 Level: beginner 401 402 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 403 404 M*/ 405 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 406 407 /*MC 408 PetscSqr - Returns the square of a number 409 410 Synopsis: 411 type sqr PetscSqr(type v1) 412 413 Not Collective 414 415 Input Parameter: 416 . v1 - the value 417 418 Notes: type can be integer or floating point value 419 420 Level: beginner 421 422 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 423 424 M*/ 425 #define PetscSqr(a) ((a)*(a)) 426 427 /* ----------------------------------------------------------------------------*/ 428 /* 429 Basic constants - These should be done much better 430 */ 431 #define PETSC_PI 3.14159265358979323846264 432 #define PETSC_DEGREES_TO_RADIANS 0.01745329251994 433 #define PETSC_MAX_INT 2147483647 434 #define PETSC_MIN_INT -2147483647 435 436 #if defined(PETSC_USE_SCALAR_SINGLE) 437 # define PETSC_MAX 1.e30 438 # define PETSC_MIN -1.e30 439 # define PETSC_MACHINE_EPSILON 1.e-7 440 # define PETSC_SQRT_MACHINE_EPSILON 3.e-4 441 # define PETSC_SMALL 1.e-5 442 #elif defined(PETSC_USE_SCALAR_INT) 443 # define PETSC_MAX PETSC_MAX_INT 444 # define PETSC_MIN PETSC_MIN_INT 445 # define PETSC_MACHINE_EPSILON 1 446 # define PETSC_SQRT_MACHINE_EPSILON 1 447 # define PETSC_SMALL 0 448 #else 449 # define PETSC_MAX 1.e300 450 # define PETSC_MIN -1.e300 451 # define PETSC_MACHINE_EPSILON 1.e-14 452 # define PETSC_SQRT_MACHINE_EPSILON 1.e-7 453 # define PETSC_SMALL 1.e-10 454 #endif 455 456 #if defined PETSC_HAVE_ADIC 457 /* Use MPI_Allreduce when ADIC is not available. */ 458 extern PetscErrorCode PetscGlobalMax(MPI_Comm, const PetscReal*,PetscReal*); 459 extern PetscErrorCode PetscGlobalMin(MPI_Comm, const PetscReal*,PetscReal*); 460 extern PetscErrorCode PetscGlobalSum(MPI_Comm, const PetscScalar*,PetscScalar*); 461 #endif 462 463 /*MC 464 PetscIsInfOrNan - Returns 1 if the input double has an infinity for Not-a-number (Nan) value, otherwise 0. 465 466 Input Parameter: 467 . a - the double 468 469 470 Notes: uses the C99 standard isinf() and isnan() on systems where they exist. 471 Otherwises uses ( (a - a) != 0.0), note that some optimizing compiles compile 472 out this form, thus removing the check. 473 474 Level: beginner 475 476 M*/ 477 #if defined(PETSC_HAVE_ISINF) && defined(PETSC_HAVE_ISNAN) 478 #define PetscIsInfOrNanScalar(a) (isinf(PetscAbsScalar(a)) || isnan(PetscAbsScalar(a))) 479 #define PetscIsInfOrNanReal(a) (isinf(a) || isnan(a)) 480 #elif defined(PETSC_HAVE__FINITE) && defined(PETSC_HAVE__ISNAN) 481 #if defined(PETSC_HAVE_FLOAT_H) 482 #include "float.h" /* windows defines _finite() in float.h */ 483 #endif 484 #if defined(PETSC_HAVE_IEEEFP_H) 485 #include "ieeefp.h" /* Solaris prototypes these here */ 486 #endif 487 #define PetscIsInfOrNanScalar(a) (!_finite(PetscAbsScalar(a)) || _isnan(PetscAbsScalar(a))) 488 #define PetscIsInfOrNanReal(a) (!_finite(a) || _isnan(a)) 489 #else 490 #define PetscIsInfOrNanScalar(a) ((a - a) != 0.0) 491 #define PetscIsInfOrNanReal(a) ((a - a) != 0.0) 492 #endif 493 494 495 /* ----------------------------------------------------------------------------*/ 496 /* 497 PetscLogDouble variables are used to contain double precision numbers 498 that are not used in the numerical computations, but rather in logging, 499 timing etc. 500 */ 501 typedef double PetscLogDouble; 502 #define MPIU_PETSCLOGDOUBLE MPI_DOUBLE 503 504 #define PassiveReal PetscReal 505 #define PassiveScalar PetscScalar 506 507 508 PETSC_EXTERN_CXX_END 509 #endif 510