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