xref: /petsc/include/petscmath.h (revision 66af8762ec03dbef0e079729eb2a1734a35ed7ff)
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 #define PETSC_INFINITY  (PETSC_MAX_REAL / 4)
786 #define PETSC_NINFINITY (-PETSC_INFINITY)
787 
788 PETSC_EXTERN PetscBool  PetscIsInfReal(PetscReal);
789 PETSC_EXTERN PetscBool  PetscIsNanReal(PetscReal);
790 PETSC_EXTERN PetscBool  PetscIsNormalReal(PetscReal);
791 static inline PetscBool PetscIsInfOrNanReal(PetscReal v)
792 {
793   return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;
794 }
795 static inline PetscBool PetscIsInfScalar(PetscScalar v)
796 {
797   return PetscIsInfReal(PetscAbsScalar(v));
798 }
799 static inline PetscBool PetscIsNanScalar(PetscScalar v)
800 {
801   return PetscIsNanReal(PetscAbsScalar(v));
802 }
803 static inline PetscBool PetscIsInfOrNanScalar(PetscScalar v)
804 {
805   return PetscIsInfOrNanReal(PetscAbsScalar(v));
806 }
807 static inline PetscBool PetscIsNormalScalar(PetscScalar v)
808 {
809   return PetscIsNormalReal(PetscAbsScalar(v));
810 }
811 
812 PETSC_EXTERN PetscBool PetscIsCloseAtTol(PetscReal, PetscReal, PetscReal, PetscReal);
813 PETSC_EXTERN PetscBool PetscEqualReal(PetscReal, PetscReal);
814 PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar, PetscScalar);
815 
816 /*@C
817   PetscIsCloseAtTolScalar - Like `PetscIsCloseAtTol()` but for `PetscScalar`
818 
819   Input Parameters:
820 + lhs  - The first number
821 . rhs  - The second number
822 . rtol - The relative tolerance
823 - atol - The absolute tolerance
824 
825   Level: beginner
826 
827   Note:
828   This routine is equivalent to `PetscIsCloseAtTol()` when PETSc is configured without complex
829   numbers.
830 
831 .seealso: `PetscIsCloseAtTol()`
832 @*/
833 static inline PetscBool PetscIsCloseAtTolScalar(PetscScalar lhs, PetscScalar rhs, PetscReal rtol, PetscReal atol)
834 {
835   PetscBool close = PetscIsCloseAtTol(PetscRealPart(lhs), PetscRealPart(rhs), rtol, atol);
836 
837   if (PetscDefined(USE_COMPLEX)) close = (PetscBool)(close && PetscIsCloseAtTol(PetscImaginaryPart(lhs), PetscImaginaryPart(rhs), rtol, atol));
838   return close;
839 }
840 
841 /*
842     These macros are currently hardwired to match the regular data types, so there is no support for a different
843     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
844  */
845 #define MPIU_MATSCALAR MPIU_SCALAR
846 typedef PetscScalar MatScalar;
847 typedef PetscReal   MatReal;
848 
849 struct petsc_mpiu_2scalar {
850   PetscScalar a, b;
851 };
852 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2scalar);
853 
854 /* MPI Datatypes for composite reductions */
855 struct petsc_mpiu_real_int {
856   PetscReal v;
857   PetscInt  i;
858 };
859 
860 struct petsc_mpiu_scalar_int {
861   PetscScalar v;
862   PetscInt    i;
863 };
864 
865 PETSC_EXTERN MPI_Datatype MPIU_REAL_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_real_int);
866 PETSC_EXTERN MPI_Datatype MPIU_SCALAR_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_scalar_int);
867 
868 #if defined(PETSC_USE_64BIT_INDICES)
869 struct /* __attribute__((packed, aligned(alignof(PetscInt *)))) */ petsc_mpiu_2int {
870   PetscInt a;
871   PetscInt b;
872 };
873 /*
874  static_assert(sizeof(struct petsc_mpiu_2int) == 2 * sizeof(PetscInt), "");
875  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt *), "");
876  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt[2]), "");
877 
878  clang generates warnings that petsc_mpiu_2int is not layout compatible with PetscInt[2] or
879  PetscInt *, even though (with everything else uncommented) both of the static_asserts above
880  pass! So we just comment it out...
881 */
882 PETSC_EXTERN MPI_Datatype MPIU_2INT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2int) */;
883 #else
884   #define MPIU_2INT MPI_2INT
885 #endif
886 PETSC_EXTERN MPI_Datatype MPI_4INT;
887 PETSC_EXTERN MPI_Datatype MPIU_4INT;
888 
889 static inline PetscInt PetscPowInt(PetscInt base, PetscInt power)
890 {
891   PetscInt result = 1;
892   while (power) {
893     if (power & 1) result *= base;
894     power >>= 1;
895     if (power) base *= base;
896   }
897   return result;
898 }
899 
900 static inline PetscInt64 PetscPowInt64(PetscInt base, PetscInt power)
901 {
902   PetscInt64 result = 1;
903   while (power) {
904     if (power & 1) result *= base;
905     power >>= 1;
906     if (power) base *= base;
907   }
908   return result;
909 }
910 
911 static inline PetscReal PetscPowRealInt(PetscReal base, PetscInt power)
912 {
913   PetscReal result = 1;
914   if (power < 0) {
915     power = -power;
916     base  = ((PetscReal)1) / base;
917   }
918   while (power) {
919     if (power & 1) result *= base;
920     power >>= 1;
921     if (power) base *= base;
922   }
923   return result;
924 }
925 
926 static inline PetscScalar PetscPowScalarInt(PetscScalar base, PetscInt power)
927 {
928   PetscScalar result = (PetscReal)1;
929   if (power < 0) {
930     power = -power;
931     base  = ((PetscReal)1) / base;
932   }
933   while (power) {
934     if (power & 1) result *= base;
935     power >>= 1;
936     if (power) base *= base;
937   }
938   return result;
939 }
940 
941 static inline PetscScalar PetscPowScalarReal(PetscScalar base, PetscReal power)
942 {
943   PetscScalar cpower = power;
944   return PetscPowScalar(base, cpower);
945 }
946 
947 /*MC
948     PetscApproximateLTE - Performs a less than or equal to on a given constant with a fudge for floating point numbers
949 
950    Synopsis:
951    #include <petscmath.h>
952    bool PetscApproximateLTE(PetscReal x,constant float)
953 
954    Not Collective
955 
956    Input Parameters:
957 +   x - the variable
958 -   b - the constant float it is checking if `x` is less than or equal to
959 
960    Level: advanced
961 
962    Notes:
963      The fudge factor is the value `PETSC_SMALL`
964 
965      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
966 
967      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
968      floating point results.
969 
970 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateGTE()`
971 M*/
972 #define PetscApproximateLTE(x, b) ((x) <= (PetscRealConstant(b) + PETSC_SMALL))
973 
974 /*MC
975     PetscApproximateGTE - Performs a greater than or equal to on a given constant with a fudge for floating point numbers
976 
977    Synopsis:
978    #include <petscmath.h>
979    bool PetscApproximateGTE(PetscReal x,constant float)
980 
981    Not Collective
982 
983    Input Parameters:
984 +   x - the variable
985 -   b - the constant float it is checking if `x` is greater than or equal to
986 
987    Level: advanced
988 
989    Notes:
990      The fudge factor is the value `PETSC_SMALL`
991 
992      The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2
993 
994      This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
995      floating point results.
996 
997 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
998 M*/
999 #define PetscApproximateGTE(x, b) ((x) >= (PetscRealConstant(b) - PETSC_SMALL))
1000 
1001 /*MC
1002     PetscCeilInt - Returns the ceiling of the quotation of two positive integers
1003 
1004    Synopsis:
1005    #include <petscmath.h>
1006    PetscInt PetscCeilInt(PetscInt x,PetscInt y)
1007 
1008    Not Collective
1009 
1010    Input Parameters:
1011 +   x - the numerator
1012 -   y - the denominator
1013 
1014    Level: advanced
1015 
1016 .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1017 M*/
1018 #define PetscCeilInt(x, y) ((((PetscInt)(x)) / ((PetscInt)(y))) + ((((PetscInt)(x)) % ((PetscInt)(y))) ? 1 : 0))
1019 
1020 #define PetscCeilInt64(x, y) ((((PetscInt64)(x)) / ((PetscInt64)(y))) + ((((PetscInt64)(x)) % ((PetscInt64)(y))) ? 1 : 0))
1021 
1022 PETSC_EXTERN PetscErrorCode PetscLinearRegression(PetscInt, const PetscReal[], const PetscReal[], PetscReal *, PetscReal *);
1023