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 { /* In both C99 and C11 (ISO/IEC 9899, Section 6.2.5), 423 424 "For each floating type there is a corresponding real type, which is always a real floating 425 type. For real floating types, it is the same type. For complex types, it is the type given 426 by deleting the keyword _Complex from the type name." 427 428 So type punning should be portable. */ 429 union { PetscComplex z; PetscReal f[2]; } uz; 430 431 uz.f[0] = x; 432 uz.f[1] = y; 433 return uz.z; 434 } 435 #endif 436 } 437 #endif 438 439 440 /*MC 441 PetscMin - Returns minimum of two numbers 442 443 Synopsis: 444 #include <petscmath.h> 445 type PetscMin(type v1,type v2) 446 447 Not Collective 448 449 Input Parameter: 450 + v1 - first value to find minimum of 451 - v2 - second value to find minimum of 452 453 Notes: type can be integer or floating point value 454 455 Level: beginner 456 457 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 458 459 M*/ 460 #define PetscMin(a,b) (((a)<(b)) ? (a) : (b)) 461 462 /*MC 463 PetscMax - Returns maxium of two numbers 464 465 Synopsis: 466 #include <petscmath.h> 467 type max PetscMax(type v1,type v2) 468 469 Not Collective 470 471 Input Parameter: 472 + v1 - first value to find maximum of 473 - v2 - second value to find maximum of 474 475 Notes: type can be integer or floating point value 476 477 Level: beginner 478 479 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 480 481 M*/ 482 #define PetscMax(a,b) (((a)<(b)) ? (b) : (a)) 483 484 /*MC 485 PetscClipInterval - Returns a number clipped to be within an interval 486 487 Synopsis: 488 #include <petscmath.h> 489 type clip PetscClipInterval(type x,type a,type b) 490 491 Not Collective 492 493 Input Parameter: 494 + x - value to use if within interval (a,b) 495 . a - lower end of interval 496 - b - upper end of interval 497 498 Notes: type can be integer or floating point value 499 500 Level: beginner 501 502 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr() 503 504 M*/ 505 #define PetscClipInterval(x,a,b) (PetscMax((a),PetscMin((x),(b)))) 506 507 /*MC 508 PetscAbsInt - Returns the absolute value of an integer 509 510 Synopsis: 511 #include <petscmath.h> 512 int abs PetscAbsInt(int v1) 513 514 Not Collective 515 516 Input Parameter: 517 . v1 - the integer 518 519 Level: beginner 520 521 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr() 522 523 M*/ 524 #define PetscAbsInt(a) (((a)<0) ? (-(a)) : (a)) 525 526 /*MC 527 PetscAbsReal - Returns the absolute value of an real number 528 529 Synopsis: 530 #include <petscmath.h> 531 Real abs PetscAbsReal(PetscReal v1) 532 533 Not Collective 534 535 Input Parameter: 536 . v1 - the double 537 538 539 Level: beginner 540 541 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr() 542 543 M*/ 544 #if defined(PETSC_USE_REAL_SINGLE) 545 #define PetscAbsReal(a) fabsf(a) 546 #elif defined(PETSC_USE_REAL_DOUBLE) 547 #define PetscAbsReal(a) fabs(a) 548 #elif defined(PETSC_USE_REAL___FLOAT128) 549 #define PetscAbsReal(a) fabsq(a) 550 #elif defined(PETSC_USE_REAL___FP16) 551 #define PetscAbsReal(a) fabsf(a) 552 #endif 553 554 /*MC 555 PetscSqr - Returns the square of a number 556 557 Synopsis: 558 #include <petscmath.h> 559 type sqr PetscSqr(type v1) 560 561 Not Collective 562 563 Input Parameter: 564 . v1 - the value 565 566 Notes: type can be integer or floating point value 567 568 Level: beginner 569 570 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal() 571 572 M*/ 573 #define PetscSqr(a) ((a)*(a)) 574 575 /* ----------------------------------------------------------------------------*/ 576 577 #if defined(PETSC_USE_REAL_SINGLE) 578 #define PetscRealConstant(constant) constant##F 579 #elif defined(PETSC_USE_REAL___FLOAT128) 580 #define PetscRealConstant(constant) constant##Q 581 #else 582 #define PetscRealConstant(constant) constant 583 #endif 584 585 /* 586 Basic constants 587 */ 588 #if defined(PETSC_USE_REAL___FLOAT128) 589 #define PETSC_PI M_PIq 590 #elif defined(M_PI) 591 #define PETSC_PI M_PI 592 #else 593 #define PETSC_PI 3.14159265358979323846264338327950288419716939937510582 594 #endif 595 #define PETSC_PHI 1.6180339887498948482 596 597 #if !defined(PETSC_USE_64BIT_INDICES) 598 #define PETSC_MAX_INT 2147483647 599 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 600 #else 601 #define PETSC_MAX_INT 9223372036854775807L 602 #define PETSC_MIN_INT (-PETSC_MAX_INT - 1) 603 #endif 604 605 #if defined(PETSC_USE_REAL_SINGLE) 606 # define PETSC_MAX_REAL 3.40282346638528860e+38F 607 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 608 # define PETSC_MACHINE_EPSILON 1.19209290e-07F 609 # define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F 610 # define PETSC_SMALL 1.e-5F 611 #elif defined(PETSC_USE_REAL_DOUBLE) 612 # define PETSC_MAX_REAL 1.7976931348623157e+308 613 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 614 # define PETSC_MACHINE_EPSILON 2.2204460492503131e-16 615 # define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08 616 # define PETSC_SMALL 1.e-10 617 #elif defined(PETSC_USE_REAL___FLOAT128) 618 # define PETSC_MAX_REAL FLT128_MAX 619 # define PETSC_MIN_REAL (-FLT128_MAX) 620 # define PETSC_MACHINE_EPSILON FLT128_EPSILON 621 # define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q 622 # define PETSC_SMALL 1.e-20Q 623 #elif defined(PETSC_USE_REAL___FP16) /* maybe should use single precision values for these? */ 624 # define PETSC_MAX_REAL 65504. 625 # define PETSC_MIN_REAL (-PETSC_MAX_REAL) 626 # define PETSC_MACHINE_EPSILON .00097656 627 # define PETSC_SQRT_MACHINE_EPSILON .0312 628 # define PETSC_SMALL 5.e-3 629 #endif 630 631 #define PETSC_INFINITY (PETSC_MAX_REAL/4) 632 #define PETSC_NINFINITY (-PETSC_INFINITY) 633 634 PETSC_EXTERN PetscBool PetscIsInfReal(PetscReal); 635 PETSC_EXTERN PetscBool PetscIsNanReal(PetscReal); 636 PETSC_EXTERN PetscBool PetscIsNormalReal(PetscReal); 637 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanReal(PetscReal v) {return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;} 638 PETSC_STATIC_INLINE PetscBool PetscIsInfScalar(PetscScalar v) {return PetscIsInfReal(PetscAbsScalar(v));} 639 PETSC_STATIC_INLINE PetscBool PetscIsNanScalar(PetscScalar v) {return PetscIsNanReal(PetscAbsScalar(v));} 640 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanScalar(PetscScalar v) {return PetscIsInfOrNanReal(PetscAbsScalar(v));} 641 PETSC_STATIC_INLINE PetscBool PetscIsNormalScalar(PetscScalar v) {return PetscIsNormalReal(PetscAbsScalar(v));} 642 643 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal,PetscReal); 644 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar,PetscScalar); 645 646 /* 647 These macros are currently hardwired to match the regular data types, so there is no support for a different 648 MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again. 649 */ 650 #define MPIU_MATSCALAR MPIU_SCALAR 651 typedef PetscScalar MatScalar; 652 typedef PetscReal MatReal; 653 654 struct petsc_mpiu_2scalar {PetscScalar a,b;}; 655 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar); 656 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT) 657 struct petsc_mpiu_2int {PetscInt a,b;}; 658 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int); 659 #else 660 #define MPIU_2INT MPI_2INT 661 #endif 662 663 PETSC_STATIC_INLINE PetscInt PetscPowInt(PetscInt base,PetscInt power) 664 { 665 PetscInt result = 1; 666 while (power) { 667 if (power & 1) result *= base; 668 power >>= 1; 669 base *= base; 670 } 671 return result; 672 } 673 674 PETSC_STATIC_INLINE PetscReal PetscPowRealInt(PetscReal base,PetscInt power) 675 { 676 PetscReal result = 1; 677 if (power < 0) { 678 power = -power; 679 base = ((PetscReal)1)/base; 680 } 681 while (power) { 682 if (power & 1) result *= base; 683 power >>= 1; 684 base *= base; 685 } 686 return result; 687 } 688 689 PETSC_STATIC_INLINE PetscScalar PetscPowScalarInt(PetscScalar base,PetscInt power) 690 { 691 PetscScalar result = 1; 692 if (power < 0) { 693 power = -power; 694 base = ((PetscReal)1)/base; 695 } 696 while (power) { 697 if (power & 1) result *= base; 698 power >>= 1; 699 base *= base; 700 } 701 return result; 702 } 703 704 PETSC_STATIC_INLINE PetscScalar PetscPowScalarReal(PetscScalar base,PetscReal power) 705 { 706 PetscScalar cpower = power; 707 return PetscPowScalar(base,cpower); 708 } 709 710 #ifndef PETSC_HAVE_LOG2 711 PETSC_STATIC_INLINE PetscReal PetscLog2Real(PetscReal n) 712 { 713 return PetscLogReal(n)/PetscLogReal(2.0); 714 } 715 #endif 716 #endif 717