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