xref: /petsc/include/petscmath.h (revision 945b2ebdc5b7ebcb622acfbe60dbc4e0db537728)
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 #pragma once
11 
12 #include <math.h>
13 #include <petscmacros.h>
14 #include <petscsystypes.h>
15 
16 /* SUBMANSEC = Sys */
17 
18 /*
19 
20    Defines operations that are different for complex and real numbers.
21    All PETSc objects in one program are built around the object
22    PetscScalar which is either always a real or a complex.
23 
24 */
25 
26 /*
27     Real number definitions
28  */
29 #if defined(PETSC_USE_REAL_SINGLE)
30   #define PetscSqrtReal(a)        sqrtf(a)
31   #define PetscCbrtReal(a)        cbrtf(a)
32   #define PetscHypotReal(a, b)    hypotf(a, b)
33   #define PetscAtan2Real(a, b)    atan2f(a, b)
34   #define PetscPowReal(a, b)      powf(a, b)
35   #define PetscExpReal(a)         expf(a)
36   #define PetscLogReal(a)         logf(a)
37   #define PetscLog10Real(a)       log10f(a)
38   #define PetscLog2Real(a)        log2f(a)
39   #define PetscSinReal(a)         sinf(a)
40   #define PetscCosReal(a)         cosf(a)
41   #define PetscTanReal(a)         tanf(a)
42   #define PetscAsinReal(a)        asinf(a)
43   #define PetscAcosReal(a)        acosf(a)
44   #define PetscAtanReal(a)        atanf(a)
45   #define PetscSinhReal(a)        sinhf(a)
46   #define PetscCoshReal(a)        coshf(a)
47   #define PetscTanhReal(a)        tanhf(a)
48   #define PetscAsinhReal(a)       asinhf(a)
49   #define PetscAcoshReal(a)       acoshf(a)
50   #define PetscAtanhReal(a)       atanhf(a)
51   #define PetscErfReal(a)         erff(a)
52   #define PetscCeilReal(a)        ceilf(a)
53   #define PetscFloorReal(a)       floorf(a)
54   #define PetscFmodReal(a, b)     fmodf(a, b)
55   #define PetscCopysignReal(a, b) copysignf(a, b)
56   #define PetscTGamma(a)          tgammaf(a)
57   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
58     #define PetscLGamma(a) gammaf(a)
59   #else
60     #define PetscLGamma(a) lgammaf(a)
61   #endif
62 
63 #elif defined(PETSC_USE_REAL_DOUBLE)
64   #define PetscSqrtReal(a)        sqrt(a)
65   #define PetscCbrtReal(a)        cbrt(a)
66   #define PetscHypotReal(a, b)    hypot(a, b)
67   #define PetscAtan2Real(a, b)    atan2(a, b)
68   #define PetscPowReal(a, b)      pow(a, b)
69   #define PetscExpReal(a)         exp(a)
70   #define PetscLogReal(a)         log(a)
71   #define PetscLog10Real(a)       log10(a)
72   #define PetscLog2Real(a)        log2(a)
73   #define PetscSinReal(a)         sin(a)
74   #define PetscCosReal(a)         cos(a)
75   #define PetscTanReal(a)         tan(a)
76   #define PetscAsinReal(a)        asin(a)
77   #define PetscAcosReal(a)        acos(a)
78   #define PetscAtanReal(a)        atan(a)
79   #define PetscSinhReal(a)        sinh(a)
80   #define PetscCoshReal(a)        cosh(a)
81   #define PetscTanhReal(a)        tanh(a)
82   #define PetscAsinhReal(a)       asinh(a)
83   #define PetscAcoshReal(a)       acosh(a)
84   #define PetscAtanhReal(a)       atanh(a)
85   #define PetscErfReal(a)         erf(a)
86   #define PetscCeilReal(a)        ceil(a)
87   #define PetscFloorReal(a)       floor(a)
88   #define PetscFmodReal(a, b)     fmod(a, b)
89   #define PetscCopysignReal(a, b) copysign(a, b)
90   #define PetscTGamma(a)          tgamma(a)
91   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
92     #define PetscLGamma(a) gamma(a)
93   #else
94     #define PetscLGamma(a) lgamma(a)
95   #endif
96 
97 #elif defined(PETSC_USE_REAL___FLOAT128)
98   #define PetscSqrtReal(a)        sqrtq(a)
99   #define PetscCbrtReal(a)        cbrtq(a)
100   #define PetscHypotReal(a, b)    hypotq(a, b)
101   #define PetscAtan2Real(a, b)    atan2q(a, b)
102   #define PetscPowReal(a, b)      powq(a, b)
103   #define PetscExpReal(a)         expq(a)
104   #define PetscLogReal(a)         logq(a)
105   #define PetscLog10Real(a)       log10q(a)
106   #define PetscLog2Real(a)        log2q(a)
107   #define PetscSinReal(a)         sinq(a)
108   #define PetscCosReal(a)         cosq(a)
109   #define PetscTanReal(a)         tanq(a)
110   #define PetscAsinReal(a)        asinq(a)
111   #define PetscAcosReal(a)        acosq(a)
112   #define PetscAtanReal(a)        atanq(a)
113   #define PetscSinhReal(a)        sinhq(a)
114   #define PetscCoshReal(a)        coshq(a)
115   #define PetscTanhReal(a)        tanhq(a)
116   #define PetscAsinhReal(a)       asinhq(a)
117   #define PetscAcoshReal(a)       acoshq(a)
118   #define PetscAtanhReal(a)       atanhq(a)
119   #define PetscErfReal(a)         erfq(a)
120   #define PetscCeilReal(a)        ceilq(a)
121   #define PetscFloorReal(a)       floorq(a)
122   #define PetscFmodReal(a, b)     fmodq(a, b)
123   #define PetscCopysignReal(a, b) copysignq(a, b)
124   #define PetscTGamma(a)          tgammaq(a)
125   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
126     #define PetscLGamma(a) gammaq(a)
127   #else
128     #define PetscLGamma(a) lgammaq(a)
129   #endif
130 
131 #elif defined(PETSC_USE_REAL___FP16)
132   #define PetscSqrtReal(a)        sqrtf(a)
133   #define PetscCbrtReal(a)        cbrtf(a)
134   #define PetscHypotReal(a, b)    hypotf(a, b)
135   #define PetscAtan2Real(a, b)    atan2f(a, b)
136   #define PetscPowReal(a, b)      powf(a, b)
137   #define PetscExpReal(a)         expf(a)
138   #define PetscLogReal(a)         logf(a)
139   #define PetscLog10Real(a)       log10f(a)
140   #define PetscLog2Real(a)        log2f(a)
141   #define PetscSinReal(a)         sinf(a)
142   #define PetscCosReal(a)         cosf(a)
143   #define PetscTanReal(a)         tanf(a)
144   #define PetscAsinReal(a)        asinf(a)
145   #define PetscAcosReal(a)        acosf(a)
146   #define PetscAtanReal(a)        atanf(a)
147   #define PetscSinhReal(a)        sinhf(a)
148   #define PetscCoshReal(a)        coshf(a)
149   #define PetscTanhReal(a)        tanhf(a)
150   #define PetscAsinhReal(a)       asinhf(a)
151   #define PetscAcoshReal(a)       acoshf(a)
152   #define PetscAtanhReal(a)       atanhf(a)
153   #define PetscErfReal(a)         erff(a)
154   #define PetscCeilReal(a)        ceilf(a)
155   #define PetscFloorReal(a)       floorf(a)
156   #define PetscFmodReal(a, b)     fmodf(a, b)
157   #define PetscCopySignReal(a, b) copysignf(a, b)
158   #define PetscTGamma(a)          tgammaf(a)
159   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
160     #define PetscLGamma(a) gammaf(a)
161   #else
162     #define PetscLGamma(a) lgammaf(a)
163   #endif
164 
165 #endif /* PETSC_USE_REAL_* */
166 
167 static inline PetscReal PetscSignReal(PetscReal a)
168 {
169   return (PetscReal)((a < (PetscReal)0) ? -1 : ((a > (PetscReal)0) ? 1 : 0));
170 }
171 
172 #if !defined(PETSC_HAVE_LOG2)
173   #undef PetscLog2Real
174 static inline PetscReal PetscLog2Real(PetscReal a)
175 {
176   return PetscLogReal(a) / PetscLogReal((PetscReal)2);
177 }
178 #endif
179 
180 #if defined(PETSC_HAVE_REAL___FLOAT128) && !defined(PETSC_SKIP_REAL___FLOAT128)
181 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128 PETSC_ATTRIBUTE_MPI_TYPE_TAG(__float128);
182 #endif
183 #if defined(PETSC_HAVE_REAL___FP16) && !defined(PETSC_SKIP_REAL___FP16)
184 PETSC_EXTERN MPI_Datatype MPIU___FP16 PETSC_ATTRIBUTE_MPI_TYPE_TAG(__fp16);
185 #endif
186 
187 /*MC
188    MPIU_REAL - Portable MPI datatype corresponding to `PetscReal` independent of what precision `PetscReal` is in
189 
190    Notes:
191    In MPI calls that require an MPI datatype that matches a `PetscReal` or array of `PetscReal` values, pass this value.
192 
193    Level: beginner
194 
195 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_SCALAR`, `MPIU_COMPLEX`, `MPIU_INT`
196 M*/
197 #if defined(PETSC_USE_REAL_SINGLE)
198   #define MPIU_REAL MPI_FLOAT
199 #elif defined(PETSC_USE_REAL_DOUBLE)
200   #define MPIU_REAL MPI_DOUBLE
201 #elif defined(PETSC_USE_REAL___FLOAT128)
202   #define MPIU_REAL MPIU___FLOAT128
203 #elif defined(PETSC_USE_REAL___FP16)
204   #define MPIU_REAL MPIU___FP16
205 #endif /* PETSC_USE_REAL_* */
206 
207 /*
208     Complex number definitions
209  */
210 #if defined(PETSC_HAVE_COMPLEX)
211   #if defined(__cplusplus) && !defined(PETSC_USE_REAL___FLOAT128)
212   /* C++ support of complex number */
213 
214     #define PetscRealPartComplex(a)      (static_cast<PetscComplex>(a)).real()
215     #define PetscImaginaryPartComplex(a) (static_cast<PetscComplex>(a)).imag()
216     #define PetscAbsComplex(a)           petsccomplexlib::abs(static_cast<PetscComplex>(a))
217     #define PetscArgComplex(a)           petsccomplexlib::arg(static_cast<PetscComplex>(a))
218     #define PetscConjComplex(a)          petsccomplexlib::conj(static_cast<PetscComplex>(a))
219     #define PetscSqrtComplex(a)          petsccomplexlib::sqrt(static_cast<PetscComplex>(a))
220     #define PetscPowComplex(a, b)        petsccomplexlib::pow(static_cast<PetscComplex>(a), static_cast<PetscComplex>(b))
221     #define PetscExpComplex(a)           petsccomplexlib::exp(static_cast<PetscComplex>(a))
222     #define PetscLogComplex(a)           petsccomplexlib::log(static_cast<PetscComplex>(a))
223     #define PetscSinComplex(a)           petsccomplexlib::sin(static_cast<PetscComplex>(a))
224     #define PetscCosComplex(a)           petsccomplexlib::cos(static_cast<PetscComplex>(a))
225     #define PetscTanComplex(a)           petsccomplexlib::tan(static_cast<PetscComplex>(a))
226     #define PetscAsinComplex(a)          petsccomplexlib::asin(static_cast<PetscComplex>(a))
227     #define PetscAcosComplex(a)          petsccomplexlib::acos(static_cast<PetscComplex>(a))
228     #define PetscAtanComplex(a)          petsccomplexlib::atan(static_cast<PetscComplex>(a))
229     #define PetscSinhComplex(a)          petsccomplexlib::sinh(static_cast<PetscComplex>(a))
230     #define PetscCoshComplex(a)          petsccomplexlib::cosh(static_cast<PetscComplex>(a))
231     #define PetscTanhComplex(a)          petsccomplexlib::tanh(static_cast<PetscComplex>(a))
232     #define PetscAsinhComplex(a)         petsccomplexlib::asinh(static_cast<PetscComplex>(a))
233     #define PetscAcoshComplex(a)         petsccomplexlib::acosh(static_cast<PetscComplex>(a))
234     #define PetscAtanhComplex(a)         petsccomplexlib::atanh(static_cast<PetscComplex>(a))
235 
236   /* TODO: Add configure tests
237 
238 #if !defined(PETSC_HAVE_CXX_TAN_COMPLEX)
239 #undef PetscTanComplex
240 static inline PetscComplex PetscTanComplex(PetscComplex z)
241 {
242   return PetscSinComplex(z)/PetscCosComplex(z);
243 }
244 #endif
245 
246 #if !defined(PETSC_HAVE_CXX_TANH_COMPLEX)
247 #undef PetscTanhComplex
248 static inline PetscComplex PetscTanhComplex(PetscComplex z)
249 {
250   return PetscSinhComplex(z)/PetscCoshComplex(z);
251 }
252 #endif
253 
254 #if !defined(PETSC_HAVE_CXX_ASIN_COMPLEX)
255 #undef PetscAsinComplex
256 static inline PetscComplex PetscAsinComplex(PetscComplex z)
257 {
258   const PetscComplex j(0,1);
259   return -j*PetscLogComplex(j*z+PetscSqrtComplex(1.0f-z*z));
260 }
261 #endif
262 
263 #if !defined(PETSC_HAVE_CXX_ACOS_COMPLEX)
264 #undef PetscAcosComplex
265 static inline PetscComplex PetscAcosComplex(PetscComplex z)
266 {
267   const PetscComplex j(0,1);
268   return j*PetscLogComplex(z-j*PetscSqrtComplex(1.0f-z*z));
269 }
270 #endif
271 
272 #if !defined(PETSC_HAVE_CXX_ATAN_COMPLEX)
273 #undef PetscAtanComplex
274 static inline PetscComplex PetscAtanComplex(PetscComplex z)
275 {
276   const PetscComplex j(0,1);
277   return 0.5f*j*PetscLogComplex((1.0f-j*z)/(1.0f+j*z));
278 }
279 #endif
280 
281 #if !defined(PETSC_HAVE_CXX_ASINH_COMPLEX)
282 #undef PetscAsinhComplex
283 static inline PetscComplex PetscAsinhComplex(PetscComplex z)
284 {
285   return PetscLogComplex(z+PetscSqrtComplex(z*z+1.0f));
286 }
287 #endif
288 
289 #if !defined(PETSC_HAVE_CXX_ACOSH_COMPLEX)
290 #undef PetscAcoshComplex
291 static inline PetscComplex PetscAcoshComplex(PetscComplex z)
292 {
293   return PetscLogComplex(z+PetscSqrtComplex(z*z-1.0f));
294 }
295 #endif
296 
297 #if !defined(PETSC_HAVE_CXX_ATANH_COMPLEX)
298 #undef PetscAtanhComplex
299 static inline PetscComplex PetscAtanhComplex(PetscComplex z)
300 {
301   return 0.5f*PetscLogComplex((1.0f+z)/(1.0f-z));
302 }
303 #endif
304 
305 */
306 
307   #else /* C99 support of complex number */
308 
309     #if defined(PETSC_USE_REAL_SINGLE)
310       #define PetscRealPartComplex(a)      crealf(a)
311       #define PetscImaginaryPartComplex(a) cimagf(a)
312       #define PetscAbsComplex(a)           cabsf(a)
313       #define PetscArgComplex(a)           cargf(a)
314       #define PetscConjComplex(a)          conjf(a)
315       #define PetscSqrtComplex(a)          csqrtf(a)
316       #define PetscPowComplex(a, b)        cpowf(a, b)
317       #define PetscExpComplex(a)           cexpf(a)
318       #define PetscLogComplex(a)           clogf(a)
319       #define PetscSinComplex(a)           csinf(a)
320       #define PetscCosComplex(a)           ccosf(a)
321       #define PetscTanComplex(a)           ctanf(a)
322       #define PetscAsinComplex(a)          casinf(a)
323       #define PetscAcosComplex(a)          cacosf(a)
324       #define PetscAtanComplex(a)          catanf(a)
325       #define PetscSinhComplex(a)          csinhf(a)
326       #define PetscCoshComplex(a)          ccoshf(a)
327       #define PetscTanhComplex(a)          ctanhf(a)
328       #define PetscAsinhComplex(a)         casinhf(a)
329       #define PetscAcoshComplex(a)         cacoshf(a)
330       #define PetscAtanhComplex(a)         catanhf(a)
331 
332     #elif defined(PETSC_USE_REAL_DOUBLE)
333       #define PetscRealPartComplex(a)      creal(a)
334       #define PetscImaginaryPartComplex(a) cimag(a)
335       #define PetscAbsComplex(a)           cabs(a)
336       #define PetscArgComplex(a)           carg(a)
337       #define PetscConjComplex(a)          conj(a)
338       #define PetscSqrtComplex(a)          csqrt(a)
339       #define PetscPowComplex(a, b)        cpow(a, b)
340       #define PetscExpComplex(a)           cexp(a)
341       #define PetscLogComplex(a)           clog(a)
342       #define PetscSinComplex(a)           csin(a)
343       #define PetscCosComplex(a)           ccos(a)
344       #define PetscTanComplex(a)           ctan(a)
345       #define PetscAsinComplex(a)          casin(a)
346       #define PetscAcosComplex(a)          cacos(a)
347       #define PetscAtanComplex(a)          catan(a)
348       #define PetscSinhComplex(a)          csinh(a)
349       #define PetscCoshComplex(a)          ccosh(a)
350       #define PetscTanhComplex(a)          ctanh(a)
351       #define PetscAsinhComplex(a)         casinh(a)
352       #define PetscAcoshComplex(a)         cacosh(a)
353       #define PetscAtanhComplex(a)         catanh(a)
354 
355     #elif defined(PETSC_USE_REAL___FLOAT128)
356       #define PetscRealPartComplex(a)      crealq(a)
357       #define PetscImaginaryPartComplex(a) cimagq(a)
358       #define PetscAbsComplex(a)           cabsq(a)
359       #define PetscArgComplex(a)           cargq(a)
360       #define PetscConjComplex(a)          conjq(a)
361       #define PetscSqrtComplex(a)          csqrtq(a)
362       #define PetscPowComplex(a, b)        cpowq(a, b)
363       #define PetscExpComplex(a)           cexpq(a)
364       #define PetscLogComplex(a)           clogq(a)
365       #define PetscSinComplex(a)           csinq(a)
366       #define PetscCosComplex(a)           ccosq(a)
367       #define PetscTanComplex(a)           ctanq(a)
368       #define PetscAsinComplex(a)          casinq(a)
369       #define PetscAcosComplex(a)          cacosq(a)
370       #define PetscAtanComplex(a)          catanq(a)
371       #define PetscSinhComplex(a)          csinhq(a)
372       #define PetscCoshComplex(a)          ccoshq(a)
373       #define PetscTanhComplex(a)          ctanhq(a)
374       #define PetscAsinhComplex(a)         casinhq(a)
375       #define PetscAcoshComplex(a)         cacoshq(a)
376       #define PetscAtanhComplex(a)         catanhq(a)
377 
378     #endif /* PETSC_USE_REAL_* */
379   #endif   /* (__cplusplus) */
380 
381 /*
382    PETSC_i is the imaginary number, i
383 */
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    Notes:
438    In MPI calls that require an MPI datatype that matches a `PetscComplex` or array of `PetscComplex` values, pass this value.
439 
440    Level: beginner
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    Notes:
464    In MPI calls that require an MPI datatype that matches a `PetscScalar` or array of `PetscScalar` values, pass this value.
465 
466    Level: beginner
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    Notes:
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
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
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 .seealso: `PetscMin()`, `PetscMax()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
674 M*/
675 #define PetscClipInterval(x, a, b) (PetscMax((a), PetscMin((x), (b))))
676 
677 /*MC
678    PetscAbsInt - Returns the absolute value of an integer
679 
680    Synopsis:
681    #include <petscmath.h>
682    int abs PetscAbsInt(int v1)
683 
684    Input Parameter:
685 .   v1 - the integer
686 
687    Level: beginner
688 
689 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsReal()`, `PetscSqr()`
690 M*/
691 #define PetscAbsInt(a) (((a) < 0) ? (-(a)) : (a))
692 
693 /*MC
694    PetscAbsReal - Returns the absolute value of an real number
695 
696    Synopsis:
697    #include <petscmath.h>
698    Real abs PetscAbsReal(PetscReal v1)
699 
700    Input Parameter:
701 .   v1 - the `PetscReal` value
702 
703    Level: beginner
704 
705 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscSqr()`
706 M*/
707 #if defined(PETSC_USE_REAL_SINGLE)
708   #define PetscAbsReal(a) fabsf(a)
709 #elif defined(PETSC_USE_REAL_DOUBLE)
710   #define PetscAbsReal(a) fabs(a)
711 #elif defined(PETSC_USE_REAL___FLOAT128)
712   #define PetscAbsReal(a) fabsq(a)
713 #elif defined(PETSC_USE_REAL___FP16)
714   #define PetscAbsReal(a) fabsf(a)
715 #endif
716 
717 /*MC
718    PetscSqr - Returns the square of a number
719 
720    Synopsis:
721    #include <petscmath.h>
722    type sqr PetscSqr(type v1)
723 
724    Not Collective
725 
726    Input Parameter:
727 .   v1 - the value
728 
729    Level: beginner
730 
731    Note:
732    The type can be integer or floating point value
733 
734 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`
735 M*/
736 #define PetscSqr(a) ((a) * (a))
737 
738 #if defined(PETSC_USE_REAL_SINGLE)
739   #define PetscRealConstant(constant) constant##F
740 #elif defined(PETSC_USE_REAL_DOUBLE)
741   #define PetscRealConstant(constant) constant
742 #elif defined(PETSC_USE_REAL___FLOAT128)
743   #define PetscRealConstant(constant) constant##Q
744 #elif defined(PETSC_USE_REAL___FP16)
745   #define PetscRealConstant(constant) constant##F
746 #endif
747 
748 /*
749      Basic constants
750 */
751 #define PETSC_PI    PetscRealConstant(3.1415926535897932384626433832795029)
752 #define PETSC_PHI   PetscRealConstant(1.6180339887498948482045868343656381)
753 #define PETSC_SQRT2 PetscRealConstant(1.4142135623730950488016887242096981)
754 
755 #if defined(PETSC_USE_REAL_SINGLE)
756   #define PETSC_MAX_REAL             3.40282346638528860e+38F
757   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
758   #define PETSC_REAL_MIN             1.1754944e-38F
759   #define PETSC_MACHINE_EPSILON      1.19209290e-07F
760   #define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F
761   #define PETSC_SMALL                1.e-5F
762 #elif defined(PETSC_USE_REAL_DOUBLE)
763   #define PETSC_MAX_REAL             1.7976931348623157e+308
764   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
765   #define PETSC_REAL_MIN             2.225073858507201e-308
766   #define PETSC_MACHINE_EPSILON      2.2204460492503131e-16
767   #define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08
768   #define PETSC_SMALL                1.e-10
769 #elif defined(PETSC_USE_REAL___FLOAT128)
770   #define PETSC_MAX_REAL             FLT128_MAX
771   #define PETSC_MIN_REAL             (-FLT128_MAX)
772   #define PETSC_REAL_MIN             FLT128_MIN
773   #define PETSC_MACHINE_EPSILON      FLT128_EPSILON
774   #define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q
775   #define PETSC_SMALL                1.e-20Q
776 #elif defined(PETSC_USE_REAL___FP16)
777   #define PETSC_MAX_REAL             65504.0F
778   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
779   #define PETSC_REAL_MIN             .00006103515625F
780   #define PETSC_MACHINE_EPSILON      .0009765625F
781   #define PETSC_SQRT_MACHINE_EPSILON .03125F
782   #define PETSC_SMALL                5.e-3F
783 #endif
784 
785 /*MC
786     PETSC_INFINITY - a finite number that represents infinity for setting certain bounds in `Tao`
787 
788    Level: intermediate
789 
790   Note:
791   This is not the IEEE infinity value
792 
793 .seealso: `PETSC_NINFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
794 M*/
795 #define PETSC_INFINITY  (PETSC_MAX_REAL / 4)
796 
797 /*MC
798     PETSC_NINFINITY - a finite number that represents negative infinity for setting certain bounds in `Tao`
799 
800    Level: intermediate
801 
802   Note:
803   This is not the negative IEEE infinity value
804 
805 .seealso: `PETSC_INFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
806 M*/
807 #define PETSC_NINFINITY (-PETSC_INFINITY)
808 
809 PETSC_EXTERN PetscBool  PetscIsInfReal(PetscReal);
810 PETSC_EXTERN PetscBool  PetscIsNanReal(PetscReal);
811 PETSC_EXTERN PetscBool  PetscIsNormalReal(PetscReal);
812 static inline PetscBool PetscIsInfOrNanReal(PetscReal v)
813 {
814   return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;
815 }
816 static inline PetscBool PetscIsInfScalar(PetscScalar v)
817 {
818   return PetscIsInfReal(PetscAbsScalar(v));
819 }
820 static inline PetscBool PetscIsNanScalar(PetscScalar v)
821 {
822   return PetscIsNanReal(PetscAbsScalar(v));
823 }
824 static inline PetscBool PetscIsInfOrNanScalar(PetscScalar v)
825 {
826   return PetscIsInfOrNanReal(PetscAbsScalar(v));
827 }
828 static inline PetscBool PetscIsNormalScalar(PetscScalar v)
829 {
830   return PetscIsNormalReal(PetscAbsScalar(v));
831 }
832 
833 PETSC_EXTERN PetscBool PetscIsCloseAtTol(PetscReal, PetscReal, PetscReal, PetscReal);
834 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal, PetscReal);
835 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar, PetscScalar);
836 
837 /*@C
838   PetscIsCloseAtTolScalar - Like `PetscIsCloseAtTol()` but for `PetscScalar`
839 
840   Input Parameters:
841 + lhs  - The first number
842 . rhs  - The second number
843 . rtol - The relative tolerance
844 - atol - The absolute tolerance
845 
846   Level: beginner
847 
848   Note:
849   This routine is equivalent to `PetscIsCloseAtTol()` when PETSc is configured without complex
850   numbers.
851 
852 .seealso: `PetscIsCloseAtTol()`
853 @*/
854 static inline PetscBool PetscIsCloseAtTolScalar(PetscScalar lhs, PetscScalar rhs, PetscReal rtol, PetscReal atol)
855 {
856   PetscBool close = PetscIsCloseAtTol(PetscRealPart(lhs), PetscRealPart(rhs), rtol, atol);
857 
858   if (PetscDefined(USE_COMPLEX)) close = (PetscBool)(close && PetscIsCloseAtTol(PetscImaginaryPart(lhs), PetscImaginaryPart(rhs), rtol, atol));
859   return close;
860 }
861 
862 /*
863     These macros are currently hardwired to match the regular data types, so there is no support for a different
864     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
865  */
866 #define MPIU_MATSCALAR MPIU_SCALAR
867 typedef PetscScalar MatScalar;
868 typedef PetscReal   MatReal;
869 
870 struct petsc_mpiu_2scalar {
871   PetscScalar a, b;
872 };
873 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2scalar);
874 
875 /* MPI Datatypes for composite reductions */
876 struct petsc_mpiu_real_int {
877   PetscReal v;
878   PetscInt  i;
879 };
880 
881 struct petsc_mpiu_scalar_int {
882   PetscScalar v;
883   PetscInt    i;
884 };
885 
886 PETSC_EXTERN MPI_Datatype MPIU_REAL_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_real_int);
887 PETSC_EXTERN MPI_Datatype MPIU_SCALAR_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_scalar_int);
888 
889 #if defined(PETSC_USE_64BIT_INDICES)
890 struct /* __attribute__((packed, aligned(alignof(PetscInt *)))) */ petsc_mpiu_2int {
891   PetscInt a;
892   PetscInt b;
893 };
894 /*
895  static_assert(sizeof(struct petsc_mpiu_2int) == 2 * sizeof(PetscInt), "");
896  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt *), "");
897  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt[2]), "");
898 
899  clang generates warnings that petsc_mpiu_2int is not layout compatible with PetscInt[2] or
900  PetscInt *, even though (with everything else uncommented) both of the static_asserts above
901  pass! So we just comment it out...
902 */
903 PETSC_EXTERN MPI_Datatype MPIU_2INT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2int) */;
904 #else
905   #define MPIU_2INT MPI_2INT
906 #endif
907 PETSC_EXTERN MPI_Datatype MPI_4INT;
908 PETSC_EXTERN MPI_Datatype MPIU_4INT;
909 
910 static inline PetscInt PetscPowInt(PetscInt base, PetscInt power)
911 {
912   PetscInt result = 1;
913   while (power) {
914     if (power & 1) result *= base;
915     power >>= 1;
916     if (power) base *= base;
917   }
918   return result;
919 }
920 
921 static inline PetscInt64 PetscPowInt64(PetscInt base, PetscInt power)
922 {
923   PetscInt64 result = 1;
924   while (power) {
925     if (power & 1) result *= base;
926     power >>= 1;
927     if (power) base *= base;
928   }
929   return result;
930 }
931 
932 static inline PetscReal PetscPowRealInt(PetscReal base, PetscInt power)
933 {
934   PetscReal result = 1;
935   if (power < 0) {
936     power = -power;
937     base  = ((PetscReal)1) / base;
938   }
939   while (power) {
940     if (power & 1) result *= base;
941     power >>= 1;
942     if (power) base *= base;
943   }
944   return result;
945 }
946 
947 static inline PetscScalar PetscPowScalarInt(PetscScalar base, PetscInt power)
948 {
949   PetscScalar result = (PetscReal)1;
950   if (power < 0) {
951     power = -power;
952     base  = ((PetscReal)1) / base;
953   }
954   while (power) {
955     if (power & 1) result *= base;
956     power >>= 1;
957     if (power) base *= base;
958   }
959   return result;
960 }
961 
962 static inline PetscScalar PetscPowScalarReal(PetscScalar base, PetscReal power)
963 {
964   PetscScalar cpower = power;
965   return PetscPowScalar(base, cpower);
966 }
967 
968 /*MC
969     PetscApproximateLTE - Performs a less than or equal to on a given constant with a fudge for floating point numbers
970 
971    Synopsis:
972    #include <petscmath.h>
973    bool PetscApproximateLTE(PetscReal x,constant float)
974 
975    Not Collective
976 
977    Input Parameters:
978 +   x - the variable
979 -   b - the constant float it is checking if `x` is less than or equal to
980 
981    Level: advanced
982 
983    Notes:
984      The fudge factor is the value `PETSC_SMALL`
985 
986      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
987 
988      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
989      floating point results.
990 
991 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateGTE()`
992 M*/
993 #define PetscApproximateLTE(x, b) ((x) <= (PetscRealConstant(b) + PETSC_SMALL))
994 
995 /*MC
996     PetscApproximateGTE - Performs a greater than or equal to on a given constant with a fudge for floating point numbers
997 
998    Synopsis:
999    #include <petscmath.h>
1000    bool PetscApproximateGTE(PetscReal x,constant float)
1001 
1002    Not Collective
1003 
1004    Input Parameters:
1005 +   x - the variable
1006 -   b - the constant float it is checking if `x` is greater than or equal to
1007 
1008    Level: advanced
1009 
1010    Notes:
1011      The fudge factor is the value `PETSC_SMALL`
1012 
1013      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
1014 
1015      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1016      floating point results.
1017 
1018 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1019 M*/
1020 #define PetscApproximateGTE(x, b) ((x) >= (PetscRealConstant(b) - PETSC_SMALL))
1021 
1022 /*MC
1023     PetscCeilInt - Returns the ceiling of the quotation of two positive integers
1024 
1025    Synopsis:
1026    #include <petscmath.h>
1027    PetscInt PetscCeilInt(PetscInt x,PetscInt y)
1028 
1029    Not Collective
1030 
1031    Input Parameters:
1032 +   x - the numerator
1033 -   y - the denominator
1034 
1035    Level: advanced
1036 
1037 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1038 M*/
1039 #define PetscCeilInt(x, y) ((((PetscInt)(x)) / ((PetscInt)(y))) + ((((PetscInt)(x)) % ((PetscInt)(y))) ? 1 : 0))
1040 
1041 #define PetscCeilInt64(x, y) ((((PetscInt64)(x)) / ((PetscInt64)(y))) + ((((PetscInt64)(x)) % ((PetscInt64)(y))) ? 1 : 0))
1042 
1043 PETSC_EXTERN PetscErrorCode PetscLinearRegression(PetscInt, const PetscReal[], const PetscReal[], PetscReal *, PetscReal *);
1044