xref: /petsc/include/petscmath.h (revision 9140fee14acbea959c6ed74f4836a1a89f708038)
1 /*
2     PETSc mathematics include file. Defines certain basic mathematical
3     constants and functions for working with single, double, and quad precision
4     floating point numbers as well as complex single and double.
5 
6     This file is included by petscsys.h and should not be used directly.
7 */
8 #pragma once
9 
10 #include <math.h>
11 #include <petscmacros.h>
12 #include <petscsystypes.h>
13 
14 /* SUBMANSEC = Sys */
15 
16 /*
17    Defines operations that are different for complex and real numbers.
18    All PETSc objects in one program are built around the object
19    PetscScalar which is either always a real or a complex.
20 */
21 
22 /*
23     Real number definitions
24  */
25 #if defined(PETSC_USE_REAL_SINGLE)
26   #define PetscSqrtReal(a)        sqrtf(a)
27   #define PetscCbrtReal(a)        cbrtf(a)
28   #define PetscHypotReal(a, b)    hypotf(a, b)
29   #define PetscAtan2Real(a, b)    atan2f(a, b)
30   #define PetscPowReal(a, b)      powf(a, b)
31   #define PetscExpReal(a)         expf(a)
32   #define PetscLogReal(a)         logf(a)
33   #define PetscLog10Real(a)       log10f(a)
34   #define PetscLog2Real(a)        log2f(a)
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 PetscSinhReal(a)        sinhf(a)
42   #define PetscCoshReal(a)        coshf(a)
43   #define PetscTanhReal(a)        tanhf(a)
44   #define PetscAsinhReal(a)       asinhf(a)
45   #define PetscAcoshReal(a)       acoshf(a)
46   #define PetscAtanhReal(a)       atanhf(a)
47   #define PetscErfReal(a)         erff(a)
48   #define PetscCeilReal(a)        ceilf(a)
49   #define PetscFloorReal(a)       floorf(a)
50   #define PetscFmodReal(a, b)     fmodf(a, b)
51   #define PetscCopysignReal(a, b) copysignf(a, b)
52   #define PetscTGamma(a)          tgammaf(a)
53   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
54     #define PetscLGamma(a) gammaf(a)
55   #else
56     #define PetscLGamma(a) lgammaf(a)
57   #endif
58 
59 #elif defined(PETSC_USE_REAL_DOUBLE)
60   #define PetscSqrtReal(a)        sqrt(a)
61   #define PetscCbrtReal(a)        cbrt(a)
62   #define PetscHypotReal(a, b)    hypot(a, b)
63   #define PetscAtan2Real(a, b)    atan2(a, b)
64   #define PetscPowReal(a, b)      pow(a, b)
65   #define PetscExpReal(a)         exp(a)
66   #define PetscLogReal(a)         log(a)
67   #define PetscLog10Real(a)       log10(a)
68   #define PetscLog2Real(a)        log2(a)
69   #define PetscSinReal(a)         sin(a)
70   #define PetscCosReal(a)         cos(a)
71   #define PetscTanReal(a)         tan(a)
72   #define PetscAsinReal(a)        asin(a)
73   #define PetscAcosReal(a)        acos(a)
74   #define PetscAtanReal(a)        atan(a)
75   #define PetscSinhReal(a)        sinh(a)
76   #define PetscCoshReal(a)        cosh(a)
77   #define PetscTanhReal(a)        tanh(a)
78   #define PetscAsinhReal(a)       asinh(a)
79   #define PetscAcoshReal(a)       acosh(a)
80   #define PetscAtanhReal(a)       atanh(a)
81   #define PetscErfReal(a)         erf(a)
82   #define PetscCeilReal(a)        ceil(a)
83   #define PetscFloorReal(a)       floor(a)
84   #define PetscFmodReal(a, b)     fmod(a, b)
85   #define PetscCopysignReal(a, b) copysign(a, b)
86   #define PetscTGamma(a)          tgamma(a)
87   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
88     #define PetscLGamma(a) gamma(a)
89   #else
90     #define PetscLGamma(a) lgamma(a)
91   #endif
92 
93 #elif defined(PETSC_USE_REAL___FLOAT128)
94   #define PetscSqrtReal(a)        sqrtq(a)
95   #define PetscCbrtReal(a)        cbrtq(a)
96   #define PetscHypotReal(a, b)    hypotq(a, b)
97   #define PetscAtan2Real(a, b)    atan2q(a, b)
98   #define PetscPowReal(a, b)      powq(a, b)
99   #define PetscExpReal(a)         expq(a)
100   #define PetscLogReal(a)         logq(a)
101   #define PetscLog10Real(a)       log10q(a)
102   #define PetscLog2Real(a)        log2q(a)
103   #define PetscSinReal(a)         sinq(a)
104   #define PetscCosReal(a)         cosq(a)
105   #define PetscTanReal(a)         tanq(a)
106   #define PetscAsinReal(a)        asinq(a)
107   #define PetscAcosReal(a)        acosq(a)
108   #define PetscAtanReal(a)        atanq(a)
109   #define PetscSinhReal(a)        sinhq(a)
110   #define PetscCoshReal(a)        coshq(a)
111   #define PetscTanhReal(a)        tanhq(a)
112   #define PetscAsinhReal(a)       asinhq(a)
113   #define PetscAcoshReal(a)       acoshq(a)
114   #define PetscAtanhReal(a)       atanhq(a)
115   #define PetscErfReal(a)         erfq(a)
116   #define PetscCeilReal(a)        ceilq(a)
117   #define PetscFloorReal(a)       floorq(a)
118   #define PetscFmodReal(a, b)     fmodq(a, b)
119   #define PetscCopysignReal(a, b) copysignq(a, b)
120   #define PetscTGamma(a)          tgammaq(a)
121   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
122     #define PetscLGamma(a) gammaq(a)
123   #else
124     #define PetscLGamma(a) lgammaq(a)
125   #endif
126 
127 #elif defined(PETSC_USE_REAL___FP16)
128   #define PetscSqrtReal(a)        sqrtf(a)
129   #define PetscCbrtReal(a)        cbrtf(a)
130   #define PetscHypotReal(a, b)    hypotf(a, b)
131   #define PetscAtan2Real(a, b)    atan2f(a, b)
132   #define PetscPowReal(a, b)      powf(a, b)
133   #define PetscExpReal(a)         expf(a)
134   #define PetscLogReal(a)         logf(a)
135   #define PetscLog10Real(a)       log10f(a)
136   #define PetscLog2Real(a)        log2f(a)
137   #define PetscSinReal(a)         sinf(a)
138   #define PetscCosReal(a)         cosf(a)
139   #define PetscTanReal(a)         tanf(a)
140   #define PetscAsinReal(a)        asinf(a)
141   #define PetscAcosReal(a)        acosf(a)
142   #define PetscAtanReal(a)        atanf(a)
143   #define PetscSinhReal(a)        sinhf(a)
144   #define PetscCoshReal(a)        coshf(a)
145   #define PetscTanhReal(a)        tanhf(a)
146   #define PetscAsinhReal(a)       asinhf(a)
147   #define PetscAcoshReal(a)       acoshf(a)
148   #define PetscAtanhReal(a)       atanhf(a)
149   #define PetscErfReal(a)         erff(a)
150   #define PetscCeilReal(a)        ceilf(a)
151   #define PetscFloorReal(a)       floorf(a)
152   #define PetscFmodReal(a, b)     fmodf(a, b)
153   #define PetscCopySignReal(a, b) copysignf(a, b)
154   #define PetscTGamma(a)          tgammaf(a)
155   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
156     #define PetscLGamma(a) gammaf(a)
157   #else
158     #define PetscLGamma(a) lgammaf(a)
159   #endif
160 
161 #endif /* PETSC_USE_REAL_* */
162 
163 static inline PetscReal PetscSignReal(PetscReal a)
164 {
165   return (PetscReal)((a < (PetscReal)0) ? -1 : ((a > (PetscReal)0) ? 1 : 0));
166 }
167 
168 #if !defined(PETSC_HAVE_LOG2)
169   #undef PetscLog2Real
170 static inline PetscReal PetscLog2Real(PetscReal a)
171 {
172   return PetscLogReal(a) / PetscLogReal((PetscReal)2);
173 }
174 #endif
175 
176 #if defined(PETSC_HAVE_REAL___FLOAT128) && !defined(PETSC_SKIP_REAL___FLOAT128)
177 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128 PETSC_ATTRIBUTE_MPI_TYPE_TAG(__float128);
178 #endif
179 #if defined(PETSC_HAVE_REAL___FP16) && !defined(PETSC_SKIP_REAL___FP16)
180 PETSC_EXTERN MPI_Datatype MPIU___FP16 PETSC_ATTRIBUTE_MPI_TYPE_TAG(__fp16);
181 #endif
182 
183 /*MC
184    MPIU_REAL - Portable MPI datatype corresponding to `PetscReal` independent of what precision `PetscReal` is in
185 
186    Notes:
187    In MPI calls that require an MPI datatype that matches a `PetscReal` or array of `PetscReal` values, pass this value.
188 
189    Level: beginner
190 
191 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_SCALAR`, `MPIU_COMPLEX`, `MPIU_INT`
192 M*/
193 #if defined(PETSC_USE_REAL_SINGLE)
194   #define MPIU_REAL MPI_FLOAT
195 #elif defined(PETSC_USE_REAL_DOUBLE)
196   #define MPIU_REAL MPI_DOUBLE
197 #elif defined(PETSC_USE_REAL___FLOAT128)
198   #define MPIU_REAL MPIU___FLOAT128
199 #elif defined(PETSC_USE_REAL___FP16)
200   #define MPIU_REAL MPIU___FP16
201 #endif /* PETSC_USE_REAL_* */
202 
203 /*
204     Complex number definitions
205  */
206 #if defined(PETSC_HAVE_COMPLEX)
207   #if defined(__cplusplus) && !defined(PETSC_USE_REAL___FLOAT128)
208   /* C++ support of complex number */
209 
210     #define PetscRealPartComplex(a)      (static_cast<PetscComplex>(a)).real()
211     #define PetscImaginaryPartComplex(a) (static_cast<PetscComplex>(a)).imag()
212     #define PetscAbsComplex(a)           petsccomplexlib::abs(static_cast<PetscComplex>(a))
213     #define PetscArgComplex(a)           petsccomplexlib::arg(static_cast<PetscComplex>(a))
214     #define PetscConjComplex(a)          petsccomplexlib::conj(static_cast<PetscComplex>(a))
215     #define PetscSqrtComplex(a)          petsccomplexlib::sqrt(static_cast<PetscComplex>(a))
216     #define PetscPowComplex(a, b)        petsccomplexlib::pow(static_cast<PetscComplex>(a), static_cast<PetscComplex>(b))
217     #define PetscExpComplex(a)           petsccomplexlib::exp(static_cast<PetscComplex>(a))
218     #define PetscLogComplex(a)           petsccomplexlib::log(static_cast<PetscComplex>(a))
219     #define PetscSinComplex(a)           petsccomplexlib::sin(static_cast<PetscComplex>(a))
220     #define PetscCosComplex(a)           petsccomplexlib::cos(static_cast<PetscComplex>(a))
221     #define PetscTanComplex(a)           petsccomplexlib::tan(static_cast<PetscComplex>(a))
222     #define PetscAsinComplex(a)          petsccomplexlib::asin(static_cast<PetscComplex>(a))
223     #define PetscAcosComplex(a)          petsccomplexlib::acos(static_cast<PetscComplex>(a))
224     #define PetscAtanComplex(a)          petsccomplexlib::atan(static_cast<PetscComplex>(a))
225     #define PetscSinhComplex(a)          petsccomplexlib::sinh(static_cast<PetscComplex>(a))
226     #define PetscCoshComplex(a)          petsccomplexlib::cosh(static_cast<PetscComplex>(a))
227     #define PetscTanhComplex(a)          petsccomplexlib::tanh(static_cast<PetscComplex>(a))
228     #define PetscAsinhComplex(a)         petsccomplexlib::asinh(static_cast<PetscComplex>(a))
229     #define PetscAcoshComplex(a)         petsccomplexlib::acosh(static_cast<PetscComplex>(a))
230     #define PetscAtanhComplex(a)         petsccomplexlib::atanh(static_cast<PetscComplex>(a))
231 
232   /* TODO: Add configure tests
233 
234 #if !defined(PETSC_HAVE_CXX_TAN_COMPLEX)
235 #undef PetscTanComplex
236 static inline PetscComplex PetscTanComplex(PetscComplex z)
237 {
238   return PetscSinComplex(z)/PetscCosComplex(z);
239 }
240 #endif
241 
242 #if !defined(PETSC_HAVE_CXX_TANH_COMPLEX)
243 #undef PetscTanhComplex
244 static inline PetscComplex PetscTanhComplex(PetscComplex z)
245 {
246   return PetscSinhComplex(z)/PetscCoshComplex(z);
247 }
248 #endif
249 
250 #if !defined(PETSC_HAVE_CXX_ASIN_COMPLEX)
251 #undef PetscAsinComplex
252 static inline PetscComplex PetscAsinComplex(PetscComplex z)
253 {
254   const PetscComplex j(0,1);
255   return -j*PetscLogComplex(j*z+PetscSqrtComplex(1.0f-z*z));
256 }
257 #endif
258 
259 #if !defined(PETSC_HAVE_CXX_ACOS_COMPLEX)
260 #undef PetscAcosComplex
261 static inline PetscComplex PetscAcosComplex(PetscComplex z)
262 {
263   const PetscComplex j(0,1);
264   return j*PetscLogComplex(z-j*PetscSqrtComplex(1.0f-z*z));
265 }
266 #endif
267 
268 #if !defined(PETSC_HAVE_CXX_ATAN_COMPLEX)
269 #undef PetscAtanComplex
270 static inline PetscComplex PetscAtanComplex(PetscComplex z)
271 {
272   const PetscComplex j(0,1);
273   return 0.5f*j*PetscLogComplex((1.0f-j*z)/(1.0f+j*z));
274 }
275 #endif
276 
277 #if !defined(PETSC_HAVE_CXX_ASINH_COMPLEX)
278 #undef PetscAsinhComplex
279 static inline PetscComplex PetscAsinhComplex(PetscComplex z)
280 {
281   return PetscLogComplex(z+PetscSqrtComplex(z*z+1.0f));
282 }
283 #endif
284 
285 #if !defined(PETSC_HAVE_CXX_ACOSH_COMPLEX)
286 #undef PetscAcoshComplex
287 static inline PetscComplex PetscAcoshComplex(PetscComplex z)
288 {
289   return PetscLogComplex(z+PetscSqrtComplex(z*z-1.0f));
290 }
291 #endif
292 
293 #if !defined(PETSC_HAVE_CXX_ATANH_COMPLEX)
294 #undef PetscAtanhComplex
295 static inline PetscComplex PetscAtanhComplex(PetscComplex z)
296 {
297   return 0.5f*PetscLogComplex((1.0f+z)/(1.0f-z));
298 }
299 #endif
300 
301 */
302 
303   #else /* C99 support of complex number */
304 
305     #if defined(PETSC_USE_REAL_SINGLE)
306       #define PetscRealPartComplex(a)      crealf(a)
307       #define PetscImaginaryPartComplex(a) cimagf(a)
308       #define PetscAbsComplex(a)           cabsf(a)
309       #define PetscArgComplex(a)           cargf(a)
310       #define PetscConjComplex(a)          conjf(a)
311       #define PetscSqrtComplex(a)          csqrtf(a)
312       #define PetscPowComplex(a, b)        cpowf(a, b)
313       #define PetscExpComplex(a)           cexpf(a)
314       #define PetscLogComplex(a)           clogf(a)
315       #define PetscSinComplex(a)           csinf(a)
316       #define PetscCosComplex(a)           ccosf(a)
317       #define PetscTanComplex(a)           ctanf(a)
318       #define PetscAsinComplex(a)          casinf(a)
319       #define PetscAcosComplex(a)          cacosf(a)
320       #define PetscAtanComplex(a)          catanf(a)
321       #define PetscSinhComplex(a)          csinhf(a)
322       #define PetscCoshComplex(a)          ccoshf(a)
323       #define PetscTanhComplex(a)          ctanhf(a)
324       #define PetscAsinhComplex(a)         casinhf(a)
325       #define PetscAcoshComplex(a)         cacoshf(a)
326       #define PetscAtanhComplex(a)         catanhf(a)
327 
328     #elif defined(PETSC_USE_REAL_DOUBLE)
329       #define PetscRealPartComplex(a)      creal(a)
330       #define PetscImaginaryPartComplex(a) cimag(a)
331       #define PetscAbsComplex(a)           cabs(a)
332       #define PetscArgComplex(a)           carg(a)
333       #define PetscConjComplex(a)          conj(a)
334       #define PetscSqrtComplex(a)          csqrt(a)
335       #define PetscPowComplex(a, b)        cpow(a, b)
336       #define PetscExpComplex(a)           cexp(a)
337       #define PetscLogComplex(a)           clog(a)
338       #define PetscSinComplex(a)           csin(a)
339       #define PetscCosComplex(a)           ccos(a)
340       #define PetscTanComplex(a)           ctan(a)
341       #define PetscAsinComplex(a)          casin(a)
342       #define PetscAcosComplex(a)          cacos(a)
343       #define PetscAtanComplex(a)          catan(a)
344       #define PetscSinhComplex(a)          csinh(a)
345       #define PetscCoshComplex(a)          ccosh(a)
346       #define PetscTanhComplex(a)          ctanh(a)
347       #define PetscAsinhComplex(a)         casinh(a)
348       #define PetscAcoshComplex(a)         cacosh(a)
349       #define PetscAtanhComplex(a)         catanh(a)
350 
351     #elif defined(PETSC_USE_REAL___FLOAT128)
352       #define PetscRealPartComplex(a)      crealq(a)
353       #define PetscImaginaryPartComplex(a) cimagq(a)
354       #define PetscAbsComplex(a)           cabsq(a)
355       #define PetscArgComplex(a)           cargq(a)
356       #define PetscConjComplex(a)          conjq(a)
357       #define PetscSqrtComplex(a)          csqrtq(a)
358       #define PetscPowComplex(a, b)        cpowq(a, b)
359       #define PetscExpComplex(a)           cexpq(a)
360       #define PetscLogComplex(a)           clogq(a)
361       #define PetscSinComplex(a)           csinq(a)
362       #define PetscCosComplex(a)           ccosq(a)
363       #define PetscTanComplex(a)           ctanq(a)
364       #define PetscAsinComplex(a)          casinq(a)
365       #define PetscAcosComplex(a)          cacosq(a)
366       #define PetscAtanComplex(a)          catanq(a)
367       #define PetscSinhComplex(a)          csinhq(a)
368       #define PetscCoshComplex(a)          ccoshq(a)
369       #define PetscTanhComplex(a)          ctanhq(a)
370       #define PetscAsinhComplex(a)         casinhq(a)
371       #define PetscAcoshComplex(a)         cacoshq(a)
372       #define PetscAtanhComplex(a)         catanhq(a)
373 
374     #endif /* PETSC_USE_REAL_* */
375   #endif   /* (__cplusplus) */
376 
377 /*MC
378     PETSC_i - the pure imaginary complex number i
379 
380    Level: intermediate
381 
382 .seealso: `PetscComplex`, `PetscScalar`
383 M*/
384 PETSC_EXTERN PetscComplex PETSC_i;
385 
386 /*
387    Try to do the right thing for complex number construction: see
388    http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1464.htm
389    for details
390 */
391 static inline PetscComplex PetscCMPLX(PetscReal x, PetscReal y)
392 {
393   #if defined(__cplusplus) && !defined(PETSC_USE_REAL___FLOAT128)
394   return PetscComplex(x, y);
395   #elif defined(_Imaginary_I)
396   return x + y * _Imaginary_I;
397   #else
398   { /* In both C99 and C11 (ISO/IEC 9899, Section 6.2.5),
399 
400        "For each floating type there is a corresponding real type, which is always a real floating
401        type. For real floating types, it is the same type. For complex types, it is the type given
402        by deleting the keyword _Complex from the type name."
403 
404        So type punning should be portable. */
405     union
406     {
407       PetscComplex z;
408       PetscReal    f[2];
409     } uz;
410 
411     uz.f[0] = x;
412     uz.f[1] = y;
413     return uz.z;
414   }
415   #endif
416 }
417 
418   #define MPIU_C_COMPLEX        MPI_C_COMPLEX PETSC_DEPRECATED_MACRO(3, 15, 0, "MPI_C_COMPLEX", )
419   #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX PETSC_DEPRECATED_MACRO(3, 15, 0, "MPI_C_DOUBLE_COMPLEX", )
420 
421   #if defined(PETSC_HAVE_REAL___FLOAT128) && !defined(PETSC_SKIP_REAL___FLOAT128)
422     // if complex is not used, then quadmath.h won't be included by petscsystypes.h
423     #if defined(PETSC_USE_COMPLEX)
424       #define MPIU___COMPLEX128_ATTR_TAG PETSC_ATTRIBUTE_MPI_TYPE_TAG(__complex128)
425     #else
426       #define MPIU___COMPLEX128_ATTR_TAG
427     #endif
428 
429 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 MPIU___COMPLEX128_ATTR_TAG;
430 
431     #undef MPIU___COMPLEX128_ATTR_TAG
432   #endif /* PETSC_HAVE_REAL___FLOAT128 */
433 
434   /*MC
435    MPIU_COMPLEX - Portable MPI datatype corresponding to `PetscComplex` independent of the precision of `PetscComplex`
436 
437    Level: beginner
438 
439    Note:
440    In MPI calls that require an MPI datatype that matches a `PetscComplex` or array of `PetscComplex` values, pass this value.
441 
442 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_REAL`, `MPIU_SCALAR`, `MPIU_COMPLEX`, `MPIU_INT`, `PETSC_i`
443 M*/
444   #if defined(PETSC_USE_REAL_SINGLE)
445     #define MPIU_COMPLEX MPI_C_COMPLEX
446   #elif defined(PETSC_USE_REAL_DOUBLE)
447     #define MPIU_COMPLEX MPI_C_DOUBLE_COMPLEX
448   #elif defined(PETSC_USE_REAL___FLOAT128)
449     #define MPIU_COMPLEX MPIU___COMPLEX128
450   #elif defined(PETSC_USE_REAL___FP16)
451     #define MPIU_COMPLEX MPI_C_COMPLEX
452   #endif /* PETSC_USE_REAL_* */
453 
454 #endif /* PETSC_HAVE_COMPLEX */
455 
456 /*
457     Scalar number definitions
458  */
459 #if defined(PETSC_USE_COMPLEX) && defined(PETSC_HAVE_COMPLEX)
460   /*MC
461    MPIU_SCALAR - Portable MPI datatype corresponding to `PetscScalar` independent of the precision of `PetscScalar`
462 
463    Level: beginner
464 
465    Note:
466    In MPI calls that require an MPI datatype that matches a `PetscScalar` or array of `PetscScalar` values, pass this value.
467 
468 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_REAL`, `MPIU_COMPLEX`, `MPIU_INT`
469 M*/
470   #define MPIU_SCALAR MPIU_COMPLEX
471 
472   /*MC
473    PetscRealPart - Returns the real part of a `PetscScalar`
474 
475    Synopsis:
476    #include <petscmath.h>
477    PetscReal PetscRealPart(PetscScalar v)
478 
479    Not Collective
480 
481    Input Parameter:
482 .  v - value to find the real part of
483 
484    Level: beginner
485 
486 .seealso: `PetscScalar`, `PetscImaginaryPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
487 M*/
488   #define PetscRealPart(a) PetscRealPartComplex(a)
489 
490   /*MC
491    PetscImaginaryPart - Returns the imaginary part of a `PetscScalar`
492 
493    Synopsis:
494    #include <petscmath.h>
495    PetscReal PetscImaginaryPart(PetscScalar v)
496 
497    Not Collective
498 
499    Input Parameter:
500 .  v - value to find the imaginary part of
501 
502    Level: beginner
503 
504    Note:
505        If PETSc was configured for real numbers then this always returns the value 0
506 
507 .seealso: `PetscScalar`, `PetscRealPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
508 M*/
509   #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a)
510 
511   #define PetscAbsScalar(a)    PetscAbsComplex(a)
512   #define PetscArgScalar(a)    PetscArgComplex(a)
513   #define PetscConj(a)         PetscConjComplex(a)
514   #define PetscSqrtScalar(a)   PetscSqrtComplex(a)
515   #define PetscPowScalar(a, b) PetscPowComplex(a, b)
516   #define PetscExpScalar(a)    PetscExpComplex(a)
517   #define PetscLogScalar(a)    PetscLogComplex(a)
518   #define PetscSinScalar(a)    PetscSinComplex(a)
519   #define PetscCosScalar(a)    PetscCosComplex(a)
520   #define PetscTanScalar(a)    PetscTanComplex(a)
521   #define PetscAsinScalar(a)   PetscAsinComplex(a)
522   #define PetscAcosScalar(a)   PetscAcosComplex(a)
523   #define PetscAtanScalar(a)   PetscAtanComplex(a)
524   #define PetscSinhScalar(a)   PetscSinhComplex(a)
525   #define PetscCoshScalar(a)   PetscCoshComplex(a)
526   #define PetscTanhScalar(a)   PetscTanhComplex(a)
527   #define PetscAsinhScalar(a)  PetscAsinhComplex(a)
528   #define PetscAcoshScalar(a)  PetscAcoshComplex(a)
529   #define PetscAtanhScalar(a)  PetscAtanhComplex(a)
530 
531 #else /* PETSC_USE_COMPLEX */
532   #define MPIU_SCALAR           MPIU_REAL
533   #define PetscRealPart(a)      (a)
534   #define PetscImaginaryPart(a) ((PetscReal)0)
535   #define PetscAbsScalar(a)     PetscAbsReal(a)
536   #define PetscArgScalar(a)     (((a) < (PetscReal)0) ? PETSC_PI : (PetscReal)0)
537   #define PetscConj(a)          (a)
538   #define PetscSqrtScalar(a)    PetscSqrtReal(a)
539   #define PetscPowScalar(a, b)  PetscPowReal(a, b)
540   #define PetscExpScalar(a)     PetscExpReal(a)
541   #define PetscLogScalar(a)     PetscLogReal(a)
542   #define PetscSinScalar(a)     PetscSinReal(a)
543   #define PetscCosScalar(a)     PetscCosReal(a)
544   #define PetscTanScalar(a)     PetscTanReal(a)
545   #define PetscAsinScalar(a)    PetscAsinReal(a)
546   #define PetscAcosScalar(a)    PetscAcosReal(a)
547   #define PetscAtanScalar(a)    PetscAtanReal(a)
548   #define PetscSinhScalar(a)    PetscSinhReal(a)
549   #define PetscCoshScalar(a)    PetscCoshReal(a)
550   #define PetscTanhScalar(a)    PetscTanhReal(a)
551   #define PetscAsinhScalar(a)   PetscAsinhReal(a)
552   #define PetscAcoshScalar(a)   PetscAcoshReal(a)
553   #define PetscAtanhScalar(a)   PetscAtanhReal(a)
554 
555 #endif /* PETSC_USE_COMPLEX */
556 
557 /*
558    Certain objects may be created using either single or double precision.
559    This is currently not used.
560 */
561 typedef enum {
562   PETSC_SCALAR_DOUBLE,
563   PETSC_SCALAR_SINGLE,
564   PETSC_SCALAR_LONG_DOUBLE,
565   PETSC_SCALAR_HALF
566 } PetscScalarPrecision;
567 
568 /*MC
569    PetscAbs - Returns the absolute value of a number
570 
571    Synopsis:
572    #include <petscmath.h>
573    type PetscAbs(type v)
574 
575    Not Collective
576 
577    Input Parameter:
578 .  v - the number
579 
580    Level: beginner
581 
582    Note:
583    The type can be integer or real floating point value, but cannot be complex
584 
585 .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`, `PetscSign()`
586 M*/
587 #define PetscAbs(a) (((a) >= 0) ? (a) : (-(a)))
588 
589 /*MC
590    PetscSign - Returns the sign of a number as an integer of value -1, 0, or 1
591 
592    Synopsis:
593    #include <petscmath.h>
594    int PetscSign(type v)
595 
596    Not Collective
597 
598    Input Parameter:
599 .  v - the number
600 
601    Level: beginner
602 
603    Note:
604    The type can be integer or real floating point value
605 
606 .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`
607 M*/
608 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1)
609 
610 /*MC
611    PetscMin - Returns minimum of two numbers
612 
613    Synopsis:
614    #include <petscmath.h>
615    type PetscMin(type v1,type v2)
616 
617    Not Collective
618 
619    Input Parameters:
620 +  v1 - first value to find minimum of
621 -  v2 - second value to find minimum of
622 
623    Level: beginner
624 
625    Note:
626    The type can be integer or floating point value, but cannot be complex
627 
628 .seealso: `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
629 M*/
630 #define PetscMin(a, b) (((a) < (b)) ? (a) : (b))
631 
632 /*MC
633    PetscMax - Returns maximum of two numbers
634 
635    Synopsis:
636    #include <petscmath.h>
637    type max PetscMax(type v1,type v2)
638 
639    Not Collective
640 
641    Input Parameters:
642 +  v1 - first value to find maximum of
643 -  v2 - second value to find maximum of
644 
645    Level: beginner
646 
647    Note:
648    The type can be integer or floating point value
649 
650 .seealso: `PetscMin()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
651 M*/
652 #define PetscMax(a, b) (((a) < (b)) ? (b) : (a))
653 
654 /*MC
655    PetscClipInterval - Returns a number clipped to be within an interval
656 
657    Synopsis:
658    #include <petscmath.h>
659    type clip PetscClipInterval(type x,type a,type b)
660 
661    Not Collective
662 
663    Input Parameters:
664 +  x - value to use if within interval [a,b]
665 .  a - lower end of interval
666 -  b - upper end of interval
667 
668    Level: beginner
669 
670    Note:
671    The type can be integer or floating point value
672 
673    Example\:
674 .vb
675   PetscInt c = PetscClipInterval(5, 2, 3); // the value of c is 3
676   PetscInt c = PetscClipInterval(5, 2, 6); // the value of c is 5
677 .ve
678 
679 .seealso: `PetscMin()`, `PetscMax()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
680 M*/
681 #define PetscClipInterval(x, a, b) (PetscMax((a), PetscMin((x), (b))))
682 
683 /*MC
684    PetscAbsInt - Returns the absolute value of an integer
685 
686    Synopsis:
687    #include <petscmath.h>
688    int abs PetscAbsInt(int v1)
689 
690    Input Parameter:
691 .   v1 - the integer
692 
693    Level: beginner
694 
695 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsReal()`, `PetscSqr()`
696 M*/
697 #define PetscAbsInt(a) (((a) < 0) ? (-(a)) : (a))
698 
699 /*MC
700    PetscAbsReal - Returns the absolute value of an real number
701 
702    Synopsis:
703    #include <petscmath.h>
704    Real abs PetscAbsReal(PetscReal v1)
705 
706    Input Parameter:
707 .   v1 - the `PetscReal` value
708 
709    Level: beginner
710 
711 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscSqr()`
712 M*/
713 #if defined(PETSC_USE_REAL_SINGLE)
714   #define PetscAbsReal(a) fabsf(a)
715 #elif defined(PETSC_USE_REAL_DOUBLE)
716   #define PetscAbsReal(a) fabs(a)
717 #elif defined(PETSC_USE_REAL___FLOAT128)
718   #define PetscAbsReal(a) fabsq(a)
719 #elif defined(PETSC_USE_REAL___FP16)
720   #define PetscAbsReal(a) fabsf(a)
721 #endif
722 
723 /*MC
724    PetscSqr - Returns the square of a number
725 
726    Synopsis:
727    #include <petscmath.h>
728    type sqr PetscSqr(type v1)
729 
730    Not Collective
731 
732    Input Parameter:
733 .   v1 - the value
734 
735    Level: beginner
736 
737    Note:
738    The type can be integer, floating point, or complex floating point
739 
740 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`
741 M*/
742 #define PetscSqr(a) ((a) * (a))
743 
744 #if defined(PETSC_USE_REAL_SINGLE)
745   #define PetscRealConstant(constant) constant##F
746 #elif defined(PETSC_USE_REAL_DOUBLE)
747   #define PetscRealConstant(constant) constant
748 #elif defined(PETSC_USE_REAL___FLOAT128)
749   #define PetscRealConstant(constant) constant##Q
750 #elif defined(PETSC_USE_REAL___FP16)
751   #define PetscRealConstant(constant) constant##F
752 #endif
753 
754 /*
755      Basic constants
756 */
757 #define PETSC_PI    PetscRealConstant(3.1415926535897932384626433832795029)
758 #define PETSC_PHI   PetscRealConstant(1.6180339887498948482045868343656381)
759 #define PETSC_SQRT2 PetscRealConstant(1.4142135623730950488016887242096981)
760 
761 #if defined(PETSC_USE_REAL_SINGLE)
762   #define PETSC_MAX_REAL             3.40282346638528860e+38F
763   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
764   #define PETSC_REAL_MIN             1.1754944e-38F
765   #define PETSC_MACHINE_EPSILON      1.19209290e-07F
766   #define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F
767   #define PETSC_SMALL                1.e-5F
768 #elif defined(PETSC_USE_REAL_DOUBLE)
769   #define PETSC_MAX_REAL             1.7976931348623157e+308
770   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
771   #define PETSC_REAL_MIN             2.225073858507201e-308
772   #define PETSC_MACHINE_EPSILON      2.2204460492503131e-16
773   #define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08
774   #define PETSC_SMALL                1.e-10
775 #elif defined(PETSC_USE_REAL___FLOAT128)
776   #define PETSC_MAX_REAL             FLT128_MAX
777   #define PETSC_MIN_REAL             (-FLT128_MAX)
778   #define PETSC_REAL_MIN             FLT128_MIN
779   #define PETSC_MACHINE_EPSILON      FLT128_EPSILON
780   #define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q
781   #define PETSC_SMALL                1.e-20Q
782 #elif defined(PETSC_USE_REAL___FP16)
783   #define PETSC_MAX_REAL             65504.0F
784   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
785   #define PETSC_REAL_MIN             .00006103515625F
786   #define PETSC_MACHINE_EPSILON      .0009765625F
787   #define PETSC_SQRT_MACHINE_EPSILON .03125F
788   #define PETSC_SMALL                5.e-3F
789 #endif
790 
791 /*MC
792     PETSC_INFINITY - a finite number that represents infinity for setting certain bounds in `Tao`
793 
794    Level: intermediate
795 
796   Note:
797   This is not the IEEE infinity value
798 
799 .seealso: `PETSC_NINFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
800 M*/
801 #define PETSC_INFINITY (PETSC_MAX_REAL / 4)
802 
803 /*MC
804     PETSC_NINFINITY - a finite number that represents negative infinity for setting certain bounds in `Tao`
805 
806    Level: intermediate
807 
808   Note:
809   This is not the negative IEEE infinity value
810 
811 .seealso: `PETSC_INFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
812 M*/
813 #define PETSC_NINFINITY (-PETSC_INFINITY)
814 
815 PETSC_EXTERN PetscBool  PetscIsInfReal(PetscReal);
816 PETSC_EXTERN PetscBool  PetscIsNanReal(PetscReal);
817 PETSC_EXTERN PetscBool  PetscIsNormalReal(PetscReal);
818 static inline PetscBool PetscIsInfOrNanReal(PetscReal v)
819 {
820   return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;
821 }
822 static inline PetscBool PetscIsInfScalar(PetscScalar v)
823 {
824   return PetscIsInfReal(PetscAbsScalar(v));
825 }
826 static inline PetscBool PetscIsNanScalar(PetscScalar v)
827 {
828   return PetscIsNanReal(PetscAbsScalar(v));
829 }
830 static inline PetscBool PetscIsInfOrNanScalar(PetscScalar v)
831 {
832   return PetscIsInfOrNanReal(PetscAbsScalar(v));
833 }
834 static inline PetscBool PetscIsNormalScalar(PetscScalar v)
835 {
836   return PetscIsNormalReal(PetscAbsScalar(v));
837 }
838 
839 PETSC_EXTERN PetscBool PetscIsCloseAtTol(PetscReal, PetscReal, PetscReal, PetscReal);
840 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal, PetscReal);
841 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar, PetscScalar);
842 
843 /*@C
844   PetscIsCloseAtTolScalar - Like `PetscIsCloseAtTol()` but for `PetscScalar`
845 
846   Input Parameters:
847 + lhs  - The first number
848 . rhs  - The second number
849 . rtol - The relative tolerance
850 - atol - The absolute tolerance
851 
852   Level: beginner
853 
854   Note:
855   This routine is equivalent to `PetscIsCloseAtTol()` when PETSc is configured without complex
856   numbers.
857 
858 .seealso: `PetscIsCloseAtTol()`
859 @*/
860 static inline PetscBool PetscIsCloseAtTolScalar(PetscScalar lhs, PetscScalar rhs, PetscReal rtol, PetscReal atol)
861 {
862   PetscBool close = PetscIsCloseAtTol(PetscRealPart(lhs), PetscRealPart(rhs), rtol, atol);
863 
864   if (PetscDefined(USE_COMPLEX)) close = (PetscBool)(close && PetscIsCloseAtTol(PetscImaginaryPart(lhs), PetscImaginaryPart(rhs), rtol, atol));
865   return close;
866 }
867 
868 /*
869     These macros are currently hardwired to match the regular data types, so there is no support for a different
870     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
871  */
872 #define MPIU_MATSCALAR MPIU_SCALAR
873 typedef PetscScalar MatScalar;
874 typedef PetscReal   MatReal;
875 
876 struct petsc_mpiu_2scalar {
877   PetscScalar a, b;
878 };
879 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2scalar);
880 
881 /* MPI Datatypes for composite reductions */
882 struct petsc_mpiu_real_int {
883   PetscReal v;
884   PetscInt  i;
885 };
886 
887 struct petsc_mpiu_scalar_int {
888   PetscScalar v;
889   PetscInt    i;
890 };
891 
892 PETSC_EXTERN MPI_Datatype MPIU_REAL_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_real_int);
893 PETSC_EXTERN MPI_Datatype MPIU_SCALAR_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_scalar_int);
894 
895 #if defined(PETSC_USE_64BIT_INDICES)
896 struct /* __attribute__((packed, aligned(alignof(PetscInt *)))) */ petsc_mpiu_2int {
897   PetscInt a;
898   PetscInt b;
899 };
900 /*
901  static_assert(sizeof(struct petsc_mpiu_2int) == 2 * sizeof(PetscInt), "");
902  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt *), "");
903  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt[2]), "");
904 
905  clang generates warnings that petsc_mpiu_2int is not layout compatible with PetscInt[2] or
906  PetscInt *, even though (with everything else uncommented) both of the static_asserts above
907  pass! So we just comment it out...
908 */
909 PETSC_EXTERN MPI_Datatype MPIU_2INT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2int) */;
910 #else
911   #define MPIU_2INT MPI_2INT
912 #endif
913 PETSC_EXTERN MPI_Datatype MPI_4INT;
914 PETSC_EXTERN MPI_Datatype MPIU_4INT;
915 
916 static inline PetscInt PetscPowInt(PetscInt base, PetscInt power)
917 {
918   PetscInt result = 1;
919   while (power) {
920     if (power & 1) result *= base;
921     power >>= 1;
922     if (power) base *= base;
923   }
924   return result;
925 }
926 
927 static inline PetscInt64 PetscPowInt64(PetscInt base, PetscInt power)
928 {
929   PetscInt64 result = 1;
930   while (power) {
931     if (power & 1) result *= base;
932     power >>= 1;
933     if (power) base *= base;
934   }
935   return result;
936 }
937 
938 static inline PetscReal PetscPowRealInt(PetscReal base, PetscInt power)
939 {
940   PetscReal result = 1;
941   if (power < 0) {
942     power = -power;
943     base  = ((PetscReal)1) / base;
944   }
945   while (power) {
946     if (power & 1) result *= base;
947     power >>= 1;
948     if (power) base *= base;
949   }
950   return result;
951 }
952 
953 static inline PetscScalar PetscPowScalarInt(PetscScalar base, PetscInt power)
954 {
955   PetscScalar result = (PetscReal)1;
956   if (power < 0) {
957     power = -power;
958     base  = ((PetscReal)1) / base;
959   }
960   while (power) {
961     if (power & 1) result *= base;
962     power >>= 1;
963     if (power) base *= base;
964   }
965   return result;
966 }
967 
968 static inline PetscScalar PetscPowScalarReal(PetscScalar base, PetscReal power)
969 {
970   PetscScalar cpower = power;
971   return PetscPowScalar(base, cpower);
972 }
973 
974 /*MC
975     PetscApproximateLTE - Performs a less than or equal to on a given constant with a fudge for floating point numbers
976 
977    Synopsis:
978    #include <petscmath.h>
979    bool PetscApproximateLTE(PetscReal x,constant float)
980 
981    Not Collective
982 
983    Input Parameters:
984 +   x - the variable
985 -   b - the constant float it is checking if `x` is less than or equal to
986 
987    Level: advanced
988 
989    Notes:
990      The fudge factor is the value `PETSC_SMALL`
991 
992      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
993 
994      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
995      floating point results.
996 
997   Example\:
998 .vb
999   PetscReal x;
1000   if (PetscApproximateLTE(x, 3.2)) { // replaces if (x <= 3.2) {
1001 .ve
1002 
1003 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateGTE()`
1004 M*/
1005 #define PetscApproximateLTE(x, b) ((x) <= (PetscRealConstant(b) + PETSC_SMALL))
1006 
1007 /*MC
1008     PetscApproximateGTE - Performs a greater than or equal to on a given constant with a fudge for floating point numbers
1009 
1010    Synopsis:
1011    #include <petscmath.h>
1012    bool PetscApproximateGTE(PetscReal x,constant float)
1013 
1014    Not Collective
1015 
1016    Input Parameters:
1017 +   x - the variable
1018 -   b - the constant float it is checking if `x` is greater than or equal to
1019 
1020    Level: advanced
1021 
1022    Notes:
1023      The fudge factor is the value `PETSC_SMALL`
1024 
1025      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
1026 
1027      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1028      floating point results.
1029 
1030   Example\:
1031 .vb
1032   PetscReal x;
1033   if (PetscApproximateGTE(x, 3.2)) {  // replaces if (x >= 3.2) {
1034 .ve
1035 
1036 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1037 M*/
1038 #define PetscApproximateGTE(x, b) ((x) >= (PetscRealConstant(b) - PETSC_SMALL))
1039 
1040 /*MC
1041     PetscCeilInt - Returns the ceiling of the quotation of two positive integers
1042 
1043    Synopsis:
1044    #include <petscmath.h>
1045    PetscInt PetscCeilInt(PetscInt x,PetscInt y)
1046 
1047    Not Collective
1048 
1049    Input Parameters:
1050 +   x - the numerator
1051 -   y - the denominator
1052 
1053    Level: advanced
1054 
1055   Example\:
1056 .vb
1057   PetscInt n = PetscCeilInt(10, 3); // n has the value of 4
1058 .ve
1059 
1060 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1061 M*/
1062 #define PetscCeilInt(x, y) ((((PetscInt)(x)) / ((PetscInt)(y))) + ((((PetscInt)(x)) % ((PetscInt)(y))) ? 1 : 0))
1063 
1064 #define PetscCeilInt64(x, y) ((((PetscInt64)(x)) / ((PetscInt64)(y))) + ((((PetscInt64)(x)) % ((PetscInt64)(y))) ? 1 : 0))
1065 
1066 PETSC_EXTERN PetscErrorCode PetscLinearRegression(PetscInt, const PetscReal[], const PetscReal[], PetscReal *, PetscReal *);
1067