1 /* 2 3 PETSc mathematics include file. Defines certain basic mathematical 4 constants and functions for working with single and double precision 5 floating point numbers as well as complex and integers. 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 PETSC_EXTERN_CXX_BEGIN 15 16 extern MPI_Datatype PETSC_DLLEXPORT MPIU_2SCALAR; 17 extern MPI_Datatype PETSC_DLLEXPORT MPIU_2INT; 18 /* 19 20 Defines operations that are different for complex and real numbers; 21 note that one cannot really 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 double or a complex. 24 25 */ 26 27 #define PetscExpPassiveScalar(a) PetscExpScalar() 28 29 #if defined(PETSC_USE_COMPLEX) 30 #if defined(PETSC_CLANGUAGE_CXX) 31 /* 32 C++ support of complex numbers: Original support 33 */ 34 #include <complex> 35 36 #if defined(PETSC_USE_SCALAR_SINGLE) 37 /* 38 For d double and c single complex defines the following operations 39 d == c 40 c == d 41 d != c 42 c != d 43 d / c 44 c /d 45 d * c 46 c * d 47 d - c 48 c - d 49 d + c 50 c + d 51 */ 52 namespace std 53 { 54 template<typename _Tp> 55 inline bool 56 operator==(const double& __x, const complex<_Tp>& __y) 57 { return __x == __y.real() && _Tp() == __y.imag(); } 58 template<typename _Tp> 59 inline bool 60 operator==(const complex<_Tp>& __x, const double& __y) 61 { return __x.real() == __y && __x.imag() == _Tp(); } 62 template<typename _Tp> 63 inline bool 64 operator!=(const complex<_Tp>& __x, const double& __y) 65 { return __x.real() != __y || __x.imag() != _Tp(); } 66 template<typename _Tp> 67 inline bool 68 operator!=(const double& __x, const complex<_Tp>& __y) 69 { return __x != __y.real() || _Tp() != __y.imag(); } 70 template<typename _Tp> 71 inline complex<_Tp> 72 operator/(const complex<_Tp>& __x, const double& __y) 73 { 74 complex<_Tp> __r = __x; 75 __r /= ((float)__y); 76 return __r; 77 } 78 template<typename _Tp> 79 inline complex<_Tp> 80 operator/(const double& __x, const complex<_Tp>& __y) 81 { 82 complex<_Tp> __r = (float)__x; 83 __r /= __y; 84 return __r; 85 } 86 template<typename _Tp> 87 inline complex<_Tp> 88 operator*(const complex<_Tp>& __x, const double& __y) 89 { 90 complex<_Tp> __r = __x; 91 __r *= ((float)__y); 92 return __r; 93 } 94 template<typename _Tp> 95 inline complex<_Tp> 96 operator*(const double& __x, const complex<_Tp>& __y) 97 { 98 complex<_Tp> __r = (float)__x; 99 __r *= __y; 100 return __r; 101 } 102 template<typename _Tp> 103 inline complex<_Tp> 104 operator-(const complex<_Tp>& __x, const double& __y) 105 { 106 complex<_Tp> __r = __x; 107 __r -= ((float)__y); 108 return __r; 109 } 110 template<typename _Tp> 111 inline complex<_Tp> 112 operator-(const double& __x, const complex<_Tp>& __y) 113 { 114 complex<_Tp> __r = (float)__x; 115 __r -= __y; 116 return __r; 117 } 118 template<typename _Tp> 119 inline complex<_Tp> 120 operator+(const complex<_Tp>& __x, const double& __y) 121 { 122 complex<_Tp> __r = __x; 123 __r += ((float)__y); 124 return __r; 125 } 126 template<typename _Tp> 127 inline complex<_Tp> 128 operator+(const double& __x, const complex<_Tp>& __y) 129 { 130 complex<_Tp> __r = (float)__x; 131 __r += __y; 132 return __r; 133 } 134 } 135 #endif 136 137 138 139 #define PetscRealPart(a) (a).real() 140 #define PetscImaginaryPart(a) (a).imag() 141 #define PetscAbsScalar(a) std::abs(a) 142 #define PetscConj(a) std::conj(a) 143 #define PetscSqrtScalar(a) std::sqrt(a) 144 #define PetscPowScalar(a,b) std::pow(a,b) 145 #define PetscExpScalar(a) std::exp(a) 146 #define PetscLogScalar(a) std::log(a) 147 #define PetscSinScalar(a) std::sin(a) 148 #define PetscCosScalar(a) std::cos(a) 149 150 #if defined(PETSC_USE_SCALAR_SINGLE) 151 typedef std::complex<float> PetscScalar; 152 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 153 typedef std::complex<long double> PetscScalar; 154 #elif defined(PETSC_USE_SCALAR_INT) 155 typedef std::complex<int> PetscScalar; 156 #else 157 typedef std::complex<double> PetscScalar; 158 #endif 159 #else 160 #include <complex.h> 161 162 /* 163 C support of complex numbers: Warning it needs a 164 C90 compliant compiler to work... 165 */ 166 167 #if defined(PETSC_USE_SCALAR_SINGLE) 168 typedef float complex PetscScalar; 169 170 #define PetscRealPart(a) crealf(a) 171 #define PetscImaginaryPart(a) cimagf(a) 172 #define PetscAbsScalar(a) cabsf(a) 173 #define PetscConj(a) conjf(a) 174 #define PetscSqrtScalar(a) csqrtf(a) 175 #define PetscPowScalar(a,b) cpowf(a,b) 176 #define PetscExpScalar(a) cexpf(a) 177 #define PetscLogScalar(a) clogf(a) 178 #define PetscSinScalar(a) csinf(a) 179 #define PetscCosScalar(a) ccosf(a) 180 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 181 typedef long double complex PetscScalar; 182 183 #define PetscRealPart(a) creall(a) 184 #define PetscImaginaryPart(a) cimagl(a) 185 #define PetscAbsScalar(a) cabsl(a) 186 #define PetscConj(a) conjl(a) 187 #define PetscSqrtScalar(a) csqrtl(a) 188 #define PetscPowScalar(a,b) cpowl(a,b) 189 #define PetscExpScalar(a) cexpl(a) 190 #define PetscLogScalar(a) clogl(a) 191 #define PetscSinScalar(a) csinl(a) 192 #define PetscCosScalar(a) ccosl(a) 193 194 #else 195 typedef double complex PetscScalar; 196 197 #define PetscRealPart(a) creal(a) 198 #define PetscImaginaryPart(a) cimag(a) 199 #define PetscAbsScalar(a) cabs(a) 200 #define PetscConj(a) conj(a) 201 #define PetscSqrtScalar(a) csqrt(a) 202 #define PetscPowScalar(a,b) cpow(a,b) 203 #define PetscExpScalar(a) cexp(a) 204 #define PetscLogScalar(a) clog(a) 205 #define PetscSinScalar(a) csin(a) 206 #define PetscCosScalar(a) ccos(a) 207 #endif 208 #endif 209 210 #if !defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX) 211 extern MPI_Datatype PETSC_DLLEXPORT MPI_C_DOUBLE_COMPLEX; 212 extern MPI_Datatype PETSC_DLLEXPORT MPI_C_COMPLEX; 213 #endif 214 215 #if defined(PETSC_USE_SCALAR_SINGLE) 216 #define MPIU_SCALAR MPI_C_COMPLEX 217 #else 218 #define MPIU_SCALAR MPI_C_DOUBLE_COMPLEX 219 #endif 220 #if defined(PETSC_USE_SCALAR_MAT_SINGLE) 221 #define MPIU_MATSCALAR ??Notdone 222 #else 223 #define MPIU_MATSCALAR MPI_C_DOUBLE_COMPLEX 224 #endif 225 226 227 /* Compiling for real numbers only */ 228 #else 229 # if defined(PETSC_USE_SCALAR_SINGLE) 230 # define MPIU_SCALAR MPI_FLOAT 231 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 232 # define MPIU_SCALAR MPI_LONG_DOUBLE 233 # elif defined(PETSC_USE_SCALAR_INT) 234 # define MPIU_SCALAR MPI_INT 235 # elif defined(PETSC_USE_SCALAR_QD_DD) 236 # define MPIU_SCALAR MPIU_QD_DD 237 # else 238 # define MPIU_SCALAR MPI_DOUBLE 239 # endif 240 # if defined(PETSC_USE_SCALAR_MAT_SINGLE) || defined(PETSC_USE_SCALAR_SINGLE) 241 # define MPIU_MATSCALAR MPI_FLOAT 242 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 243 # define MPIU_MATSCALAR MPI_LONG_DOUBLE 244 # elif defined(PETSC_USE_SCALAR_INT) 245 # define MPIU_MATSCALAR MPI_INT 246 # elif defined(PETSC_USE_SCALAR_QD_DD) 247 # define MPIU_MATSCALAR MPIU_QD_DD 248 # else 249 # define MPIU_MATSCALAR MPI_DOUBLE 250 # endif 251 # define PetscRealPart(a) (a) 252 # define PetscImaginaryPart(a) (0.) 253 # define PetscAbsScalar(a) (((a)<0.0) ? -(a) : (a)) 254 # define PetscConj(a) (a) 255 # define PetscSqrtScalar(a) sqrt(a) 256 # define PetscPowScalar(a,b) pow(a,b) 257 # define PetscExpScalar(a) exp(a) 258 # define PetscLogScalar(a) log(a) 259 # define PetscSinScalar(a) sin(a) 260 # define PetscCosScalar(a) cos(a) 261 262 # if defined(PETSC_USE_SCALAR_SINGLE) 263 typedef float PetscScalar; 264 # elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 265 typedef long double PetscScalar; 266 # elif defined(PETSC_USE_SCALAR_INT) 267 typedef int PetscScalar; 268 # elif defined(PETSC_USE_SCALAR_QD_DD) 269 # include "qd/dd_real.h" 270 typedef dd_real PetscScalar; 271 # else 272 typedef double PetscScalar; 273 # endif 274 #endif 275 276 #if defined(PETSC_USE_SCALAR_SINGLE) 277 # define MPIU_REAL MPI_FLOAT 278 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 279 # define MPIU_REAL MPI_LONG_DOUBLE 280 #elif defined(PETSC_USE_SCALAR_INT) 281 # define MPIU_REAL MPI_INT 282 #elif defined(PETSC_USE_SCALAR_QD_DD) 283 # define MPIU_REAL MPIU_QD_DD 284 #else 285 # define MPIU_REAL MPI_DOUBLE 286 #endif 287 288 #if defined(PETSC_USE_SCALAR_QD_DD) 289 extern MPI_Datatype PETSC_DLLEXPORT MPIU_QD_DD; 290 #endif 291 292 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1) 293 #define PetscAbs(a) (((a) >= 0) ? (a) : -(a)) 294 /* 295 Allows compiling PETSc so that matrix values are stored in 296 single precision but all other objects still use double 297 precision. This does not work for complex numbers in that case 298 it remains double 299 300 EXPERIMENTAL! NOT YET COMPLETELY WORKING 301 */ 302 303 #if defined(PETSC_USE_SCALAR_MAT_SINGLE) 304 typedef float MatScalar; 305 #else 306 typedef PetscScalar MatScalar; 307 #endif 308 309 #if defined(PETSC_USE_SCALAR_SINGLE) 310 typedef float PetscReal; 311 #elif defined(PETSC_USE_SCALAR_LONG_DOUBLE) 312 typedef long double PetscReal; 313 #elif defined(PETSC_USE_SCALAR_INT) 314 typedef int PetscReal; 315 #elif defined(PETSC_USE_SCALAR_QD_DD) 316 typedef dd_real PetscReal; 317 #else 318 typedef double PetscReal; 319 #endif 320 321 #if defined(PETSC_USE_COMPLEX) 322 typedef PetscReal MatReal; 323 #elif defined(PETSC_USE_SCALAR_MAT_SINGLE) || defined(PETSC_USE_SCALAR_SINGLE) 324 typedef float MatReal; 325 #else 326 typedef PetscReal MatReal; 327 #endif 328 329 330 /* --------------------------------------------------------------------------*/ 331 332 /* 333 Certain objects may be created using either single or double precision. 334 This is currently not used. 335 */ 336 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE, PETSC_SCALAR_QD_DD } PetscScalarPrecision; 337 338 /* PETSC_i is the imaginary number, i */ 339 extern PetscScalar PETSC_DLLEXPORT PETSC_i; 340 341 /*MC 342 PetscMin - Returns minimum of two numbers 343 344 Synopsis: 345 type PetscMin(type v1,type v2) 346 347 Not Collective 348 349 Input Parameter: 350 + v1 - first value to find minimum of 351 - v2 - second value to find minimum of 352 353 354 Notes: type can be integer or floating point value 355 356 Level: beginner 357 358 359 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 360 361 M*/ 362 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 363 364 /*MC 365 PetscMax - Returns maxium of two numbers 366 367 Synopsis: 368 type max PetscMax(type v1,type v2) 369 370 Not Collective 371 372 Input Parameter: 373 + v1 - first value to find maximum of 374 - v2 - second value to find maximum of 375 376 Notes: type can be integer or floating point value 377 378 Level: beginner 379 380 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 381 382 M*/ 383 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 384 385 /*MC 386 PetscAbsInt - Returns the absolute value of an integer 387 388 Synopsis: 389 int abs PetscAbsInt(int v1) 390 391 Not Collective 392 393 Input Parameter: 394 . v1 - the integer 395 396 Level: beginner 397 398 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 399 400 M*/ 401 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 402 403 /*MC 404 PetscAbsReal - Returns the absolute value of an real number 405 406 Synopsis: 407 Real abs PetscAbsReal(PetscReal v1) 408 409 Not Collective 410 411 Input Parameter: 412 . v1 - the double 413 414 415 Level: beginner 416 417 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 418 419 M*/ 420 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 421 422 /*MC 423 PetscSqr - Returns the square of a number 424 425 Synopsis: 426 type sqr PetscSqr(type v1) 427 428 Not Collective 429 430 Input Parameter: 431 . v1 - the value 432 433 Notes: type can be integer or floating point value 434 435 Level: beginner 436 437 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 438 439 M*/ 440 #define PetscSqr(a) ((a)*(a)) 441 442 /* ----------------------------------------------------------------------------*/ 443 /* 444 Basic constants - These should be done much better 445 */ 446 #define PETSC_PI 3.14159265358979323846264 447 #define PETSC_DEGREES_TO_RADIANS 0.01745329251994 448 #define PETSC_MAX_INT 2147483647 449 #define PETSC_MIN_INT -2147483647 450 451 #if defined(PETSC_USE_SCALAR_SINGLE) 452 # define PETSC_MAX 1.e30 453 # define PETSC_MIN -1.e30 454 # define PETSC_MACHINE_EPSILON 1.e-7 455 # define PETSC_SQRT_MACHINE_EPSILON 3.e-4 456 # define PETSC_SMALL 1.e-5 457 #elif defined(PETSC_USE_SCALAR_INT) 458 # define PETSC_MAX PETSC_MAX_INT 459 # define PETSC_MIN PETSC_MIN_INT 460 # define PETSC_MACHINE_EPSILON 1 461 # define PETSC_SQRT_MACHINE_EPSILON 1 462 # define PETSC_SMALL 0 463 #elif defined(PETSC_USE_SCALAR_QD_DD) 464 # define PETSC_MAX 1.e300 465 # define PETSC_MIN -1.e300 466 # define PETSC_MACHINE_EPSILON 1.e-30 467 # define PETSC_SQRT_MACHINE_EPSILON 1.e-15 468 # define PETSC_SMALL 1.e-25 469 #else 470 # define PETSC_MAX 1.e300 471 # define PETSC_MIN -1.e300 472 # define PETSC_MACHINE_EPSILON 1.e-14 473 # define PETSC_SQRT_MACHINE_EPSILON 1.e-7 474 # define PETSC_SMALL 1.e-10 475 #endif 476 477 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalMax(PetscReal*,PetscReal*,MPI_Comm); 478 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalMin(PetscReal*,PetscReal*,MPI_Comm); 479 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalSum(PetscScalar*,PetscScalar*,MPI_Comm); 480 481 /*MC 482 PetscIsInfOrNan - Returns 1 if the input double has an infinity for Not-a-number (Nan) value, otherwise 0. 483 484 Input Parameter: 485 . a - the double 486 487 488 Notes: uses the C99 standard isinf() and isnan() on systems where they exist. 489 Otherwises uses ( (a - a) != 0.0), note that some optimizing compiles compile 490 out this form, thus removing the check. 491 492 Level: beginner 493 494 M*/ 495 #if defined(PETSC_HAVE_ISINF) && defined(PETSC_HAVE_ISNAN) 496 #define PetscIsInfOrNanScalar(a) (isinf(PetscAbsScalar(a)) || isnan(PetscAbsScalar(a))) 497 #define PetscIsInfOrNanReal(a) (isinf(a) || isnan(a)) 498 #elif defined(PETSC_HAVE__FINITE) && defined(PETSC_HAVE__ISNAN) 499 #if defined(PETSC_HAVE_FLOAT_H) 500 #include "float.h" /* windows defines _finite() in float.h */ 501 #endif 502 #if defined(PETSC_HAVE_IEEEFP_H) 503 #include "ieeefp.h" /* Solaris prototypes these here */ 504 #endif 505 #define PetscIsInfOrNanScalar(a) (!_finite(PetscAbsScalar(a)) || _isnan(PetscAbsScalar(a))) 506 #define PetscIsInfOrNanReal(a) (!_finite(a) || _isnan(a)) 507 #else 508 #define PetscIsInfOrNanScalar(a) ((a - a) != 0.0) 509 #define PetscIsInfOrNanReal(a) ((a - a) != 0.0) 510 #endif 511 512 513 /* ----------------------------------------------------------------------------*/ 514 /* 515 PetscLogDouble variables are used to contain double precision numbers 516 that are not used in the numerical computations, but rather in logging, 517 timing etc. 518 */ 519 typedef double PetscLogDouble; 520 #define MPIU_PETSCLOGDOUBLE MPI_DOUBLE 521 522 #define PassiveReal PetscReal 523 #define PassiveScalar PetscScalar 524 525 526 PETSC_EXTERN_CXX_END 527 #endif 528