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