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 #if defined(PETSC_USE_REAL_SINGLE) 25 #define MPIU_REAL MPI_FLOAT 26 typedef float PetscReal; 27 #define PetscSqrtReal(a) sqrt(a) 28 #define PetscExpReal(a) exp(a) 29 #define PetscLogReal(a) log(a) 30 #define PetscLog10Real(a) log10(a) 31 #ifdef PETSC_HAVE_LOG2 32 #define PetscLog2Real(a) log2(a) 33 #endif 34 #define PetscSinReal(a) sin(a) 35 #define PetscCosReal(a) cos(a) 36 #define PetscTanReal(a) tan(a) 37 #define PetscAsinReal(a) asin(a) 38 #define PetscAcosReal(a) acos(a) 39 #define PetscAtanReal(a) atan(a) 40 #define PetscAtan2Real(a,b) atan2(a,b) 41 #define PetscSinhReal(a) sinh(a) 42 #define PetscCoshReal(a) cosh(a) 43 #define PetscTanhReal(a) tanh(a) 44 #define PetscPowReal(a,b) pow(a,b) 45 #define PetscCeilReal(a) ceil(a) 46 #define PetscFloorReal(a) floor(a) 47 #define PetscFmodReal(a,b) fmod(a,b) 48 #define PetscTGamma(a) tgammaf(a) 49 #elif defined(PETSC_USE_REAL_DOUBLE) 50 #define MPIU_REAL MPI_DOUBLE 51 typedef double PetscReal; 52 #define PetscSqrtReal(a) sqrt(a) 53 #define PetscExpReal(a) exp(a) 54 #define PetscLogReal(a) log(a) 55 #define PetscLog10Real(a) log10(a) 56 #ifdef PETSC_HAVE_LOG2 57 #define PetscLog2Real(a) log2(a) 58 #endif 59 #define PetscSinReal(a) sin(a) 60 #define PetscCosReal(a) cos(a) 61 #define PetscTanReal(a) tan(a) 62 #define PetscAsinReal(a) asin(a) 63 #define PetscAcosReal(a) acos(a) 64 #define PetscAtanReal(a) atan(a) 65 #define PetscAtan2Real(a,b) atan2(a,b) 66 #define PetscSinhReal(a) sinh(a) 67 #define PetscCoshReal(a) cosh(a) 68 #define PetscTanhReal(a) tanh(a) 69 #define PetscPowReal(a,b) pow(a,b) 70 #define PetscCeilReal(a) ceil(a) 71 #define PetscFloorReal(a) floor(a) 72 #define PetscFmodReal(a,b) fmod(a,b) 73 #define PetscTGamma(a) tgamma(a) 74 #elif defined(PETSC_USE_REAL___FLOAT128) 75 #if defined(__cplusplus) 76 extern "C" { 77 #endif 78 #include <quadmath.h> 79 #if defined(__cplusplus) 80 } 81 #endif 82 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128 PetscAttrMPITypeTag(__float128); 83 #define MPIU_REAL MPIU___FLOAT128 84 typedef __float128 PetscReal; 85 #define PetscSqrtReal(a) sqrtq(a) 86 #define PetscExpReal(a) expq(a) 87 #define PetscLogReal(a) logq(a) 88 #define PetscLog10Real(a) log10q(a) 89 #ifdef PETSC_HAVE_LOG2 90 #define PetscLog2Real(a) log2q(a) 91 #endif 92 #define PetscSinReal(a) sinq(a) 93 #define PetscCosReal(a) cosq(a) 94 #define PetscTanReal(a) tanq(a) 95 #define PetscAsinReal(a) asinq(a) 96 #define PetscAcosReal(a) acosq(a) 97 #define PetscAtanReal(a) atanq(a) 98 #define PetscAtan2Real(a,b) atan2q(a,b) 99 #define PetscSinhReal(a) sinhq(a) 100 #define PetscCoshReal(a) coshq(a) 101 #define PetscTanhReal(a) tanhq(a) 102 #define PetscPowReal(a,b) powq(a,b) 103 #define PetscCeilReal(a) ceilq(a) 104 #define PetscFloorReal(a) floorq(a) 105 #define PetscFmodReal(a,b) fmodq(a,b) 106 #define PetscTGamma(a) tgammaq(a) 107 #elif defined(PETSC_USE_REAL___FP16) 108 PETSC_EXTERN MPI_Datatype MPIU___FP16 PetscAttrMPITypeTag(__fp16); 109 #define MPIU_REAL MPIU___FP16 110 typedef __fp16 PetscReal; 111 #define PetscSqrtReal(a) sqrtf(a) 112 #define PetscExpReal(a) expf(a) 113 #define PetscLogReal(a) logf(a) 114 #define PetscLog10Real(a) log10f(a) 115 #ifdef PETSC_HAVE_LOG2 116 #define PetscLog2Real(a) log2f(a) 117 #endif 118 #define PetscSinReal(a) sinf(a) 119 #define PetscCosReal(a) cosf(a) 120 #define PetscTanReal(a) tanf(a) 121 #define PetscAsinReal(a) asinf(a) 122 #define PetscAcosReal(a) acosf(a) 123 #define PetscAtanReal(a) atanf(a) 124 #define PetscAtan2Real(a,b) atan2f(a,b) 125 #define PetscSinhReal(a) sinhf(a) 126 #define PetscCoshReal(a) coshf(a) 127 #define PetscTanhReal(a) tanhf(a) 128 #define PetscPowReal(a,b) powf(a,b) 129 #define PetscCeilReal(a) ceilf(a) 130 #define PetscFloorReal(a) floorf(a) 131 #define PetscFmodReal(a,b) fmodf(a,b) 132 #define PetscTGamma(a) tgammaf(a) 133 #endif /* PETSC_USE_REAL_* */ 134 135 /* 136 Complex number definitions 137 */ 138 #if defined(__cplusplus) && defined(PETSC_HAVE_CXX_COMPLEX) && !defined(PETSC_USE_REAL___FLOAT128) 139 #if !defined(PETSC_SKIP_COMPLEX) 140 #define PETSC_HAVE_COMPLEX 1 141 /* C++ support of complex number */ 142 #if defined(PETSC_HAVE_CUSP) 143 #define complexlib cusp 144 #include <cusp/complex.h> 145 #elif defined(PETSC_HAVE_VECCUDA) && __CUDACC_VER_MAJOR__ > 6 146 /* complex headers in thrust only available in CUDA 7.0 and above */ 147 #define complexlib thrust 148 #include <thrust/complex.h> 149 #else 150 #define complexlib std 151 #include <complex> 152 #endif 153 154 #define PetscRealPartComplex(a) (a).real() 155 #define PetscImaginaryPartComplex(a) (a).imag() 156 #define PetscAbsComplex(a) complexlib::abs(a) 157 #define PetscConjComplex(a) complexlib::conj(a) 158 #define PetscSqrtComplex(a) complexlib::sqrt(a) 159 #define PetscPowComplex(a,b) complexlib::pow(a,b) 160 #define PetscExpComplex(a) complexlib::exp(a) 161 #define PetscLogComplex(a) complexlib::log(a) 162 #define PetscSinComplex(a) complexlib::sin(a) 163 #define PetscCosComplex(a) complexlib::cos(a) 164 #define PetscAsinComplex(a) complexlib::asin(a) 165 #define PetscAcosComplex(a) complexlib::acos(a) 166 #if defined(PETSC_HAVE_TANCOMPLEX) 167 #define PetscTanComplex(a) complexlib::tan(a) 168 #else 169 #define PetscTanComplex(a) PetscSinComplex(a)/PetscCosComplex(a) 170 #endif 171 #define PetscSinhComplex(a) complexlib::sinh(a) 172 #define PetscCoshComplex(a) complexlib::cosh(a) 173 #if defined(PETSC_HAVE_TANHCOMPLEX) 174 #define PetscTanhComplex(a) complexlib::tanh(a) 175 #else 176 #define PetscTanhComplex(a) PetscSinhComplex(a)/PetscCoshComplex(a) 177 #endif 178 179 #if defined(PETSC_USE_REAL_SINGLE) 180 typedef complexlib::complex<float> PetscComplex; 181 #if defined(PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND) 182 static inline PetscComplex operator+(const PetscComplex& lhs, const double& rhs) { return lhs + float(rhs); } 183 static inline PetscComplex operator+(const double& lhs, const PetscComplex& rhs) { return float(lhs) + rhs; } 184 static inline PetscComplex operator-(const PetscComplex& lhs, const double& rhs) { return lhs - float(rhs); } 185 static inline PetscComplex operator-(const double& lhs, const PetscComplex& rhs) { return float(lhs) - rhs; } 186 static inline PetscComplex operator*(const PetscComplex& lhs, const double& rhs) { return lhs * float(rhs); } 187 static inline PetscComplex operator*(const double& lhs, const PetscComplex& rhs) { return float(lhs) * rhs; } 188 static inline PetscComplex operator/(const PetscComplex& lhs, const double& rhs) { return lhs / float(rhs); } 189 static inline PetscComplex operator/(const double& lhs, const PetscComplex& rhs) { return float(lhs) / rhs; } 190 static inline bool operator==(const PetscComplex& lhs, const double& rhs) { return lhs.imag() == float(0) && lhs.real() == float(rhs); } 191 static inline bool operator==(const double& lhs, const PetscComplex& rhs) { return rhs.imag() == float(0) && rhs.real() == float(lhs); } 192 static inline bool operator!=(const PetscComplex& lhs, const double& rhs) { return lhs.imag() != float(0) || lhs.real() != float(rhs); } 193 static inline bool operator!=(const double& lhs, const PetscComplex& rhs) { return rhs.imag() != float(0) || rhs.real() != float(lhs); } 194 #endif /* PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND */ 195 #elif defined(PETSC_USE_REAL_DOUBLE) 196 typedef complexlib::complex<double> PetscComplex; 197 #if defined(PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND) 198 static inline PetscComplex operator+(const PetscComplex& lhs, const PetscInt& rhs) { return lhs + double(rhs); } 199 static inline PetscComplex operator+(const PetscInt& lhs, const PetscComplex& rhs) { return double(lhs) + rhs; } 200 static inline PetscComplex operator-(const PetscComplex& lhs, const PetscInt& rhs) { return lhs - double(rhs); } 201 static inline PetscComplex operator-(const PetscInt& lhs, const PetscComplex& rhs) { return double(lhs) - rhs; } 202 static inline PetscComplex operator*(const PetscComplex& lhs, const PetscInt& rhs) { return lhs * double(rhs); } 203 static inline PetscComplex operator*(const PetscInt& lhs, const PetscComplex& rhs) { return double(lhs) * rhs; } 204 static inline PetscComplex operator/(const PetscComplex& lhs, const PetscInt& rhs) { return lhs / double(rhs); } 205 static inline PetscComplex operator/(const PetscInt& lhs, const PetscComplex& rhs) { return double(lhs) / rhs; } 206 static inline bool operator==(const PetscComplex& lhs, const PetscInt& rhs) { return lhs.imag() == double(0) && lhs.real() == double(rhs); } 207 static inline bool operator==(const PetscInt& lhs, const PetscComplex& rhs) { return rhs.imag() == double(0) && rhs.real() == double(lhs); } 208 static inline bool operator!=(const PetscComplex& lhs, const PetscInt& rhs) { return lhs.imag() != double(0) || lhs.real() != double(rhs); } 209 static inline bool operator!=(const PetscInt& lhs, const PetscComplex& rhs) { return rhs.imag() != double(0) || rhs.real() != double(lhs); } 210 #endif /* PETSC_USE_CXX_COMPLEX_FLOAT_WORKAROUND */ 211 #elif defined(PETSC_USE_REAL___FLOAT128) 212 typedef complexlib::complex<__float128> PetscComplex; /* Notstandard and not expected to work, use __complex128 */ 213 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128; 214 #endif /* PETSC_USE_REAL_ */ 215 #endif /* ! PETSC_SKIP_COMPLEX */ 216 217 #elif !defined(__cplusplus) && defined(PETSC_HAVE_C99_COMPLEX) && !defined(PETSC_USE_REAL___FP16) 218 #if !defined(PETSC_SKIP_COMPLEX) 219 #define PETSC_HAVE_COMPLEX 1 220 #include <complex.h> 221 222 #if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16) 223 typedef float _Complex PetscComplex; 224 225 #define PetscRealPartComplex(a) crealf(a) 226 #define PetscImaginaryPartComplex(a) cimagf(a) 227 #define PetscAbsComplex(a) cabsf(a) 228 #define PetscConjComplex(a) conjf(a) 229 #define PetscSqrtComplex(a) csqrtf(a) 230 #define PetscPowComplex(a,b) cpowf(a,b) 231 #define PetscExpComplex(a) cexpf(a) 232 #define PetscLogComplex(a) clogf(a) 233 #define PetscSinComplex(a) csinf(a) 234 #define PetscCosComplex(a) ccosf(a) 235 #define PetscAsinComplex(a) casinf(a) 236 #define PetscAcosComplex(a) cacosf(a) 237 #define PetscTanComplex(a) ctanf(a) 238 #define PetscSinhComplex(a) csinhf(a) 239 #define PetscCoshComplex(a) ccoshf(a) 240 #define PetscTanhComplex(a) ctanhf(a) 241 242 #elif defined(PETSC_USE_REAL_DOUBLE) 243 typedef double _Complex PetscComplex; 244 245 #define PetscRealPartComplex(a) creal(a) 246 #define PetscImaginaryPartComplex(a) cimag(a) 247 #define PetscAbsComplex(a) cabs(a) 248 #define PetscConjComplex(a) conj(a) 249 #define PetscSqrtComplex(a) csqrt(a) 250 #define PetscPowComplex(a,b) cpow(a,b) 251 #define PetscExpComplex(a) cexp(a) 252 #define PetscLogComplex(a) clog(a) 253 #define PetscSinComplex(a) csin(a) 254 #define PetscCosComplex(a) ccos(a) 255 #define PetscAsinComplex(a) casin(a) 256 #define PetscAcosComplex(a) cacos(a) 257 #define PetscTanComplex(a) ctan(a) 258 #define PetscSinhComplex(a) csinh(a) 259 #define PetscCoshComplex(a) ccosh(a) 260 #define PetscTanhComplex(a) ctanh(a) 261 262 #elif defined(PETSC_USE_REAL___FLOAT128) 263 typedef __complex128 PetscComplex; 264 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 PetscAttrMPITypeTag(__complex128); 265 266 #define PetscRealPartComplex(a) crealq(a) 267 #define PetscImaginaryPartComplex(a) cimagq(a) 268 #define PetscAbsComplex(a) cabsq(a) 269 #define PetscConjComplex(a) conjq(a) 270 #define PetscSqrtComplex(a) csqrtq(a) 271 #define PetscPowComplex(a,b) cpowq(a,b) 272 #define PetscExpComplex(a) cexpq(a) 273 #define PetscLogComplex(a) clogq(a) 274 #define PetscSinComplex(a) csinq(a) 275 #define PetscCosComplex(a) ccosq(a) 276 #define PetscAsinComplex(a) casinq(a) 277 #define PetscAcosComplex(a) cacosq(a) 278 #define PetscTanComplex(a) ctanq(a) 279 #define PetscSinhComplex(a) csinhq(a) 280 #define PetscCoshComplex(a) ccoshq(a) 281 #define PetscTanhComplex(a) ctanhq(a) 282 283 #endif /* PETSC_USE_REAL_* */ 284 #elif (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) 285 #error "PETSc was configured --with-scalar-type=complex, but a language-appropriate complex library is not available" 286 #endif /* !PETSC_SKIP_COMPLEX */ 287 #endif /* (__cplusplus && PETSC_HAVE_CXX_COMPLEX) else-if (!__cplusplus && PETSC_HAVE_C99_COMPLEX) */ 288 289 #if defined(PETSC_HAVE_COMPLEX) 290 #if defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX) 291 #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX 292 #define MPIU_C_COMPLEX MPI_C_COMPLEX 293 #else 294 # if defined(__cplusplus) && defined(PETSC_HAVE_CXX_COMPLEX) 295 typedef complexlib::complex<double> petsc_mpiu_c_double_complex; 296 typedef complexlib::complex<float> petsc_mpiu_c_complex; 297 # elif !defined(__cplusplus) && defined(PETSC_HAVE_C99_COMPLEX) 298 typedef double _Complex petsc_mpiu_c_double_complex; 299 typedef float _Complex petsc_mpiu_c_complex; 300 # else 301 typedef struct {double real,imag;} petsc_mpiu_c_double_complex; 302 typedef struct {float real,imag;} petsc_mpiu_c_complex; 303 # endif 304 PETSC_EXTERN MPI_Datatype MPIU_C_DOUBLE_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_double_complex); 305 PETSC_EXTERN MPI_Datatype MPIU_C_COMPLEX PetscAttrMPITypeTagLayoutCompatible(petsc_mpiu_c_complex); 306 #endif /* PETSC_HAVE_MPI_C_DOUBLE_COMPLEX */ 307 #endif /* PETSC_HAVE_COMPLEX */ 308 309 #if defined(PETSC_HAVE_COMPLEX) 310 # if defined(PETSC_USE_REAL_SINGLE) || defined(PETSC_USE_REAL___FP16) 311 # define MPIU_COMPLEX MPIU_C_COMPLEX 312 # elif defined(PETSC_USE_REAL_DOUBLE) 313 # define MPIU_COMPLEX MPIU_C_DOUBLE_COMPLEX 314 # elif defined(PETSC_USE_REAL___FLOAT128) 315 # define MPIU_COMPLEX MPIU___COMPLEX128 316 # endif /* PETSC_USE_REAL_* */ 317 #endif 318 319 #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX)) 320 typedef PetscComplex PetscScalar; 321 #define PetscRealPart(a) PetscRealPartComplex(a) 322 #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a) 323 #define PetscAbsScalar(a) PetscAbsComplex(a) 324 #define PetscConj(a) PetscConjComplex(a) 325 #define PetscSqrtScalar(a) PetscSqrtComplex(a) 326 #define PetscPowScalar(a,b) PetscPowComplex(a,b) 327 #define PetscExpScalar(a) PetscExpComplex(a) 328 #define PetscLogScalar(a) PetscLogComplex(a) 329 #define PetscSinScalar(a) PetscSinComplex(a) 330 #define PetscCosScalar(a) PetscCosComplex(a) 331 #define PetscAsinScalar(a) PetscAsinComplex(a) 332 #define PetscAcosScalar(a) PetscAcosComplex(a) 333 #define PetscTanScalar(a) PetscTanComplex(a) 334 #define PetscSinhScalar(a) PetscSinhComplex(a) 335 #define PetscCoshScalar(a) PetscCoshComplex(a) 336 #define PetscTanhScalar(a) PetscTanhComplex(a) 337 #define MPIU_SCALAR MPIU_COMPLEX 338 339 /* 340 real number definitions 341 */ 342 #else /* PETSC_USE_COMPLEX */ 343 typedef PetscReal PetscScalar; 344 #define MPIU_SCALAR MPIU_REAL 345 346 #define PetscRealPart(a) (a) 347 #define PetscImaginaryPart(a) ((PetscReal)0.) 348 PETSC_STATIC_INLINE PetscReal PetscAbsScalar(PetscScalar a) {return a < 0.0 ? -a : a;} 349 #define PetscConj(a) (a) 350 #if !defined(PETSC_USE_REAL___FLOAT128) && !defined(PETSC_USE_REAL___FP16) 351 #define PetscSqrtScalar(a) sqrt(a) 352 #define PetscPowScalar(a,b) pow(a,b) 353 #define PetscExpScalar(a) exp(a) 354 #define PetscLogScalar(a) log(a) 355 #define PetscSinScalar(a) sin(a) 356 #define PetscCosScalar(a) cos(a) 357 #define PetscAsinScalar(a) asin(a) 358 #define PetscAcosScalar(a) acos(a) 359 #define PetscTanScalar(a) tan(a) 360 #define PetscSinhScalar(a) sinh(a) 361 #define PetscCoshScalar(a) cosh(a) 362 #define PetscTanhScalar(a) tanh(a) 363 #elif defined(PETSC_USE_REAL___FP16) 364 #define PetscSqrtScalar(a) sqrtf(a) 365 #define PetscPowScalar(a,b) powf(a,b) 366 #define PetscExpScalar(a) expf(a) 367 #define PetscLogScalar(a) logf(a) 368 #define PetscSinScalar(a) sinf(a) 369 #define PetscCosScalar(a) cosf(a) 370 #define PetscAsinScalar(a) asinf(a) 371 #define PetscAcosScalar(a) acosf(a) 372 #define PetscTanScalar(a) tanf(a) 373 #define PetscSinhScalar(a) sinhf(a) 374 #define PetscCoshScalar(a) coshf(a) 375 #define PetscTanhScalar(a) tanhf(a) 376 #else /* PETSC_USE_REAL___FLOAT128 */ 377 #define PetscSqrtScalar(a) sqrtq(a) 378 #define PetscPowScalar(a,b) powq(a,b) 379 #define PetscExpScalar(a) expq(a) 380 #define PetscLogScalar(a) logq(a) 381 #define PetscSinScalar(a) sinq(a) 382 #define PetscCosScalar(a) cosq(a) 383 #define PetscAsinScalar(a) asinq(a) 384 #define PetscAcosScalar(a) acosq(a) 385 #define PetscTanScalar(a) tanq(a) 386 #define PetscSinhScalar(a) sinhq(a) 387 #define PetscCoshScalar(a) coshq(a) 388 #define PetscTanhScalar(a) tanhq(a) 389 #endif /* PETSC_USE_REAL___FLOAT128 */ 390 391 #endif /* PETSC_USE_COMPLEX */ 392 393 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1) 394 #define PetscSignReal(a) (((a) >= 0.0) ? ((a) == 0.0 ? 0.0 : 1.0) : -1.0) 395 #define PetscAbs(a) (((a) >= 0) ? (a) : (-(a))) 396 397 /* --------------------------------------------------------------------------*/ 398 399 /* 400 Certain objects may be created using either single or double precision. 401 This is currently not used. 402 */ 403 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE, PETSC_SCALAR_HALF } PetscScalarPrecision; 404 405 #if defined(PETSC_HAVE_COMPLEX) 406 /* PETSC_i is the imaginary number, i */ 407 PETSC_EXTERN PetscComplex PETSC_i; 408 409 /* Try to do the right thing for complex number construction: see 410 411 http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1464.htm 412 413 for details 414 */ 415 PETSC_STATIC_INLINE PetscComplex PetscCMPLX(PetscReal x, PetscReal y) 416 { 417 #if defined(__cplusplus) 418 return PetscComplex(x,y); 419 #elif defined(_Imaginary_I) 420 return x + y * _Imaginary_I; 421 #else 422 #if defined(PETSC_USE_REAL_SINGLE) && defined(CMPLXF) 423 return CMPLXF(x,y); 424 #elif defined(PETSC_USE_REAL_DOUBLE) && defined(CMPLX) 425 return CMPLX(x,y); 426 #else 427 { /* In both C99 and C11 (ISO/IEC 9899, Section 6.2.5), 428 429 "For each floating type there is a corresponding real type, which is always a real floating 430 type. For real floating types, it is the same type. For complex types, it is the type given 431 by deleting the keyword _Complex from the type name." 432 433 So type punning should be portable. */ 434 union { PetscComplex z; PetscReal f[2]; } uz; 435 436 uz.f[0] = x; 437 uz.f[1] = y; 438 return uz.z; 439 } 440 #endif 441 #endif 442 } 443 #endif 444 445 446 /*MC 447 PetscMin - Returns minimum of two numbers 448 449 Synopsis: 450 #include <petscmath.h> 451 type PetscMin(type v1,type v2) 452 453 Not Collective 454 455 Input Parameter: 456 + v1 - first value to find minimum of 457 - v2 - second value to find minimum of 458 459 Notes: type can be integer or floating point value 460 461 Level: beginner 462 463 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 464 465 M*/ 466 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 467 468 /*MC 469 PetscMax - Returns maxium of two numbers 470 471 Synopsis: 472 #include <petscmath.h> 473 type max PetscMax(type v1,type v2) 474 475 Not Collective 476 477 Input Parameter: 478 + v1 - first value to find maximum of 479 - v2 - second value to find maximum of 480 481 Notes: type can be integer or floating point value 482 483 Level: beginner 484 485 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 486 487 M*/ 488 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 489 490 /*MC 491 PetscClipInterval - Returns a number clipped to be within an interval 492 493 Synopsis: 494 #include <petscmath.h> 495 type clip PetscClipInterval(type x,type a,type b) 496 497 Not Collective 498 499 Input Parameter: 500 + x - value to use if within interval (a,b) 501 . a - lower end of interval 502 - b - upper end of interval 503 504 Notes: type can be integer or floating point value 505 506 Level: beginner 507 508 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 509 510 M*/ 511 #define PetscClipInterval(x,a,b) (PetscMax((a),PetscMin((x),(b)))) 512 513 /*MC 514 PetscAbsInt - Returns the absolute value of an integer 515 516 Synopsis: 517 #include <petscmath.h> 518 int abs PetscAbsInt(int v1) 519 520 Not Collective 521 522 Input Parameter: 523 . v1 - the integer 524 525 Level: beginner 526 527 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 528 529 M*/ 530 #define PetscAbsInt(a) (((a)<0) ? (-(a)) : (a)) 531 532 /*MC 533 PetscAbsReal - Returns the absolute value of an real number 534 535 Synopsis: 536 #include <petscmath.h> 537 Real abs PetscAbsReal(PetscReal v1) 538 539 Not Collective 540 541 Input Parameter: 542 . v1 - the double 543 544 545 Level: beginner 546 547 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 548 549 M*/ 550 #if defined(PETSC_USE_REAL_SINGLE) 551 #define PetscAbsReal(a) fabsf(a) 552 #elif defined(PETSC_USE_REAL_DOUBLE) 553 #define PetscAbsReal(a) fabs(a) 554 #elif defined(PETSC_USE_REAL___FLOAT128) 555 #define PetscAbsReal(a) fabsq(a) 556 #elif defined(PETSC_USE_REAL___FP16) 557 #define PetscAbsReal(a) fabsf(a) 558 #endif 559 560 /*MC 561 PetscSqr - Returns the square of a number 562 563 Synopsis: 564 #include <petscmath.h> 565 type sqr PetscSqr(type v1) 566 567 Not Collective 568 569 Input Parameter: 570 . v1 - the value 571 572 Notes: type can be integer or floating point value 573 574 Level: beginner 575 576 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 577 578 M*/ 579 #define PetscSqr(a) ((a)*(a)) 580 581 /* ----------------------------------------------------------------------------*/ 582 583 #if defined(PETSC_USE_REAL_SINGLE) 584 #define PetscRealConstant(constant) constant##F 585 #elif defined(PETSC_USE_REAL___FLOAT128) 586 #define PetscRealConstant(constant) constant##Q 587 #else 588 #define PetscRealConstant(constant) constant 589 #endif 590 591 /* 592 Basic constants 593 */ 594 #if defined(PETSC_USE_REAL___FLOAT128) 595 #define PETSC_PI M_PIq 596 #elif defined(M_PI) 597 #define PETSC_PI M_PI 598 #else 599 #define PETSC_PI 3.14159265358979323846264338327950288419716939937510582 600 #endif 601 #define PETSC_PHI 1.6180339887498948482 602 603 #if !defined(PETSC_USE_64BIT_INDICES) 604 #define PETSC_MAX_INT 2147483647 605 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 606 #else 607 #define PETSC_MAX_INT 9223372036854775807L 608 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 609 #endif 610 611 #if defined(PETSC_USE_REAL_SINGLE) 612 # define PETSC_MAX_REAL 3.40282346638528860e+38F 613 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 614 # define PETSC_MACHINE_EPSILON 1.19209290e-07F 615 # define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F 616 # define PETSC_SMALL 1.e-5F 617 #elif defined(PETSC_USE_REAL_DOUBLE) 618 # define PETSC_MAX_REAL 1.7976931348623157e+308 619 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 620 # define PETSC_MACHINE_EPSILON 2.2204460492503131e-16 621 # define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08 622 # define PETSC_SMALL 1.e-10 623 #elif defined(PETSC_USE_REAL___FLOAT128) 624 # define PETSC_MAX_REAL FLT128_MAX 625 # define PETSC_MIN_REAL (-FLT128_MAX) 626 # define PETSC_MACHINE_EPSILON FLT128_EPSILON 627 # define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q 628 # define PETSC_SMALL 1.e-20Q 629 #elif defined(PETSC_USE_REAL___FP16) /* maybe should use single precision values for these? */ 630 # define PETSC_MAX_REAL 65504. 631 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 632 # define PETSC_MACHINE_EPSILON .00097656 633 # define PETSC_SQRT_MACHINE_EPSILON .0312 634 # define PETSC_SMALL 5.e-3 635 #endif 636 637 #define PETSC_INFINITY (PETSC_MAX_REAL/4) 638 #define PETSC_NINFINITY (-PETSC_INFINITY) 639 640 PETSC_EXTERN PetscBool PetscIsInfReal(PetscReal); 641 PETSC_EXTERN PetscBool PetscIsNanReal(PetscReal); 642 PETSC_EXTERN PetscBool PetscIsNormalReal(PetscReal); 643 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanReal(PetscReal v) {return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;} 644 PETSC_STATIC_INLINE PetscBool PetscIsInfScalar(PetscScalar v) {return PetscIsInfReal(PetscAbsScalar(v));} 645 PETSC_STATIC_INLINE PetscBool PetscIsNanScalar(PetscScalar v) {return PetscIsNanReal(PetscAbsScalar(v));} 646 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanScalar(PetscScalar v) {return PetscIsInfOrNanReal(PetscAbsScalar(v));} 647 PETSC_STATIC_INLINE PetscBool PetscIsNormalScalar(PetscScalar v) {return PetscIsNormalReal(PetscAbsScalar(v));} 648 649 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal,PetscReal); 650 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar,PetscScalar); 651 652 /* 653 These macros are currently hardwired to match the regular data types, so there is no support for a different 654 MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again. 655 */ 656 #define MPIU_MATSCALAR MPIU_SCALAR 657 typedef PetscScalar MatScalar; 658 typedef PetscReal MatReal; 659 660 struct petsc_mpiu_2scalar {PetscScalar a,b;}; 661 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar); 662 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT) 663 struct petsc_mpiu_2int {PetscInt a,b;}; 664 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int); 665 #else 666 #define MPIU_2INT MPI_2INT 667 #endif 668 669 PETSC_STATIC_INLINE PetscInt PetscPowInt(PetscInt base,PetscInt power) 670 { 671 PetscInt result = 1; 672 while (power) { 673 if (power & 1) result *= base; 674 power >>= 1; 675 base *= base; 676 } 677 return result; 678 } 679 680 PETSC_STATIC_INLINE PetscReal PetscPowRealInt(PetscReal base,PetscInt power) 681 { 682 PetscReal result = 1; 683 if (power < 0) { 684 power = -power; 685 base = ((PetscReal)1)/base; 686 } 687 while (power) { 688 if (power & 1) result *= base; 689 power >>= 1; 690 base *= base; 691 } 692 return result; 693 } 694 695 PETSC_STATIC_INLINE PetscScalar PetscPowScalarInt(PetscScalar base,PetscInt power) 696 { 697 PetscScalar result = 1; 698 if (power < 0) { 699 power = -power; 700 base = ((PetscReal)1)/base; 701 } 702 while (power) { 703 if (power & 1) result *= base; 704 power >>= 1; 705 base *= base; 706 } 707 return result; 708 } 709 710 PETSC_STATIC_INLINE PetscScalar PetscPowScalarReal(PetscScalar base,PetscReal power) 711 { 712 PetscScalar cpower = power; 713 return PetscPowScalar(base,cpower); 714 } 715 716 #ifndef PETSC_HAVE_LOG2 717 PETSC_STATIC_INLINE PetscReal PetscLog2Real(PetscReal n) 718 { 719 return PetscLogReal(n)/PetscLogReal(2.0); 720 } 721 #endif 722 #endif 723