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