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