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