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