xref: /petsc/include/petscmath.h (revision 6dd63270497ad23dcf16ae500a87ff2b2a0b7474)
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 PetscRintReal(a)        rintf(a)
51   #define PetscFmodReal(a, b)     fmodf(a, b)
52   #define PetscCopysignReal(a, b) copysignf(a, b)
53   #define PetscTGamma(a)          tgammaf(a)
54   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
55     #define PetscLGamma(a) gammaf(a)
56   #else
57     #define PetscLGamma(a) lgammaf(a)
58   #endif
59 
60 #elif defined(PETSC_USE_REAL_DOUBLE)
61   #define PetscSqrtReal(a)        sqrt(a)
62   #define PetscCbrtReal(a)        cbrt(a)
63   #define PetscHypotReal(a, b)    hypot(a, b)
64   #define PetscAtan2Real(a, b)    atan2(a, b)
65   #define PetscPowReal(a, b)      pow(a, b)
66   #define PetscExpReal(a)         exp(a)
67   #define PetscLogReal(a)         log(a)
68   #define PetscLog10Real(a)       log10(a)
69   #define PetscLog2Real(a)        log2(a)
70   #define PetscSinReal(a)         sin(a)
71   #define PetscCosReal(a)         cos(a)
72   #define PetscTanReal(a)         tan(a)
73   #define PetscAsinReal(a)        asin(a)
74   #define PetscAcosReal(a)        acos(a)
75   #define PetscAtanReal(a)        atan(a)
76   #define PetscSinhReal(a)        sinh(a)
77   #define PetscCoshReal(a)        cosh(a)
78   #define PetscTanhReal(a)        tanh(a)
79   #define PetscAsinhReal(a)       asinh(a)
80   #define PetscAcoshReal(a)       acosh(a)
81   #define PetscAtanhReal(a)       atanh(a)
82   #define PetscErfReal(a)         erf(a)
83   #define PetscCeilReal(a)        ceil(a)
84   #define PetscFloorReal(a)       floor(a)
85   #define PetscRintReal(a)        rint(a)
86   #define PetscFmodReal(a, b)     fmod(a, b)
87   #define PetscCopysignReal(a, b) copysign(a, b)
88   #define PetscTGamma(a)          tgamma(a)
89   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
90     #define PetscLGamma(a) gamma(a)
91   #else
92     #define PetscLGamma(a) lgamma(a)
93   #endif
94 
95 #elif defined(PETSC_USE_REAL___FLOAT128)
96   #define PetscSqrtReal(a)        sqrtq(a)
97   #define PetscCbrtReal(a)        cbrtq(a)
98   #define PetscHypotReal(a, b)    hypotq(a, b)
99   #define PetscAtan2Real(a, b)    atan2q(a, b)
100   #define PetscPowReal(a, b)      powq(a, b)
101   #define PetscExpReal(a)         expq(a)
102   #define PetscLogReal(a)         logq(a)
103   #define PetscLog10Real(a)       log10q(a)
104   #define PetscLog2Real(a)        log2q(a)
105   #define PetscSinReal(a)         sinq(a)
106   #define PetscCosReal(a)         cosq(a)
107   #define PetscTanReal(a)         tanq(a)
108   #define PetscAsinReal(a)        asinq(a)
109   #define PetscAcosReal(a)        acosq(a)
110   #define PetscAtanReal(a)        atanq(a)
111   #define PetscSinhReal(a)        sinhq(a)
112   #define PetscCoshReal(a)        coshq(a)
113   #define PetscTanhReal(a)        tanhq(a)
114   #define PetscAsinhReal(a)       asinhq(a)
115   #define PetscAcoshReal(a)       acoshq(a)
116   #define PetscAtanhReal(a)       atanhq(a)
117   #define PetscErfReal(a)         erfq(a)
118   #define PetscCeilReal(a)        ceilq(a)
119   #define PetscFloorReal(a)       floorq(a)
120   #define PetscRintReal(a)        rintq(a)
121   #define PetscFmodReal(a, b)     fmodq(a, b)
122   #define PetscCopysignReal(a, b) copysignq(a, b)
123   #define PetscTGamma(a)          tgammaq(a)
124   #if defined(PETSC_HAVE_LGAMMA_IS_GAMMA)
125     #define PetscLGamma(a) gammaq(a)
126   #else
127     #define PetscLGamma(a) lgammaq(a)
128   #endif
129 
130 #elif defined(PETSC_USE_REAL___FP16)
131   #define PetscSqrtReal(a)        sqrtf(a)
132   #define PetscCbrtReal(a)        cbrtf(a)
133   #define PetscHypotReal(a, b)    hypotf(a, b)
134   #define PetscAtan2Real(a, b)    atan2f(a, b)
135   #define PetscPowReal(a, b)      powf(a, b)
136   #define PetscExpReal(a)         expf(a)
137   #define PetscLogReal(a)         logf(a)
138   #define PetscLog10Real(a)       log10f(a)
139   #define PetscLog2Real(a)        log2f(a)
140   #define PetscSinReal(a)         sinf(a)
141   #define PetscCosReal(a)         cosf(a)
142   #define PetscTanReal(a)         tanf(a)
143   #define PetscAsinReal(a)        asinf(a)
144   #define PetscAcosReal(a)        acosf(a)
145   #define PetscAtanReal(a)        atanf(a)
146   #define PetscSinhReal(a)        sinhf(a)
147   #define PetscCoshReal(a)        coshf(a)
148   #define PetscTanhReal(a)        tanhf(a)
149   #define PetscAsinhReal(a)       asinhf(a)
150   #define PetscAcoshReal(a)       acoshf(a)
151   #define PetscAtanhReal(a)       atanhf(a)
152   #define PetscErfReal(a)         erff(a)
153   #define PetscCeilReal(a)        ceilf(a)
154   #define PetscFloorReal(a)       floorf(a)
155   #define PetscRintReal(a)        rintf(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    Level: beginner
191 
192    Note:
193    In MPI calls that require an MPI datatype that matches a `PetscReal` or array of `PetscReal` values, pass this value.
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 /*MC
382    PETSC_i - the pure imaginary complex number i
383 
384    Level: intermediate
385 
386 .seealso: `PetscComplex`, `PetscScalar`
387 M*/
388 PETSC_EXTERN PetscComplex PETSC_i;
389 
390 /*
391    Try to do the right thing for complex number construction: see
392    http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1464.htm
393    for details
394 */
395 static inline PetscComplex PetscCMPLX(PetscReal x, PetscReal y)
396 {
397   #if defined(__cplusplus) && !defined(PETSC_USE_REAL___FLOAT128)
398   return PetscComplex(x, y);
399   #elif defined(_Imaginary_I)
400   return x + y * _Imaginary_I;
401   #else
402   { /* In both C99 and C11 (ISO/IEC 9899, Section 6.2.5),
403 
404        "For each floating type there is a corresponding real type, which is always a real floating
405        type. For real floating types, it is the same type. For complex types, it is the type given
406        by deleting the keyword _Complex from the type name."
407 
408        So type punning should be portable. */
409     union
410     {
411       PetscComplex z;
412       PetscReal    f[2];
413     } uz;
414 
415     uz.f[0] = x;
416     uz.f[1] = y;
417     return uz.z;
418   }
419   #endif
420 }
421 
422   #define MPIU_C_COMPLEX        MPI_C_COMPLEX PETSC_DEPRECATED_MACRO(3, 15, 0, "MPI_C_COMPLEX", )
423   #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX PETSC_DEPRECATED_MACRO(3, 15, 0, "MPI_C_DOUBLE_COMPLEX", )
424 
425   #if defined(PETSC_HAVE_REAL___FLOAT128) && !defined(PETSC_SKIP_REAL___FLOAT128)
426     // if complex is not used, then quadmath.h won't be included by petscsystypes.h
427     #if defined(PETSC_USE_COMPLEX)
428       #define MPIU___COMPLEX128_ATTR_TAG PETSC_ATTRIBUTE_MPI_TYPE_TAG(__complex128)
429     #else
430       #define MPIU___COMPLEX128_ATTR_TAG
431     #endif
432 
433 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 MPIU___COMPLEX128_ATTR_TAG;
434 
435     #undef MPIU___COMPLEX128_ATTR_TAG
436   #endif /* PETSC_HAVE_REAL___FLOAT128 */
437 
438   /*MC
439    MPIU_COMPLEX - Portable MPI datatype corresponding to `PetscComplex` independent of the precision of `PetscComplex`
440 
441    Level: beginner
442 
443    Note:
444    In MPI calls that require an MPI datatype that matches a `PetscComplex` or array of `PetscComplex` values, pass this value.
445 
446 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_REAL`, `MPIU_SCALAR`, `MPIU_COMPLEX`, `MPIU_INT`, `PETSC_i`
447 M*/
448   #if defined(PETSC_USE_REAL_SINGLE)
449     #define MPIU_COMPLEX MPI_C_COMPLEX
450   #elif defined(PETSC_USE_REAL_DOUBLE)
451     #define MPIU_COMPLEX MPI_C_DOUBLE_COMPLEX
452   #elif defined(PETSC_USE_REAL___FLOAT128)
453     #define MPIU_COMPLEX MPIU___COMPLEX128
454   #elif defined(PETSC_USE_REAL___FP16)
455     #define MPIU_COMPLEX MPI_C_COMPLEX
456   #endif /* PETSC_USE_REAL_* */
457 
458 #endif /* PETSC_HAVE_COMPLEX */
459 
460 /*
461     Scalar number definitions
462  */
463 #if defined(PETSC_USE_COMPLEX) && defined(PETSC_HAVE_COMPLEX)
464   /*MC
465    MPIU_SCALAR - Portable MPI datatype corresponding to `PetscScalar` independent of the precision of `PetscScalar`
466 
467    Level: beginner
468 
469    Note:
470    In MPI calls that require an MPI datatype that matches a `PetscScalar` or array of `PetscScalar` values, pass this value.
471 
472 .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_REAL`, `MPIU_COMPLEX`, `MPIU_INT`
473 M*/
474   #define MPIU_SCALAR MPIU_COMPLEX
475 
476   /*MC
477    PetscRealPart - Returns the real part of a `PetscScalar`
478 
479    Synopsis:
480    #include <petscmath.h>
481    PetscReal PetscRealPart(PetscScalar v)
482 
483    Not Collective
484 
485    Input Parameter:
486 .  v - value to find the real part of
487 
488    Level: beginner
489 
490 .seealso: `PetscScalar`, `PetscImaginaryPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
491 M*/
492   #define PetscRealPart(a) PetscRealPartComplex(a)
493 
494   /*MC
495    PetscImaginaryPart - Returns the imaginary part of a `PetscScalar`
496 
497    Synopsis:
498    #include <petscmath.h>
499    PetscReal PetscImaginaryPart(PetscScalar v)
500 
501    Not Collective
502 
503    Input Parameter:
504 .  v - value to find the imaginary part of
505 
506    Level: beginner
507 
508    Note:
509    If PETSc was configured for real numbers then this always returns the value 0
510 
511 .seealso: `PetscScalar`, `PetscRealPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
512 M*/
513   #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a)
514 
515   #define PetscAbsScalar(a)    PetscAbsComplex(a)
516   #define PetscArgScalar(a)    PetscArgComplex(a)
517   #define PetscConj(a)         PetscConjComplex(a)
518   #define PetscSqrtScalar(a)   PetscSqrtComplex(a)
519   #define PetscPowScalar(a, b) PetscPowComplex(a, b)
520   #define PetscExpScalar(a)    PetscExpComplex(a)
521   #define PetscLogScalar(a)    PetscLogComplex(a)
522   #define PetscSinScalar(a)    PetscSinComplex(a)
523   #define PetscCosScalar(a)    PetscCosComplex(a)
524   #define PetscTanScalar(a)    PetscTanComplex(a)
525   #define PetscAsinScalar(a)   PetscAsinComplex(a)
526   #define PetscAcosScalar(a)   PetscAcosComplex(a)
527   #define PetscAtanScalar(a)   PetscAtanComplex(a)
528   #define PetscSinhScalar(a)   PetscSinhComplex(a)
529   #define PetscCoshScalar(a)   PetscCoshComplex(a)
530   #define PetscTanhScalar(a)   PetscTanhComplex(a)
531   #define PetscAsinhScalar(a)  PetscAsinhComplex(a)
532   #define PetscAcoshScalar(a)  PetscAcoshComplex(a)
533   #define PetscAtanhScalar(a)  PetscAtanhComplex(a)
534 
535 #else /* PETSC_USE_COMPLEX */
536   #define MPIU_SCALAR           MPIU_REAL
537   #define PetscRealPart(a)      (a)
538   #define PetscImaginaryPart(a) ((PetscReal)0)
539   #define PetscAbsScalar(a)     PetscAbsReal(a)
540   #define PetscArgScalar(a)     (((a) < (PetscReal)0) ? PETSC_PI : (PetscReal)0)
541   #define PetscConj(a)          (a)
542   #define PetscSqrtScalar(a)    PetscSqrtReal(a)
543   #define PetscPowScalar(a, b)  PetscPowReal(a, b)
544   #define PetscExpScalar(a)     PetscExpReal(a)
545   #define PetscLogScalar(a)     PetscLogReal(a)
546   #define PetscSinScalar(a)     PetscSinReal(a)
547   #define PetscCosScalar(a)     PetscCosReal(a)
548   #define PetscTanScalar(a)     PetscTanReal(a)
549   #define PetscAsinScalar(a)    PetscAsinReal(a)
550   #define PetscAcosScalar(a)    PetscAcosReal(a)
551   #define PetscAtanScalar(a)    PetscAtanReal(a)
552   #define PetscSinhScalar(a)    PetscSinhReal(a)
553   #define PetscCoshScalar(a)    PetscCoshReal(a)
554   #define PetscTanhScalar(a)    PetscTanhReal(a)
555   #define PetscAsinhScalar(a)   PetscAsinhReal(a)
556   #define PetscAcoshScalar(a)   PetscAcoshReal(a)
557   #define PetscAtanhScalar(a)   PetscAtanhReal(a)
558 
559 #endif /* PETSC_USE_COMPLEX */
560 
561 /*
562    Certain objects may be created using either single or double precision.
563    This is currently not used.
564 */
565 typedef enum {
566   PETSC_SCALAR_DOUBLE,
567   PETSC_SCALAR_SINGLE,
568   PETSC_SCALAR_LONG_DOUBLE,
569   PETSC_SCALAR_HALF
570 } PetscScalarPrecision;
571 
572 /*MC
573    PetscAbs - Returns the absolute value of a number
574 
575    Synopsis:
576    #include <petscmath.h>
577    type PetscAbs(type v)
578 
579    Not Collective
580 
581    Input Parameter:
582 .  v - the number
583 
584    Level: beginner
585 
586    Note:
587    The type can be integer or real floating point value, but cannot be complex
588 
589 .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`, `PetscSign()`
590 M*/
591 #define PetscAbs(a) (((a) >= 0) ? (a) : (-(a)))
592 
593 /*MC
594    PetscSign - Returns the sign of a number as an integer of value -1, 0, or 1
595 
596    Synopsis:
597    #include <petscmath.h>
598    int PetscSign(type v)
599 
600    Not Collective
601 
602    Input Parameter:
603 .  v - the number
604 
605    Level: beginner
606 
607    Note:
608    The type can be integer or real floating point value
609 
610 .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`
611 M*/
612 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1)
613 
614 /*MC
615    PetscMin - Returns minimum of two numbers
616 
617    Synopsis:
618    #include <petscmath.h>
619    type PetscMin(type v1,type v2)
620 
621    Not Collective
622 
623    Input Parameters:
624 +  v1 - first value to find minimum of
625 -  v2 - second value to find minimum of
626 
627    Level: beginner
628 
629    Note:
630    The type can be integer or floating point value, but cannot be complex
631 
632 .seealso: `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
633 M*/
634 #define PetscMin(a, b) (((a) < (b)) ? (a) : (b))
635 
636 /*MC
637    PetscMax - Returns maximum of two numbers
638 
639    Synopsis:
640    #include <petscmath.h>
641    type max PetscMax(type v1,type v2)
642 
643    Not Collective
644 
645    Input Parameters:
646 +  v1 - first value to find maximum of
647 -  v2 - second value to find maximum of
648 
649    Level: beginner
650 
651    Note:
652    The type can be integer or floating point value
653 
654 .seealso: `PetscMin()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
655 M*/
656 #define PetscMax(a, b) (((a) < (b)) ? (b) : (a))
657 
658 /*MC
659    PetscClipInterval - Returns a number clipped to be within an interval
660 
661    Synopsis:
662    #include <petscmath.h>
663    type clip PetscClipInterval(type x,type a,type b)
664 
665    Not Collective
666 
667    Input Parameters:
668 +  x - value to use if within interval [a,b]
669 .  a - lower end of interval
670 -  b - upper end of interval
671 
672    Level: beginner
673 
674    Note:
675    The type can be integer or floating point value
676 
677    Example\:
678 .vb
679   PetscInt c = PetscClipInterval(5, 2, 3); // the value of c is 3
680   PetscInt c = PetscClipInterval(5, 2, 6); // the value of c is 5
681 .ve
682 
683 .seealso: `PetscMin()`, `PetscMax()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
684 M*/
685 #define PetscClipInterval(x, a, b) (PetscMax((a), PetscMin((x), (b))))
686 
687 /*MC
688    PetscAbsInt - Returns the absolute value of an integer
689 
690    Synopsis:
691    #include <petscmath.h>
692    int abs PetscAbsInt(int v1)
693 
694    Input Parameter:
695 .   v1 - the integer
696 
697    Level: beginner
698 
699 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsReal()`, `PetscSqr()`
700 M*/
701 #define PetscAbsInt(a) (((a) < 0) ? (-(a)) : (a))
702 
703 /*MC
704    PetscAbsReal - Returns the absolute value of a real number
705 
706    Synopsis:
707    #include <petscmath.h>
708    Real abs PetscAbsReal(PetscReal v1)
709 
710    Input Parameter:
711 .   v1 - the `PetscReal` value
712 
713    Level: beginner
714 
715 .seealso: `PetscReal`, `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscSqr()`
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    Level: beginner
740 
741    Note:
742    The type can be integer, floating point, or complex floating point
743 
744 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`
745 M*/
746 #define PetscSqr(a) ((a) * (a))
747 
748 /*MC
749    PetscRealConstant - a compile time macro that ensures a given constant real number is properly represented in the configured
750    precision of `PetscReal` be it half, single, double or 128-bit representation
751 
752    Synopsis:
753    #include <petscmath.h>
754    PetscReal PetscRealConstant(real_number)
755 
756    Not Collective
757 
758    Input Parameter:
759 .   v1 - the real number, for example 1.5
760 
761    Level: beginner
762 
763    Note:
764    For example, if PETSc is configured with `--with-precision=__float128` and one writes
765 .vb
766    PetscReal d = 1.5;
767 .ve
768    the result is 1.5 in double precision extended to 128 bit representation, meaning it is very far from the correct value. Hence, one should write
769 .vb
770    PetscReal d = PetscRealConstant(1.5);
771 .ve
772 
773 .seealso: `PetscReal`
774 M*/
775 #if defined(PETSC_USE_REAL_SINGLE)
776   #define PetscRealConstant(constant) constant##F
777 #elif defined(PETSC_USE_REAL_DOUBLE)
778   #define PetscRealConstant(constant) constant
779 #elif defined(PETSC_USE_REAL___FLOAT128)
780   #define PetscRealConstant(constant) constant##Q
781 #elif defined(PETSC_USE_REAL___FP16)
782   #define PetscRealConstant(constant) constant##F
783 #endif
784 
785 /*
786      Basic constants
787 */
788 /*MC
789   PETSC_PI - the value of $ \pi$ to the correct precision of `PetscReal`.
790 
791   Level: beginner
792 
793 .seealso: `PetscReal`, `PETSC_PHI`, `PETSC_SQRT2`
794 M*/
795 
796 /*MC
797   PETSC_PHI - the value of $ \phi$, the Golden Ratio, to the correct precision of `PetscReal`.
798 
799   Level: beginner
800 
801 .seealso: `PetscReal`, `PETSC_PI`, `PETSC_SQRT2`
802 M*/
803 
804 /*MC
805   PETSC_SQRT2 - the value of $ \sqrt{2} $ to the correct precision of `PetscReal`.
806 
807   Level: beginner
808 
809 .seealso: `PetscReal`, `PETSC_PI`, `PETSC_PHI`
810 M*/
811 
812 #define PETSC_PI    PetscRealConstant(3.1415926535897932384626433832795029)
813 #define PETSC_PHI   PetscRealConstant(1.6180339887498948482045868343656381)
814 #define PETSC_SQRT2 PetscRealConstant(1.4142135623730950488016887242096981)
815 
816 /*MC
817   PETSC_MAX_REAL - the largest real value that can be stored in a `PetscReal`
818 
819   Level: beginner
820 
821 .seealso: `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
822 M*/
823 
824 /*MC
825   PETSC_MIN_REAL - the smallest real value that can be stored in a `PetscReal`, generally this is - `PETSC_MAX_REAL`
826 
827   Level: beginner
828 
829 .seealso `PETSC_MAX_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
830 M*/
831 
832 /*MC
833   PETSC_REAL_MIN - the smallest positive normalized real value that can be stored in a `PetscReal`.
834 
835   Level: beginner
836 
837   Note:
838   See <https://en.wikipedia.org/wiki/Subnormal_number> for a discussion of normalized and subnormal floating point numbers
839 
840   Developer Note:
841   The naming is confusing as there is both a `PETSC_REAL_MIN` and `PETSC_MIN_REAL` with different meanings.
842 
843 .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
844 M*/
845 
846 /*MC
847   PETSC_MACHINE_EPSILON - the machine epsilon for the precision of `PetscReal`
848 
849   Level: beginner
850 
851   Note:
852   See <https://en.wikipedia.org/wiki/Machine_epsilon>
853 
854 .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
855 M*/
856 
857 /*MC
858   PETSC_SQRT_MACHINE_EPSILON - the square root of the machine epsilon for the precision of `PetscReal`
859 
860   Level: beginner
861 
862   Note:
863   See `PETSC_MACHINE_EPSILON`
864 
865 .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SMALL`
866 M*/
867 
868 /*MC
869   PETSC_SMALL - an arbitrary "small" number which depends on the precision of `PetscReal` used in some PETSc examples
870   and in `PetscApproximateLTE()` and `PetscApproximateGTE()` to determine if a computation was successful.
871 
872   Level: beginner
873 
874   Note:
875   See `PETSC_MACHINE_EPSILON`
876 
877 .seealso `PetscApproximateLTE()`, `PetscApproximateGTE()`, `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`,
878          `PETSC_SQRT_MACHINE_EPSILON`
879 M*/
880 
881 #if defined(PETSC_USE_REAL_SINGLE)
882   #define PETSC_MAX_REAL             3.40282346638528860e+38F
883   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
884   #define PETSC_REAL_MIN             1.1754944e-38F
885   #define PETSC_MACHINE_EPSILON      1.19209290e-07F
886   #define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F
887   #define PETSC_SMALL                1.e-5F
888 #elif defined(PETSC_USE_REAL_DOUBLE)
889   #define PETSC_MAX_REAL             1.7976931348623157e+308
890   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
891   #define PETSC_REAL_MIN             2.225073858507201e-308
892   #define PETSC_MACHINE_EPSILON      2.2204460492503131e-16
893   #define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08
894   #define PETSC_SMALL                1.e-10
895 #elif defined(PETSC_USE_REAL___FLOAT128)
896   #define PETSC_MAX_REAL             FLT128_MAX
897   #define PETSC_MIN_REAL             (-FLT128_MAX)
898   #define PETSC_REAL_MIN             FLT128_MIN
899   #define PETSC_MACHINE_EPSILON      FLT128_EPSILON
900   #define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q
901   #define PETSC_SMALL                1.e-20Q
902 #elif defined(PETSC_USE_REAL___FP16)
903   #define PETSC_MAX_REAL             65504.0F
904   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
905   #define PETSC_REAL_MIN             .00006103515625F
906   #define PETSC_MACHINE_EPSILON      .0009765625F
907   #define PETSC_SQRT_MACHINE_EPSILON .03125F
908   #define PETSC_SMALL                5.e-3F
909 #endif
910 
911 /*MC
912   PETSC_INFINITY - a finite number that represents infinity for setting certain bounds in `Tao`
913 
914   Level: intermediate
915 
916   Note:
917   This is not the IEEE infinity value
918 
919 .seealso: `PETSC_NINFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
920 M*/
921 #define PETSC_INFINITY (PETSC_MAX_REAL / 4)
922 
923 /*MC
924   PETSC_NINFINITY - a finite number that represents negative infinity for setting certain bounds in `Tao`
925 
926   Level: intermediate
927 
928   Note:
929   This is not the negative IEEE infinity value
930 
931 .seealso: `PETSC_INFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
932 M*/
933 #define PETSC_NINFINITY (-PETSC_INFINITY)
934 
935 PETSC_EXTERN PetscBool  PetscIsInfReal(PetscReal);
936 PETSC_EXTERN PetscBool  PetscIsNanReal(PetscReal);
937 PETSC_EXTERN PetscBool  PetscIsNormalReal(PetscReal);
938 static inline PetscBool PetscIsInfOrNanReal(PetscReal v)
939 {
940   return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;
941 }
942 static inline PetscBool PetscIsInfScalar(PetscScalar v)
943 {
944   return PetscIsInfReal(PetscAbsScalar(v));
945 }
946 static inline PetscBool PetscIsNanScalar(PetscScalar v)
947 {
948   return PetscIsNanReal(PetscAbsScalar(v));
949 }
950 static inline PetscBool PetscIsInfOrNanScalar(PetscScalar v)
951 {
952   return PetscIsInfOrNanReal(PetscAbsScalar(v));
953 }
954 static inline PetscBool PetscIsNormalScalar(PetscScalar v)
955 {
956   return PetscIsNormalReal(PetscAbsScalar(v));
957 }
958 
959 PETSC_EXTERN PetscBool PetscIsCloseAtTol(PetscReal, PetscReal, PetscReal, PetscReal);
960 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal, PetscReal);
961 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar, PetscScalar);
962 
963 /*@C
964   PetscIsCloseAtTolScalar - Like `PetscIsCloseAtTol()` but for `PetscScalar`
965 
966   Input Parameters:
967 + lhs  - The first number
968 . rhs  - The second number
969 . rtol - The relative tolerance
970 - atol - The absolute tolerance
971 
972   Level: beginner
973 
974   Note:
975   This routine is equivalent to `PetscIsCloseAtTol()` when PETSc is configured without complex
976   numbers.
977 
978 .seealso: `PetscIsCloseAtTol()`
979 @*/
980 static inline PetscBool PetscIsCloseAtTolScalar(PetscScalar lhs, PetscScalar rhs, PetscReal rtol, PetscReal atol)
981 {
982   PetscBool close = PetscIsCloseAtTol(PetscRealPart(lhs), PetscRealPart(rhs), rtol, atol);
983 
984   if (PetscDefined(USE_COMPLEX)) close = (PetscBool)(close && PetscIsCloseAtTol(PetscImaginaryPart(lhs), PetscImaginaryPart(rhs), rtol, atol));
985   return close;
986 }
987 
988 /*
989     These macros are currently hardwired to match the regular data types, so there is no support for a different
990     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
991  */
992 #define MPIU_MATSCALAR MPIU_SCALAR
993 typedef PetscScalar MatScalar;
994 typedef PetscReal   MatReal;
995 
996 struct petsc_mpiu_2scalar {
997   PetscScalar a, b;
998 };
999 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2scalar);
1000 
1001 /* MPI Datatypes for composite reductions */
1002 struct petsc_mpiu_real_int {
1003   PetscReal v;
1004   PetscInt  i;
1005 };
1006 
1007 struct petsc_mpiu_scalar_int {
1008   PetscScalar v;
1009   PetscInt    i;
1010 };
1011 
1012 PETSC_EXTERN MPI_Datatype MPIU_REAL_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_real_int);
1013 PETSC_EXTERN MPI_Datatype MPIU_SCALAR_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_scalar_int);
1014 
1015 #if defined(PETSC_USE_64BIT_INDICES)
1016 struct /* __attribute__((packed, aligned(alignof(PetscInt *)))) */ petsc_mpiu_2int {
1017   PetscInt a;
1018   PetscInt b;
1019 };
1020 struct __attribute__((packed)) petsc_mpiu_int_mpiint {
1021   PetscInt    a;
1022   PetscMPIInt b;
1023 };
1024 /*
1025  static_assert(sizeof(struct petsc_mpiu_2int) == 2 * sizeof(PetscInt), "");
1026  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt *), "");
1027  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt[2]), "");
1028 
1029  clang generates warnings that petsc_mpiu_2int is not layout compatible with PetscInt[2] or
1030  PetscInt *, even though (with everything else uncommented) both of the static_asserts above
1031  pass! So we just comment it out...
1032 */
1033 PETSC_EXTERN MPI_Datatype MPIU_2INT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2int) */;
1034 PETSC_EXTERN MPI_Datatype MPIU_INT_MPIINT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_int_mpiint) */;
1035 #else
1036   #define MPIU_2INT       MPI_2INT
1037   #define MPIU_INT_MPIINT MPI_2INT
1038 #endif
1039 PETSC_EXTERN MPI_Datatype MPI_4INT;
1040 PETSC_EXTERN MPI_Datatype MPIU_4INT;
1041 
1042 static inline PetscInt PetscPowInt(PetscInt base, PetscInt power)
1043 {
1044   PetscInt result = 1;
1045   while (power) {
1046     if (power & 1) result *= base;
1047     power >>= 1;
1048     if (power) base *= base;
1049   }
1050   return result;
1051 }
1052 
1053 static inline PetscInt64 PetscPowInt64(PetscInt base, PetscInt power)
1054 {
1055   PetscInt64 result = 1;
1056   while (power) {
1057     if (power & 1) result *= base;
1058     power >>= 1;
1059     if (power) base *= base;
1060   }
1061   return result;
1062 }
1063 
1064 static inline PetscReal PetscPowRealInt(PetscReal base, PetscInt power)
1065 {
1066   PetscReal result = 1;
1067   if (power < 0) {
1068     power = -power;
1069     base  = ((PetscReal)1) / base;
1070   }
1071   while (power) {
1072     if (power & 1) result *= base;
1073     power >>= 1;
1074     if (power) base *= base;
1075   }
1076   return result;
1077 }
1078 
1079 static inline PetscScalar PetscPowScalarInt(PetscScalar base, PetscInt power)
1080 {
1081   PetscScalar result = (PetscReal)1;
1082   if (power < 0) {
1083     power = -power;
1084     base  = ((PetscReal)1) / base;
1085   }
1086   while (power) {
1087     if (power & 1) result *= base;
1088     power >>= 1;
1089     if (power) base *= base;
1090   }
1091   return result;
1092 }
1093 
1094 static inline PetscScalar PetscPowScalarReal(PetscScalar base, PetscReal power)
1095 {
1096   PetscScalar cpower = power;
1097   return PetscPowScalar(base, cpower);
1098 }
1099 
1100 /*MC
1101    PetscApproximateLTE - Performs a less than or equal to on a given constant with a fudge for floating point numbers
1102 
1103    Synopsis:
1104    #include <petscmath.h>
1105    bool PetscApproximateLTE(PetscReal x,constant float)
1106 
1107    Not Collective
1108 
1109    Input Parameters:
1110 +   x - the variable
1111 -   b - the constant float it is checking if `x` is less than or equal to
1112 
1113    Level: advanced
1114 
1115    Notes:
1116    The fudge factor is the value `PETSC_SMALL`
1117 
1118    The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
1119 
1120    This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1121    floating point results.
1122 
1123    Example\:
1124 .vb
1125   PetscReal x;
1126   if (PetscApproximateLTE(x, 3.2)) { // replaces if (x <= 3.2) {
1127 .ve
1128 
1129 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateGTE()`
1130 M*/
1131 #define PetscApproximateLTE(x, b) ((x) <= (PetscRealConstant(b) + PETSC_SMALL))
1132 
1133 /*MC
1134    PetscApproximateGTE - Performs a greater than or equal to on a given constant with a fudge for floating point numbers
1135 
1136    Synopsis:
1137    #include <petscmath.h>
1138    bool PetscApproximateGTE(PetscReal x,constant float)
1139 
1140    Not Collective
1141 
1142    Input Parameters:
1143 +   x - the variable
1144 -   b - the constant float it is checking if `x` is greater than or equal to
1145 
1146    Level: advanced
1147 
1148    Notes:
1149    The fudge factor is the value `PETSC_SMALL`
1150 
1151    The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
1152 
1153    This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1154    floating point results.
1155 
1156    Example\:
1157 .vb
1158   PetscReal x;
1159   if (PetscApproximateGTE(x, 3.2)) {  // replaces if (x >= 3.2) {
1160 .ve
1161 
1162 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1163 M*/
1164 #define PetscApproximateGTE(x, b) ((x) >= (PetscRealConstant(b) - PETSC_SMALL))
1165 
1166 /*@C
1167    PetscCeilInt - Returns the ceiling of the quotation of two positive integers
1168 
1169    Not Collective
1170 
1171    Input Parameters:
1172 +   x - the numerator
1173 -   y - the denominator
1174 
1175    Level: advanced
1176 
1177   Example\:
1178 .vb
1179   PetscInt n = PetscCeilInt(10, 3); // n has the value of 4
1180 .ve
1181 
1182 .seealso: `PetscCeilInt64()`, `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1183 @*/
1184 static inline PetscInt PetscCeilInt(PetscInt x, PetscInt y)
1185 {
1186   return x / y + (x % y ? 1 : 0);
1187 }
1188 
1189 /*@C
1190    PetscCeilInt64 - Returns the ceiling of the quotation of two positive integers
1191 
1192    Not Collective
1193 
1194    Input Parameters:
1195 +   x - the numerator
1196 -   y - the denominator
1197 
1198    Level: advanced
1199 
1200   Example\:
1201 .vb
1202   PetscInt64 n = PetscCeilInt64(10, 3); // n has the value of 4
1203 .ve
1204 
1205 .seealso: `PetscCeilInt()`, `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1206 @*/
1207 static inline PetscInt64 PetscCeilInt64(PetscInt64 x, PetscInt64 y)
1208 {
1209   return x / y + (x % y ? 1 : 0);
1210 }
1211 
1212 PETSC_EXTERN PetscErrorCode PetscLinearRegression(PetscInt, const PetscReal[], const PetscReal[], PetscReal *, PetscReal *);
1213