xref: /petsc/include/petscmath.h (revision 487a658c8b32ba712a1dc8280daad2fd70c1dcd9)
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 PetscRoundReal(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:
495     type can be integer or floating point value
496 
497    Level: beginner
498 
499 .seealso: PetscMax(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
500 
501 M*/
502 #define PetscMin(a,b)   (((a)<(b)) ?  (a) : (b))
503 
504 /*MC
505    PetscMax - Returns maxium of two numbers
506 
507    Synopsis:
508    #include <petscmath.h>
509    type max PetscMax(type v1,type v2)
510 
511    Not Collective
512 
513    Input Parameter:
514 +  v1 - first value to find maximum of
515 -  v2 - second value to find maximum of
516 
517    Notes:
518     type can be integer or floating point value
519 
520    Level: beginner
521 
522 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
523 
524 M*/
525 #define PetscMax(a,b)   (((a)<(b)) ?  (b) : (a))
526 
527 /*MC
528    PetscClipInterval - Returns a number clipped to be within an interval
529 
530    Synopsis:
531    #include <petscmath.h>
532    type clip PetscClipInterval(type x,type a,type b)
533 
534    Not Collective
535 
536    Input Parameter:
537 +  x - value to use if within interval (a,b)
538 .  a - lower end of interval
539 -  b - upper end of interval
540 
541    Notes:
542     type can be integer or floating point value
543 
544    Level: beginner
545 
546 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
547 
548 M*/
549 #define PetscClipInterval(x,a,b)   (PetscMax((a),PetscMin((x),(b))))
550 
551 /*MC
552    PetscAbsInt - Returns the absolute value of an integer
553 
554    Synopsis:
555    #include <petscmath.h>
556    int abs PetscAbsInt(int v1)
557 
558    Not Collective
559 
560    Input Parameter:
561 .   v1 - the integer
562 
563    Level: beginner
564 
565 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr()
566 
567 M*/
568 #define PetscAbsInt(a)  (((a)<0)   ? (-(a)) : (a))
569 
570 /*MC
571    PetscAbsReal - Returns the absolute value of an real number
572 
573    Synopsis:
574    #include <petscmath.h>
575    Real abs PetscAbsReal(PetscReal v1)
576 
577    Not Collective
578 
579    Input Parameter:
580 .   v1 - the double
581 
582 
583    Level: beginner
584 
585 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr()
586 
587 M*/
588 #if defined(PETSC_USE_REAL_SINGLE)
589 #define PetscAbsReal(a) fabsf(a)
590 #elif defined(PETSC_USE_REAL_DOUBLE)
591 #define PetscAbsReal(a) fabs(a)
592 #elif defined(PETSC_USE_REAL___FLOAT128)
593 #define PetscAbsReal(a) fabsq(a)
594 #elif defined(PETSC_USE_REAL___FP16)
595 #define PetscAbsReal(a) fabsf(a)
596 #endif
597 
598 /*MC
599    PetscSqr - Returns the square of a number
600 
601    Synopsis:
602    #include <petscmath.h>
603    type sqr PetscSqr(type v1)
604 
605    Not Collective
606 
607    Input Parameter:
608 .   v1 - the value
609 
610    Notes:
611     type can be integer or floating point value
612 
613    Level: beginner
614 
615 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal()
616 
617 M*/
618 #define PetscSqr(a)     ((a)*(a))
619 
620 /* ----------------------------------------------------------------------------*/
621 
622 #if defined(PETSC_USE_REAL_SINGLE)
623 #define PetscRealConstant(constant) constant##F
624 #elif defined(PETSC_USE_REAL___FLOAT128)
625 #define PetscRealConstant(constant) constant##Q
626 #else
627 #define PetscRealConstant(constant) constant
628 #endif
629 
630 /*
631      Basic constants
632 */
633 #define PETSC_PI   PetscRealConstant(3.1415926535897932384626433832795029)
634 #define PETSC_PHI  PetscRealConstant(1.6180339887498948482045868343656381)
635 
636 #if !defined(PETSC_USE_64BIT_INDICES)
637 #define PETSC_MAX_INT            2147483647
638 #define PETSC_MIN_INT            (-PETSC_MAX_INT - 1)
639 #else
640 #define PETSC_MAX_INT            9223372036854775807L
641 #define PETSC_MIN_INT            (-PETSC_MAX_INT - 1)
642 #endif
643 
644 #if defined(PETSC_USE_REAL_SINGLE)
645 #  define PETSC_MAX_REAL                3.40282346638528860e+38F
646 #  define PETSC_MIN_REAL                (-PETSC_MAX_REAL)
647 #  define PETSC_MACHINE_EPSILON         1.19209290e-07F
648 #  define PETSC_SQRT_MACHINE_EPSILON    3.45266983e-04F
649 #  define PETSC_SMALL                   1.e-5F
650 #elif defined(PETSC_USE_REAL_DOUBLE)
651 #  define PETSC_MAX_REAL                1.7976931348623157e+308
652 #  define PETSC_MIN_REAL                (-PETSC_MAX_REAL)
653 #  define PETSC_MACHINE_EPSILON         2.2204460492503131e-16
654 #  define PETSC_SQRT_MACHINE_EPSILON    1.490116119384766e-08
655 #  define PETSC_SMALL                   1.e-10
656 #elif defined(PETSC_USE_REAL___FLOAT128)
657 #  define PETSC_MAX_REAL                FLT128_MAX
658 #  define PETSC_MIN_REAL                (-FLT128_MAX)
659 #  define PETSC_MACHINE_EPSILON         FLT128_EPSILON
660 #  define PETSC_SQRT_MACHINE_EPSILON    1.38777878078144567552953958511352539e-17Q
661 #  define PETSC_SMALL                   1.e-20Q
662 #elif defined(PETSC_USE_REAL___FP16)  /* maybe should use single precision values for these? */
663 #  define PETSC_MAX_REAL                65504.
664 #  define PETSC_MIN_REAL                (-PETSC_MAX_REAL)
665 #  define PETSC_MACHINE_EPSILON         .00097656
666 #  define PETSC_SQRT_MACHINE_EPSILON    .0312
667 #  define PETSC_SMALL                   5.e-3
668 #endif
669 
670 #define PETSC_INFINITY               (PETSC_MAX_REAL/4)
671 #define PETSC_NINFINITY              (-PETSC_INFINITY)
672 
673 PETSC_EXTERN PetscBool PetscIsInfReal(PetscReal);
674 PETSC_EXTERN PetscBool PetscIsNanReal(PetscReal);
675 PETSC_EXTERN PetscBool PetscIsNormalReal(PetscReal);
676 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanReal(PetscReal v) {return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;}
677 PETSC_STATIC_INLINE PetscBool PetscIsInfScalar(PetscScalar v) {return PetscIsInfReal(PetscAbsScalar(v));}
678 PETSC_STATIC_INLINE PetscBool PetscIsNanScalar(PetscScalar v) {return PetscIsNanReal(PetscAbsScalar(v));}
679 PETSC_STATIC_INLINE PetscBool PetscIsInfOrNanScalar(PetscScalar v) {return PetscIsInfOrNanReal(PetscAbsScalar(v));}
680 PETSC_STATIC_INLINE PetscBool PetscIsNormalScalar(PetscScalar v) {return PetscIsNormalReal(PetscAbsScalar(v));}
681 
682 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal,PetscReal);
683 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar,PetscScalar);
684 
685 /*
686     These macros are currently hardwired to match the regular data types, so there is no support for a different
687     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
688  */
689 #define MPIU_MATSCALAR MPIU_SCALAR
690 typedef PetscScalar MatScalar;
691 typedef PetscReal MatReal;
692 
693 struct petsc_mpiu_2scalar {PetscScalar a,b;};
694 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2scalar);
695 #if defined(PETSC_USE_64BIT_INDICES) || !defined(MPI_2INT)
696 struct petsc_mpiu_2int {PetscInt a,b;};
697 PETSC_EXTERN MPI_Datatype MPIU_2INT PetscAttrMPITypeTagLayoutCompatible(struct petsc_mpiu_2int);
698 #else
699 #define MPIU_2INT MPI_2INT
700 #endif
701 
702 PETSC_STATIC_INLINE PetscInt PetscPowInt(PetscInt base,PetscInt power)
703 {
704   PetscInt result = 1;
705   while (power) {
706     if (power & 1) result *= base;
707     power >>= 1;
708     base *= base;
709   }
710   return result;
711 }
712 
713 PETSC_STATIC_INLINE PetscReal PetscPowRealInt(PetscReal base,PetscInt power)
714 {
715   PetscReal result = 1;
716   if (power < 0) {
717     power = -power;
718     base  = ((PetscReal)1)/base;
719   }
720   while (power) {
721     if (power & 1) result *= base;
722     power >>= 1;
723     base *= base;
724   }
725   return result;
726 }
727 
728 PETSC_STATIC_INLINE PetscScalar PetscPowScalarInt(PetscScalar base,PetscInt power)
729 {
730   PetscScalar result = 1;
731   if (power < 0) {
732     power = -power;
733     base  = ((PetscReal)1)/base;
734   }
735   while (power) {
736     if (power & 1) result *= base;
737     power >>= 1;
738     base *= base;
739   }
740   return result;
741 }
742 
743 PETSC_STATIC_INLINE PetscScalar PetscPowScalarReal(PetscScalar base,PetscReal power)
744 {
745   PetscScalar cpower = power;
746   return PetscPowScalar(base,cpower);
747 }
748 
749 #ifndef PETSC_HAVE_LOG2
750 PETSC_STATIC_INLINE PetscReal PetscLog2Real(PetscReal n)
751 {
752   return PetscLogReal(n)/PetscLogReal(2.0);
753 }
754 #endif
755 #endif
756