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