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