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 334 or double precision. 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 Input Parameter: 345 + v1 - first value to find minimum of 346 - v2 - second value to find minimum of 347 348 Synopsis: 349 type PetscMin(type v1,type v2) 350 351 Notes: type can be integer or floating point value 352 353 Level: beginner 354 355 356 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 357 358 M*/ 359 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 360 361 /*MC 362 PetscMax - Returns maxium of two numbers 363 364 Input Parameter: 365 + v1 - first value to find maximum of 366 - v2 - second value to find maximum of 367 368 Synopsis: 369 type max PetscMax(type v1,type v2) 370 371 Notes: type can be integer or floating point value 372 373 Level: beginner 374 375 .seealso: PetscMin(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 376 377 M*/ 378 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 379 380 /*MC 381 PetscAbsInt - Returns the absolute value of an integer 382 383 Input Parameter: 384 . v1 - the integer 385 386 Synopsis: 387 int abs PetscAbsInt(int v1) 388 389 390 Level: beginner 391 392 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 393 394 M*/ 395 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 396 397 /*MC 398 PetscAbsReal - Returns the absolute value of an real number 399 400 Input Parameter: 401 . v1 - the double 402 403 Synopsis: 404 int abs PetscAbsReal(PetscReal v1) 405 406 407 Level: beginner 408 409 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 410 411 M*/ 412 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 413 414 /*MC 415 PetscSqr - Returns the square of a number 416 417 Input Parameter: 418 . v1 - the value 419 420 Synopsis: 421 type sqr PetscSqr(type v1) 422 423 Notes: type can be integer or floating point value 424 425 Level: beginner 426 427 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 428 429 M*/ 430 #define PetscSqr(a) ((a)*(a)) 431 432 /* ----------------------------------------------------------------------------*/ 433 /* 434 Basic constants - These should be done much better 435 */ 436 #define PETSC_PI 3.14159265358979323846264 437 #define PETSC_DEGREES_TO_RADIANS 0.01745329251994 438 #define PETSC_MAX_INT 2147483647 439 #define PETSC_MIN_INT -2147483647 440 441 #if defined(PETSC_USE_SCALAR_SINGLE) 442 # define PETSC_MAX 1.e30 443 # define PETSC_MIN -1.e30 444 # define PETSC_MACHINE_EPSILON 1.e-7 445 # define PETSC_SQRT_MACHINE_EPSILON 3.e-4 446 # define PETSC_SMALL 1.e-5 447 #elif defined(PETSC_USE_SCALAR_INT) 448 # define PETSC_MAX PETSC_MAX_INT 449 # define PETSC_MIN PETSC_MIN_INT 450 # define PETSC_MACHINE_EPSILON 1 451 # define PETSC_SQRT_MACHINE_EPSILON 1 452 # define PETSC_SMALL 0 453 #elif defined(PETSC_USE_SCALAR_QD_DD) 454 # define PETSC_MAX 1.e300 455 # define PETSC_MIN -1.e300 456 # define PETSC_MACHINE_EPSILON 1.e-30 457 # define PETSC_SQRT_MACHINE_EPSILON 1.e-15 458 # define PETSC_SMALL 1.e-25 459 #else 460 # define PETSC_MAX 1.e300 461 # define PETSC_MIN -1.e300 462 # define PETSC_MACHINE_EPSILON 1.e-14 463 # define PETSC_SQRT_MACHINE_EPSILON 1.e-7 464 # define PETSC_SMALL 1.e-10 465 #endif 466 467 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalMax(PetscReal*,PetscReal*,MPI_Comm); 468 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalMin(PetscReal*,PetscReal*,MPI_Comm); 469 EXTERN PetscErrorCode PETSC_DLLEXPORT PetscGlobalSum(PetscScalar*,PetscScalar*,MPI_Comm); 470 471 /*MC 472 PetscIsInfOrNan - Returns 1 if the input double has an infinity for Not-a-number (Nan) value, otherwise 0. 473 474 Input Parameter: 475 . a - the double 476 477 478 Notes: uses the C99 standard isinf() and isnan() on systems where they exist. 479 Otherwises uses ( (a - a) != 0.0), note that some optimizing compiles compile 480 out this form, thus removing the check. 481 482 Level: beginner 483 484 M*/ 485 #if defined(PETSC_HAVE_ISINF) && defined(PETSC_HAVE_ISNAN) 486 #define PetscIsInfOrNanScalar(a) (isinf(PetscAbsScalar(a)) || isnan(PetscAbsScalar(a))) 487 #define PetscIsInfOrNanReal(a) (isinf(a) || isnan(a)) 488 #elif defined(PETSC_HAVE__FINITE) && defined(PETSC_HAVE__ISNAN) 489 #if defined(PETSC_HAVE_FLOAT_H) 490 #include "float.h" /* windows defines _finite() in float.h */ 491 #endif 492 #if defined(PETSC_HAVE_IEEEFP_H) 493 #include "ieeefp.h" /* Solaris prototypes these here */ 494 #endif 495 #define PetscIsInfOrNanScalar(a) (!_finite(PetscAbsScalar(a)) || _isnan(PetscAbsScalar(a))) 496 #define PetscIsInfOrNanReal(a) (!_finite(a) || _isnan(a)) 497 #else 498 #define PetscIsInfOrNanScalar(a) ((a - a) != 0.0) 499 #define PetscIsInfOrNanReal(a) ((a - a) != 0.0) 500 #endif 501 502 503 /* ----------------------------------------------------------------------------*/ 504 /* 505 PetscLogDouble variables are used to contain double precision numbers 506 that are not used in the numerical computations, but rather in logging, 507 timing etc. 508 */ 509 typedef double PetscLogDouble; 510 #define MPIU_PETSCLOGDOUBLE MPI_DOUBLE 511 512 #define PassiveReal PetscReal 513 #define PassiveScalar PetscScalar 514 515 516 PETSC_EXTERN_CXX_END 517 #endif 518