1 /* 2 3 PETSc mathematics include file. Defines certain basic mathematical 4 constants and functions for working with single, double, and quad precision 5 floating point numbers as well as complex single and double. 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 15 /* 16 17 Defines operations that are different for complex and real numbers; 18 note that one cannot mix the use of complex and real in the same 19 PETSc program. All PETSc objects in one program are built around the object 20 PetscScalar which is either always a real or a complex. 21 22 */ 23 24 #define PetscExpPassiveScalar(a) PetscExpScalar() 25 #if defined(PETSC_USE_REAL_SINGLE) 26 #define MPIU_REAL MPI_FLOAT 27 typedef float PetscReal; 28 #define PetscSqrtReal(a) sqrt(a) 29 #define PetscExpReal(a) exp(a) 30 #define PetscLogReal(a) log(a) 31 #define PetscLog10Real(a) log10(a) 32 #define PetscSinReal(a) sin(a) 33 #define PetscCosReal(a) cos(a) 34 #define PetscTanReal(a) tan(a) 35 #define PetscAsinReal(a) asin(a) 36 #define PetscAcosReal(a) acos(a) 37 #define PetscAtanReal(a) atan(a) 38 #define PetscSinhReal(a) sinh(a) 39 #define PetscCoshReal(a) cosh(a) 40 #define PetscTanhReal(a) tanh(a) 41 #define PetscPowReal(a,b) pow(a,b) 42 #define PetscCeilReal(a) ceil(a) 43 #define PetscFloorReal(a) floor(a) 44 #define PetscFmodReal(a,b) fmod(a,b) 45 #elif defined(PETSC_USE_REAL_DOUBLE) 46 #define MPIU_REAL MPI_DOUBLE 47 typedef double PetscReal; 48 #define PetscSqrtReal(a) sqrt(a) 49 #define PetscExpReal(a) exp(a) 50 #define PetscLogReal(a) log(a) 51 #define PetscLog10Real(a) log10(a) 52 #define PetscSinReal(a) sin(a) 53 #define PetscCosReal(a) cos(a) 54 #define PetscTanReal(a) tan(a) 55 #define PetscSinhReal(a) sinh(a) 56 #define PetscCoshReal(a) cosh(a) 57 #define PetscTanhReal(a) tanh(a) 58 #define PetscAsinReal(a) asin(a) 59 #define PetscAcosReal(a) acos(a) 60 #define PetscAtanReal(a) atan(a) 61 #define PetscPowReal(a,b) pow(a,b) 62 #define PetscCeilReal(a) ceil(a) 63 #define PetscFloorReal(a) floor(a) 64 #define PetscFmodReal(a,b) fmod(a,b) 65 #elif defined(PETSC_USE_REAL___FLOAT128) 66 #if defined(__cplusplus) 67 extern "C" { 68 #endif 69 #include <quadmath.h> 70 #if defined(__cplusplus) 71 } 72 #endif 73 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128 PetscAttrMPITypeTag(__float128); 74 #define MPIU_REAL MPIU___FLOAT128 75 typedef __float128 PetscReal; 76 #define PetscSqrtReal(a) sqrtq(a) 77 #define PetscExpReal(a) expq(a) 78 #define PetscLogReal(a) logq(a) 79 #define PetscLog10Real(a) log10q(a) 80 #define PetscSinReal(a) sinq(a) 81 #define PetscCosReal(a) cosq(a) 82 #define PetscTanReal(a) tanq(a) 83 #define PetscSinhReal(a) sinhq(a) 84 #define PetscCoshReal(a) coshq(a) 85 #define PetscTanhReal(a) tanhq(a) 86 #define PetscAsinReal(a) asinq(a) 87 #define PetscAcosReal(a) acosq(a) 88 #define PetscAtanReal(a) atanq(a) 89 #define PetscAtan2Real(a) atan2q(a) 90 #define PetscPowReal(a,b) powq(a,b) 91 #define PetscCeilReal(a) ceilq(a) 92 #define PetscFloorReal(a) floorq(a) 93 #define PetscFmodReal(a,b) fmodq(a,b) 94 #endif /* PETSC_USE_REAL_* */ 95 96 /* 97 Complex number definitions 98 */ 99 #if defined(PETSC_CLANGUAGE_CXX) && defined(PETSC_HAVE_CXX_COMPLEX) 100 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_DESIRE_COMPLEX) 101 #define PETSC_HAVE_COMPLEX 1 102 /* C++ support of complex number */ 103 #if defined(PETSC_HAVE_CUSP) 104 #define complexlib cusp 105 #include <cusp/complex.h> 106 #else 107 #define complexlib std 108 #include <complex> 109 #endif 110 111 #define PetscRealPartComplex(a) (a).real() 112 #define PetscImaginaryPartComplex(a) (a).imag() 113 #define PetscAbsComplex(a) complexlib::abs(a) 114 #define PetscConjComplex(a) complexlib::conj(a) 115 #define PetscSqrtComplex(a) complexlib::sqrt(a) 116 #define PetscPowComplex(a,b) complexlib::pow(a,b) 117 #define PetscExpComplex(a) complexlib::exp(a) 118 #define PetscLogComplex(a) complexlib::log(a) 119 #define PetscSinComplex(a) complexlib::sin(a) 120 #define PetscCosComplex(a) complexlib::cos(a) 121 #define PetscTanComplex(a) complexlib::tan(a) 122 #define PetscSinhComplex(a) complexlib::sinh(a) 123 #define PetscCoshComplex(a) complexlib::cosh(a) 124 #define PetscTanhComplex(a) complexlib::tanh(a) 125 126 #if defined(PETSC_USE_REAL_SINGLE) 127 typedef complexlib::complex<float> PetscComplex; 128 #elif defined(PETSC_USE_REAL_DOUBLE) 129 typedef complexlib::complex<double> PetscComplex; 130 #elif defined(PETSC_USE_REAL___FLOAT128) 131 typedef complexlib::complex<__float128> PetscComplex; /* Notstandard and not expected to work, use __complex128 */ 132 #endif /* PETSC_USE_REAL_ */ 133 #endif /* PETSC_USE_COMPLEX && PETSC_DESIRE_COMPLEX */ 134 135 #elif defined(PETSC_CLANGUAGE_C) && defined(PETSC_HAVE_C99_COMPLEX) 136 /* Use C99 _Complex for the type. Do not include complex.h by default to define "complex" because of symbol conflicts in Hypre. */ 137 /* Compilation units that can safely use complex should define PETSC_DESIRE_COMPLEX before including any headers */ 138 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_DESIRE_COMPLEX) 139 #define PETSC_HAVE_COMPLEX 1 140 #include <complex.h> 141 142 #if defined(PETSC_USE_REAL_SINGLE) 143 typedef float _Complex PetscComplex; 144 145 #define PetscRealPartComplex(a) crealf(a) 146 #define PetscImaginaryPartComplex(a) cimagf(a) 147 #define PetscAbsComplex(a) cabsf(a) 148 #define PetscConjComplex(a) conjf(a) 149 #define PetscSqrtComplex(a) csqrtf(a) 150 #define PetscPowComplex(a,b) cpowf(a,b) 151 #define PetscExpComplex(a) cexpf(a) 152 #define PetscLogComplex(a) clogf(a) 153 #define PetscSinComplex(a) csinf(a) 154 #define PetscCosComplex(a) ccosf(a) 155 #define PetscTanComplex(a) ctanf(a) 156 #define PetscSinhComplex(a) csinhf(a) 157 #define PetscCoshComplex(a) ccoshf(a) 158 #define PetscTanhComplex(a) ctanhf(a) 159 160 #elif defined(PETSC_USE_REAL_DOUBLE) 161 typedef double _Complex PetscComplex; 162 163 #define PetscRealPartComplex(a) creal(a) 164 #define PetscImaginaryPartComplex(a) cimag(a) 165 #define PetscAbsComplex(a) cabs(a) 166 #define PetscConjComplex(a) conj(a) 167 #define PetscSqrtComplex(a) csqrt(a) 168 #define PetscPowComplex(a,b) cpow(a,b) 169 #define PetscExpComplex(a) cexp(a) 170 #define PetscLogComplex(a) clog(a) 171 #define PetscSinComplex(a) csin(a) 172 #define PetscCosComplex(a) ccos(a) 173 #define PetscTanComplex(a) ctan(a) 174 #define PetscSinhComplex(a) csinh(a) 175 #define PetscCoshComplex(a) ccosh(a) 176 #define PetscTanhComplex(a) ctanh(a) 177 178 #elif defined(PETSC_USE_REAL___FLOAT128) 179 typedef __complex128 PetscComplex; 180 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 PetscAttrMPITypeTag(__complex128); 181 182 #define PetscRealPartComplex(a) crealq(a) 183 #define PetscImaginaryPartComplex(a) cimagq(a) 184 #define PetscAbsComplex(a) cabsq(a) 185 #define PetscConjComplex(a) conjq(a) 186 #define PetscSqrtComplex(a) csqrtq(a) 187 #define PetscPowComplex(a,b) cpowq(a,b) 188 #define PetscExpComplex(a) cexpq(a) 189 #define PetscLogComplex(a) clogq(a) 190 #define PetscSinComplex(a) csinq(a) 191 #define PetscCosComplex(a) ccosq(a) 192 #define PetscTanComplex(a) ctanq(a) 193 #define PetscSinhComplex(a) csinhq(a) 194 #define PetscCoshComplex(a) ccoshq(a) 195 #define PetscTanhComplex(a) ctanhq(a) 196 197 #endif /* PETSC_USE_REAL_* */ 198 #elif defined(PETSC_USE_COMPLEX) 199 #error "PETSc was configured --with-scalar-type=complex, but a language-appropriate complex library is not available" 200 #endif /* PETSC_USE_COMPLEX || PETSC_DESIRE_COMPLEX */ 201 #endif /* (PETSC_CLANGUAGE_CXX && PETSC_HAVE_CXX_COMPLEX) else-if (PETSC_CLANGUAGE_C && PETSC_HAVE_C99_COMPLEX) */ 202 203 #if defined(PETSC_HAVE_COMPLEX) 204 #if defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX) 205 #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX 206 #define MPIU_C_COMPLEX MPI_C_COMPLEX 207 #else 208 # if defined(PETSC_CLANGUAGE_CXX) && defined(PETSC_HAVE_CXX_COMPLEX) 209 typedef complexlib::complex<double> petsc_mpiu_c_double_complex; 210 typedef complexlib::complex<float> petsc_mpiu_c_complex; 211 # elif defined(PETSC_CLANGUAGE_C) && defined(PETSC_HAVE_C99_COMPLEX) 212 typedef double _Complex petsc_mpiu_c_double_complex; 213 typedef float _Complex petsc_mpiu_c_complex; 214 # else 215 typedef struct {double real,imag;} petsc_mpiu_c_double_complex; 216 typedef struct {float real,imag;} petsc_mpiu_c_complex; 217 # endif 218 PETSC_EXTERN MPI_Datatype MPIU_C_DOUBLE_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_double_complex); 219 PETSC_EXTERN MPI_Datatype MPIU_C_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_complex); 220 #endif /* PETSC_HAVE_MPI_C_DOUBLE_COMPLEX */ 221 #endif /* PETSC_HAVE_COMPLEX */ 222 223 #if defined(PETSC_HAVE_COMPLEX) 224 # if defined(PETSC_USE_REAL_SINGLE) 225 # define MPIU_COMPLEX MPIU_C_COMPLEX 226 # elif defined(PETSC_USE_REAL_DOUBLE) 227 # define MPIU_COMPLEX MPIU_C_DOUBLE_COMPLEX 228 # elif defined(PETSC_USE_REAL___FLOAT128) 229 # define MPIU_COMPLEX MPIU___COMPLEX128 230 # endif /* PETSC_USE_REAL_* */ 231 #endif 232 233 #if defined(PETSC_USE_COMPLEX) 234 typedef PetscComplex PetscScalar; 235 #define PetscRealPart(a) PetscRealPartComplex(a) 236 #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a) 237 #define PetscAbsScalar(a) PetscAbsComplex(a) 238 #define PetscConj(a) PetscConjComplex(a) 239 #define PetscSqrtScalar(a) PetscSqrtComplex(a) 240 #define PetscPowScalar(a,b) PetscPowComplex(a,b) 241 #define PetscExpScalar(a) PetscExpComplex(a) 242 #define PetscLogScalar(a) PetscLogComplex(a) 243 #define PetscSinScalar(a) PetscSinComplex(a) 244 #define PetscCosScalar(a) PetscCosComplex(a) 245 #define PetscTanScalar(a) PetscTanComplex(a) 246 #define PetscSinhScalar(a) PetscSinhComplex(a) 247 #define PetscCoshScalar(a) PetscCoshComplex(a) 248 #define PetscTanhScalar(a) PetscTanhComplex(a) 249 #define MPIU_SCALAR MPIU_COMPLEX 250 251 /* 252 real number definitions 253 */ 254 #else /* PETSC_USE_COMPLEX */ 255 typedef PetscReal PetscScalar; 256 #define MPIU_SCALAR MPIU_REAL 257 258 #define PetscRealPart(a) (a) 259 #define PetscImaginaryPart(a) ((PetscReal)0.) 260 PETSC_STATIC_INLINE PetscReal PetscAbsScalar(PetscScalar a) {return a < 0.0 ? -a : a;} 261 #define PetscConj(a) (a) 262 #if !defined(PETSC_USE_REAL___FLOAT128) 263 #define PetscSqrtScalar(a) sqrt(a) 264 #define PetscPowScalar(a,b) pow(a,b) 265 #define PetscExpScalar(a) exp(a) 266 #define PetscLogScalar(a) log(a) 267 #define PetscSinScalar(a) sin(a) 268 #define PetscCosScalar(a) cos(a) 269 #define PetscTanScalar(a) tan(a) 270 #define PetscSinhScalar(a) sinh(a) 271 #define PetscCoshScalar(a) cosh(a) 272 #define PetscTanhScalar(a) tanh(a) 273 #else /* PETSC_USE_REAL___FLOAT128 */ 274 #define PetscSqrtScalar(a) sqrtq(a) 275 #define PetscPowScalar(a,b) powq(a,b) 276 #define PetscExpScalar(a) expq(a) 277 #define PetscLogScalar(a) logq(a) 278 #define PetscSinScalar(a) sinq(a) 279 #define PetscCosScalar(a) cosq(a) 280 #define PetscTanScalar(a) tanq(a) 281 #define PetscSinhScalar(a) sinhq(a) 282 #define PetscCoshScalar(a) coshq(a) 283 #define PetscTanhScalar(a) tanhq(a) 284 #endif /* PETSC_USE_REAL___FLOAT128 */ 285 286 #endif /* PETSC_USE_COMPLEX */ 287 288 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1) 289 #define PetscAbs(a) (((a) >= 0) ? (a) : -(a)) 290 291 /* --------------------------------------------------------------------------*/ 292 293 /* 294 Certain objects may be created using either single or double precision. 295 This is currently not used. 296 */ 297 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE } PetscScalarPrecision; 298 299 #if defined(PETSC_HAVE_COMPLEX) 300 /* PETSC_i is the imaginary number, i */ 301 PETSC_EXTERN PetscComplex PETSC_i; 302 #endif 303 304 /*MC 305 PetscMin - Returns minimum of two numbers 306 307 Synopsis: 308 #include "petscmath.h" 309 type PetscMin(type v1,type v2) 310 311 Not Collective 312 313 Input Parameter: 314 + v1 - first value to find minimum of 315 - v2 - second value to find minimum of 316 317 318 Notes: type can be integer or floating point value 319 320 Level: beginner 321 322 323 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 324 325 M*/ 326 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 327 328 /*MC 329 PetscMax - Returns maxium of two numbers 330 331 Synopsis: 332 #include "petscmath.h" 333 type max PetscMax(type v1,type v2) 334 335 Not Collective 336 337 Input Parameter: 338 + v1 - first value to find maximum of 339 - v2 - second value to find maximum of 340 341 Notes: type can be integer or floating point value 342 343 Level: beginner 344 345 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 346 347 M*/ 348 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 349 350 /*MC 351 PetscClipInterval - Returns a number clipped to be within an interval 352 353 Synopsis: 354 #include "petscmath.h" 355 type clip PetscClipInterval(type x,type a,type b) 356 357 Not Collective 358 359 Input Parameter: 360 + x - value to use if within interval (a,b) 361 . a - lower end of interval 362 - b - upper end of interval 363 364 Notes: type can be integer or floating point value 365 366 Level: beginner 367 368 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 369 370 M*/ 371 #define PetscClipInterval(x,a,b) (PetscMax((a),PetscMin((x),(b)))) 372 373 /*MC 374 PetscAbsInt - Returns the absolute value of an integer 375 376 Synopsis: 377 #include "petscmath.h" 378 int abs PetscAbsInt(int v1) 379 380 Not Collective 381 382 Input Parameter: 383 . v1 - the integer 384 385 Level: beginner 386 387 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 388 389 M*/ 390 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 391 392 /*MC 393 PetscAbsReal - Returns the absolute value of an real number 394 395 Synopsis: 396 #include "petscmath.h" 397 Real abs PetscAbsReal(PetscReal v1) 398 399 Not Collective 400 401 Input Parameter: 402 . v1 - the double 403 404 405 Level: beginner 406 407 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 408 409 M*/ 410 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 411 412 /*MC 413 PetscSqr - Returns the square of a number 414 415 Synopsis: 416 #include "petscmath.h" 417 type sqr PetscSqr(type v1) 418 419 Not Collective 420 421 Input Parameter: 422 . v1 - the value 423 424 Notes: type can be integer or floating point value 425 426 Level: beginner 427 428 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 429 430 M*/ 431 #define PetscSqr(a) ((a)*(a)) 432 433 /* ----------------------------------------------------------------------------*/ 434 /* 435 Basic constants 436 */ 437 #if defined(PETSC_USE_REAL___FLOAT128) 438 #define PETSC_PI M_PIq 439 #elif defined(M_PI) 440 #define PETSC_PI M_PI 441 #else 442 #define PETSC_PI 3.14159265358979323846264338327950288419716939937510582 443 #endif 444 445 #if !defined(PETSC_USE_64BIT_INDICES) 446 #define PETSC_MAX_INT 2147483647 447 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 448 #else 449 #define PETSC_MAX_INT 9223372036854775807L 450 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 451 #endif 452 453 #if defined(PETSC_USE_REAL_SINGLE) 454 # define PETSC_MAX_REAL 3.40282346638528860e+38F 455 # define PETSC_MIN_REAL -PETSC_MAX_REAL 456 # define PETSC_MACHINE_EPSILON 1.19209290e-07F 457 # define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F 458 # define PETSC_SMALL 1.e-5 459 #elif defined(PETSC_USE_REAL_DOUBLE) 460 # define PETSC_MAX_REAL 1.7976931348623157e+308 461 # define PETSC_MIN_REAL -PETSC_MAX_REAL 462 # define PETSC_MACHINE_EPSILON 2.2204460492503131e-16 463 # define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08 464 # define PETSC_SMALL 1.e-10 465 #elif defined(PETSC_USE_REAL___FLOAT128) 466 # define PETSC_MAX_REAL FLT128_MAX 467 # define PETSC_MIN_REAL -FLT128_MAX 468 # define PETSC_MACHINE_EPSILON FLT128_EPSILON 469 # define PETSC_SQRT_MACHINE_EPSILON 1.38777878078e-17 470 # define PETSC_SMALL 1.e-20 471 #endif 472 473 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanScalar(PetscScalar); 474 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanReal(PetscReal); 475 476 /* ----------------------------------------------------------------------------*/ 477 #define PassiveReal PetscReal 478 #define PassiveScalar PetscScalar 479 480 /* 481 These macros are currently hardwired to match the regular data types, so there is no support for a different 482 MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again. 483 */ 484 #define MPIU_MATSCALAR MPIU_SCALAR 485 typedef PetscScalar MatScalar; 486 typedef PetscReal MatReal; 487 488 struct petsc_mpiu_2scalar {PetscScalar a,b;}; 489 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar); 490 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT) 491 struct petsc_mpiu_2int {PetscInt a,b;}; 492 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int); 493 #else 494 #define MPIU_2INT MPI_2INT 495 #endif 496 497 PETSC_STATIC_INLINE PetscInt PetscPowInt(PetscInt base,PetscInt power) { 498 PetscInt result = 1; 499 while (power) { 500 if (power & 1) result *= base; 501 power >>= 1; 502 base *= base; 503 } 504 return result; 505 } 506 PETSC_STATIC_INLINE PetscReal PetscPowRealInt(PetscReal base,PetscInt power) { 507 PetscReal result = 1; 508 while (power) { 509 if (power & 1) result *= base; 510 power >>= 1; 511 base *= base; 512 } 513 return result; 514 } 515 516 #endif 517