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