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 type PetscMin(type v1,type v2) 244 245 Not Collective 246 247 Input Parameter: 248 + v1 - first value to find minimum of 249 - v2 - second value to find minimum of 250 251 252 Notes: type can be integer or floating point value 253 254 Level: beginner 255 256 257 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 258 259 M*/ 260 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 261 262 /*MC 263 PetscMax - Returns maxium of two numbers 264 265 Synopsis: 266 type max PetscMax(type v1,type v2) 267 268 Not Collective 269 270 Input Parameter: 271 + v1 - first value to find maximum of 272 - v2 - second value to find maximum of 273 274 Notes: type can be integer or floating point value 275 276 Level: beginner 277 278 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 279 280 M*/ 281 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 282 283 /*MC 284 PetscClipInterval - Returns a number clipped to be within an interval 285 286 Synopsis: 287 type clip PetscClipInterval(type x,type a,type b) 288 289 Not Collective 290 291 Input Parameter: 292 + x - value to use if within interval (a,b) 293 . a - lower end of interval 294 - b - upper end of interval 295 296 Notes: type can be integer or floating point value 297 298 Level: beginner 299 300 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 301 302 M*/ 303 #define PetscClipInterval(x,a,b) (PetscMax((a),PetscMin((x),(b)))) 304 305 /*MC 306 PetscAbsInt - Returns the absolute value of an integer 307 308 Synopsis: 309 int abs PetscAbsInt(int v1) 310 311 Not Collective 312 313 Input Parameter: 314 . v1 - the integer 315 316 Level: beginner 317 318 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 319 320 M*/ 321 #define PetscAbsInt(a) (((a)<0) ? -(a) : (a)) 322 323 /*MC 324 PetscAbsReal - Returns the absolute value of an real number 325 326 Synopsis: 327 Real abs PetscAbsReal(PetscReal v1) 328 329 Not Collective 330 331 Input Parameter: 332 . v1 - the double 333 334 335 Level: beginner 336 337 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 338 339 M*/ 340 #define PetscAbsReal(a) (((a)<0) ? -(a) : (a)) 341 342 /*MC 343 PetscSqr - Returns the square of a number 344 345 Synopsis: 346 type sqr PetscSqr(type v1) 347 348 Not Collective 349 350 Input Parameter: 351 . v1 - the value 352 353 Notes: type can be integer or floating point value 354 355 Level: beginner 356 357 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 358 359 M*/ 360 #define PetscSqr(a) ((a)*(a)) 361 362 /* ----------------------------------------------------------------------------*/ 363 /* 364 Basic constants 365 */ 366 #if defined(PETSC_USE_REAL___FLOAT128) 367 #define PETSC_PI M_PIq 368 #elif defined(M_PI) 369 #define PETSC_PI M_PI 370 #else 371 #define PETSC_PI 3.14159265358979323846264338327950288419716939937510582 372 #endif 373 374 #if !defined(PETSC_USE_64BIT_INDICES) 375 #define PETSC_MAX_INT 2147483647 376 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 377 #else 378 #define PETSC_MAX_INT 9223372036854775807L 379 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 380 #endif 381 382 #if defined(PETSC_USE_REAL_SINGLE) 383 # define PETSC_MAX_REAL 3.40282346638528860e+38F 384 # define PETSC_MIN_REAL -PETSC_MAX_REAL 385 # define PETSC_MACHINE_EPSILON 1.19209290e-07F 386 # define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F 387 # define PETSC_SMALL 1.e-5 388 #elif defined(PETSC_USE_REAL_DOUBLE) 389 # define PETSC_MAX_REAL 1.7976931348623157e+308 390 # define PETSC_MIN_REAL -PETSC_MAX_REAL 391 # define PETSC_MACHINE_EPSILON 2.2204460492503131e-16 392 # define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08 393 # define PETSC_SMALL 1.e-10 394 #elif defined(PETSC_USE_REAL___FLOAT128) 395 # define PETSC_MAX_REAL FLT128_MAX 396 # define PETSC_MIN_REAL -FLT128_MAX 397 # define PETSC_MACHINE_EPSILON FLT128_EPSILON 398 # define PETSC_SQRT_MACHINE_EPSILON 1.38777878078e-17 399 # define PETSC_SMALL 1.e-20 400 #endif 401 402 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanScalar(PetscScalar); 403 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanReal(PetscReal); 404 405 /* ----------------------------------------------------------------------------*/ 406 #define PassiveReal PetscReal 407 #define PassiveScalar PetscScalar 408 409 /* 410 These macros are currently hardwired to match the regular data types, so there is no support for a different 411 MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again. 412 */ 413 #define MPIU_MATSCALAR MPIU_SCALAR 414 typedef PetscScalar MatScalar; 415 typedef PetscReal MatReal; 416 417 struct petsc_mpiu_2scalar {PetscScalar a,b;}; 418 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar); 419 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT) 420 struct petsc_mpiu_2int {PetscInt a,b;}; 421 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int); 422 #else 423 #define MPIU_2INT MPI_2INT 424 #endif 425 426 #endif 427