Actual source code: petscmath.h

  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.

  6:     This file is included by petscsys.h and should not be used directly.
  7: */
  8: #pragma once

 10: #include <math.h>
 11: #include <petscmacros.h>
 12: #include <petscsystypes.h>

 14: /* SUBMANSEC = Sys */

 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: */

 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

 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

 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

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

165: #endif /* PETSC_USE_REAL_* */

167: static inline PetscReal PetscSignReal(PetscReal a)
168: {
169:   return (PetscReal)((a < (PetscReal)0) ? -1 : ((a > (PetscReal)0) ? 1 : 0));
170: }

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

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

187: /*MC
188:    MPIU_REAL - Portable MPI datatype corresponding to `PetscReal` independent of what precision `PetscReal` is in

190:    Level: beginner

192:    Note:
193:    In MPI calls that require an MPI datatype that matches a `PetscReal` or array of `PetscReal` values, pass this value.

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_* */

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 */

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))

236:   /* TODO: Add configure tests

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

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

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

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

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

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

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

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

305: */

307:   #else /* C99 support of complex number */

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)

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)

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)

378:     #endif /* PETSC_USE_REAL_* */
379:   #endif   /* (__cplusplus) */

381: /*MC
382:    PETSC_i - the pure imaginary complex number i

384:    Level: intermediate

386: .seealso: `PetscComplex`, `PetscScalar`
387: M*/
388: PETSC_EXTERN PetscComplex PETSC_i;

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),

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."

408:        So type punning should be portable. */
409:     union
410:     {
411:       PetscComplex z;
412:       PetscReal    f[2];
413:     } uz;

415:     uz.f[0] = x;
416:     uz.f[1] = y;
417:     return uz.z;
418:   }
419:   #endif
420: }

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", )

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

433: PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128 MPIU___COMPLEX128_ATTR_TAG;

435:     #undef MPIU___COMPLEX128_ATTR_TAG
436:   #endif /* PETSC_HAVE_REAL___FLOAT128 */

438:   /*MC
439:    MPIU_COMPLEX - Portable MPI datatype corresponding to `PetscComplex` independent of the precision of `PetscComplex`

441:    Level: beginner

443:    Note:
444:    In MPI calls that require an MPI datatype that matches a `PetscComplex` or array of `PetscComplex` values, pass this value.

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_* */

458: #endif /* PETSC_HAVE_COMPLEX */

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`

467:    Level: beginner

469:    Note:
470:    In MPI calls that require an MPI datatype that matches a `PetscScalar` or array of `PetscScalar` values, pass this value.

472: .seealso: `PetscReal`, `PetscScalar`, `PetscComplex`, `PetscInt`, `MPIU_REAL`, `MPIU_COMPLEX`, `MPIU_INT`
473: M*/
474:   #define MPIU_SCALAR MPIU_COMPLEX

476:   /*MC
477:    PetscRealPart - Returns the real part of a `PetscScalar`

479:    Synopsis:
480: #include <petscmath.h>
481:    PetscReal PetscRealPart(PetscScalar v)

483:    Not Collective

485:    Input Parameter:
486: .  v - value to find the real part of

488:    Level: beginner

490: .seealso: `PetscScalar`, `PetscImaginaryPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
491: M*/
492:   #define PetscRealPart(a) PetscRealPartComplex(a)

494:   /*MC
495:    PetscImaginaryPart - Returns the imaginary part of a `PetscScalar`

497:    Synopsis:
498: #include <petscmath.h>
499:    PetscReal PetscImaginaryPart(PetscScalar v)

501:    Not Collective

503:    Input Parameter:
504: .  v - value to find the imaginary part of

506:    Level: beginner

508:    Note:
509:    If PETSc was configured for real numbers then this always returns the value 0

511: .seealso: `PetscScalar`, `PetscRealPart()`, `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
512: M*/
513:   #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a)

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)

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)

559: #endif /* PETSC_USE_COMPLEX */

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;

572: /*MC
573:    PetscAbs - Returns the absolute value of a number

575:    Synopsis:
576: #include <petscmath.h>
577:    type PetscAbs(type v)

579:    Not Collective

581:    Input Parameter:
582: .  v - the number

584:    Level: beginner

586:    Note:
587:    The type can be integer or real floating point value, but cannot be complex

589: .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`, `PetscSign()`
590: M*/
591: #define PetscAbs(a) (((a) >= 0) ? (a) : (-(a)))

593: /*MC
594:    PetscSign - Returns the sign of a number as an integer of value -1, 0, or 1

596:    Synopsis:
597: #include <petscmath.h>
598:    int PetscSign(type v)

600:    Not Collective

602:    Input Parameter:
603: .  v - the number

605:    Level: beginner

607:    Note:
608:    The type can be integer or real floating point value

610: .seealso: `PetscAbsInt()`, `PetscAbsReal()`, `PetscAbsScalar()`
611: M*/
612: #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1)

614: /*MC
615:    PetscMin - Returns minimum of two numbers

617:    Synopsis:
618: #include <petscmath.h>
619:    type PetscMin(type v1,type v2)

621:    Not Collective

623:    Input Parameters:
624: +  v1 - first value to find minimum of
625: -  v2 - second value to find minimum of

627:    Level: beginner

629:    Note:
630:    The type can be integer or floating point value, but cannot be complex

632: .seealso: `PetscMax()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
633: M*/
634: #define PetscMin(a, b) (((a) < (b)) ? (a) : (b))

636: /*MC
637:    PetscMax - Returns maximum of two numbers

639:    Synopsis:
640: #include <petscmath.h>
641:    type max PetscMax(type v1,type v2)

643:    Not Collective

645:    Input Parameters:
646: +  v1 - first value to find maximum of
647: -  v2 - second value to find maximum of

649:    Level: beginner

651:    Note:
652:    The type can be integer or floating point value

654: .seealso: `PetscMin()`, `PetscClipInterval()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
655: M*/
656: #define PetscMax(a, b) (((a) < (b)) ? (b) : (a))

658: /*MC
659:    PetscClipInterval - Returns a number clipped to be within an interval

661:    Synopsis:
662: #include <petscmath.h>
663:    type clip PetscClipInterval(type x,type a,type b)

665:    Not Collective

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

672:    Level: beginner

674:    Note:
675:    The type can be integer or floating point value

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

683: .seealso: `PetscMin()`, `PetscMax()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscSqr()`
684: M*/
685: #define PetscClipInterval(x, a, b) (PetscMax((a), PetscMin((x), (b))))

687: /*MC
688:    PetscAbsInt - Returns the absolute value of an integer

690:    Synopsis:
691: #include <petscmath.h>
692:    int abs PetscAbsInt(int v1)

694:    Input Parameter:
695: .   v1 - the integer

697:    Level: beginner

699: .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsReal()`, `PetscSqr()`
700: M*/
701: #define PetscAbsInt(a) (((a) < 0) ? (-(a)) : (a))

703: /*MC
704:    PetscAbsReal - Returns the absolute value of a real number

706:    Synopsis:
707: #include <petscmath.h>
708:    Real abs PetscAbsReal(PetscReal v1)

710:    Input Parameter:
711: .   v1 - the `PetscReal` value

713:    Level: beginner

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

727: /*MC
728:    PetscSqr - Returns the square of a number

730:    Synopsis:
731: #include <petscmath.h>
732:    type sqr PetscSqr(type v1)

734:    Not Collective

736:    Input Parameter:
737: .   v1 - the value

739:    Level: beginner

741:    Note:
742:    The type can be integer, floating point, or complex floating point

744: .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`
745: M*/
746: #define PetscSqr(a) ((a) * (a))

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

752:    Synopsis:
753: #include <petscmath.h>
754:    PetscReal PetscRealConstant(real_number)

756:    Not Collective

758:    Input Parameter:
759: .   v1 - the real number, for example 1.5

761:    Level: beginner

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

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

785: /*
786:      Basic constants
787: */
788: /*MC
789:   PETSC_PI - the value of $ \pi$ to the correct precision of `PetscReal`.

791:   Level: beginner

793: .seealso: `PetscReal`, `PETSC_PHI`, `PETSC_SQRT2`
794: M*/

796: /*MC
797:   PETSC_PHI - the value of $ \phi$, the Golden Ratio, to the correct precision of `PetscReal`.

799:   Level: beginner

801: .seealso: `PetscReal`, `PETSC_PI`, `PETSC_SQRT2`
802: M*/

804: /*MC
805:   PETSC_SQRT2 - the value of $ \sqrt{2} $ to the correct precision of `PetscReal`.

807:   Level: beginner

809: .seealso: `PetscReal`, `PETSC_PI`, `PETSC_PHI`
810: M*/

812: /*MC
813:   PETSC_E - the value of Euler's constant $ e $ to the correct precision of `PetscReal`.

815:   Level: beginner

817: .seealso: `PetscReal`, `PETSC_PI`, `PETSC_PHI`
818: M*/

820: #define PETSC_PI    PetscRealConstant(3.1415926535897932384626433832795029)
821: #define PETSC_PHI   PetscRealConstant(1.6180339887498948482045868343656381)
822: #define PETSC_SQRT2 PetscRealConstant(1.4142135623730950488016887242096981)
823: #define PETSC_E     PetscRealConstant(2.7182818284590452353602874713526625)

825: /*MC
826:   PETSC_MAX_REAL - the largest real value that can be stored in a `PetscReal`

828:   Level: beginner

830: .seealso: `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
831: M*/

833: /*MC
834:   PETSC_MIN_REAL - the smallest real value that can be stored in a `PetscReal`, generally this is - `PETSC_MAX_REAL`

836:   Level: beginner

838: .seealso `PETSC_MAX_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
839: M*/

841: /*MC
842:   PETSC_REAL_MIN - the smallest positive normalized real value that can be stored in a `PetscReal`.

844:   Level: beginner

846:   Note:
847:   See <https://en.wikipedia.org/wiki/Subnormal_number> for a discussion of normalized and subnormal floating point numbers

849:   Developer Note:
850:   The naming is confusing as there is both a `PETSC_REAL_MIN` and `PETSC_MIN_REAL` with different meanings.

852: .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_MACHINE_EPSILON`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
853: M*/

855: /*MC
856:   PETSC_MACHINE_EPSILON - the machine epsilon for the precision of `PetscReal`

858:   Level: beginner

860:   Note:
861:   See <https://en.wikipedia.org/wiki/Machine_epsilon>

863: .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_SQRT_MACHINE_EPSILON`, `PETSC_SMALL`
864: M*/

866: /*MC
867:   PETSC_SQRT_MACHINE_EPSILON - the square root of the machine epsilon for the precision of `PetscReal`

869:   Level: beginner

871:   Note:
872:   See `PETSC_MACHINE_EPSILON`

874: .seealso `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`, `PETSC_SMALL`
875: M*/

877: /*MC
878:   PETSC_SMALL - an arbitrary "small" number which depends on the precision of `PetscReal` used in some PETSc examples
879:   and in `PetscApproximateLTE()` and `PetscApproximateGTE()` to determine if a computation was successful.

881:   Level: beginner

883:   Note:
884:   See `PETSC_MACHINE_EPSILON`

886: .seealso `PetscApproximateLTE()`, `PetscApproximateGTE()`, `PETSC_MAX_REAL`, `PETSC_MIN_REAL`, `PETSC_REAL_MIN`, `PETSC_MACHINE_EPSILON`,
887:          `PETSC_SQRT_MACHINE_EPSILON`
888: M*/

890: #if defined(PETSC_USE_REAL_SINGLE)
891:   #define PETSC_MAX_REAL             3.40282346638528860e+38F
892:   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
893:   #define PETSC_REAL_MIN             1.1754944e-38F
894:   #define PETSC_MACHINE_EPSILON      1.19209290e-07F
895:   #define PETSC_SQRT_MACHINE_EPSILON 3.45266983e-04F
896:   #define PETSC_SMALL                1.e-5F
897: #elif defined(PETSC_USE_REAL_DOUBLE)
898:   #define PETSC_MAX_REAL             1.7976931348623157e+308
899:   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
900:   #define PETSC_REAL_MIN             2.225073858507201e-308
901:   #define PETSC_MACHINE_EPSILON      2.2204460492503131e-16
902:   #define PETSC_SQRT_MACHINE_EPSILON 1.490116119384766e-08
903:   #define PETSC_SMALL                1.e-10
904: #elif defined(PETSC_USE_REAL___FLOAT128)
905:   #define PETSC_MAX_REAL             FLT128_MAX
906:   #define PETSC_MIN_REAL             (-FLT128_MAX)
907:   #define PETSC_REAL_MIN             FLT128_MIN
908:   #define PETSC_MACHINE_EPSILON      FLT128_EPSILON
909:   #define PETSC_SQRT_MACHINE_EPSILON 1.38777878078144567552953958511352539e-17Q
910:   #define PETSC_SMALL                1.e-20Q
911: #elif defined(PETSC_USE_REAL___FP16)
912:   #define PETSC_MAX_REAL             65504.0F
913:   #define PETSC_MIN_REAL             (-PETSC_MAX_REAL)
914:   #define PETSC_REAL_MIN             .00006103515625F
915:   #define PETSC_MACHINE_EPSILON      .0009765625F
916:   #define PETSC_SQRT_MACHINE_EPSILON .03125F
917:   #define PETSC_SMALL                5.e-3F
918: #endif

920: /*MC
921:   PETSC_INFINITY - a finite number that represents infinity for setting certain bounds in `Tao`

923:   Level: intermediate

925:   Note:
926:   This is not the IEEE infinity value

928: .seealso: `PETSC_NINFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
929: M*/
930: #define PETSC_INFINITY (PETSC_MAX_REAL / 4)

932: /*MC
933:   PETSC_NINFINITY - a finite number that represents negative infinity for setting certain bounds in `Tao`

935:   Level: intermediate

937:   Note:
938:   This is not the negative IEEE infinity value

940: .seealso: `PETSC_INFINITY`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESVISetVariableBounds()`
941: M*/
942: #define PETSC_NINFINITY (-PETSC_INFINITY)

944: PETSC_EXTERN PetscBool  PetscIsInfReal(PetscReal);
945: PETSC_EXTERN PetscBool  PetscIsNanReal(PetscReal);
946: PETSC_EXTERN PetscBool  PetscIsNormalReal(PetscReal);
947: static inline PetscBool PetscIsInfOrNanReal(PetscReal v)
948: {
949:   return PetscIsInfReal(v) || PetscIsNanReal(v) ? PETSC_TRUE : PETSC_FALSE;
950: }
951: static inline PetscBool PetscIsInfScalar(PetscScalar v)
952: {
953:   return PetscIsInfReal(PetscAbsScalar(v));
954: }
955: static inline PetscBool PetscIsNanScalar(PetscScalar v)
956: {
957:   return PetscIsNanReal(PetscAbsScalar(v));
958: }
959: static inline PetscBool PetscIsInfOrNanScalar(PetscScalar v)
960: {
961:   return PetscIsInfOrNanReal(PetscAbsScalar(v));
962: }
963: static inline PetscBool PetscIsNormalScalar(PetscScalar v)
964: {
965:   return PetscIsNormalReal(PetscAbsScalar(v));
966: }

968: PETSC_EXTERN PetscBool PetscIsCloseAtTol(PetscReal, PetscReal, PetscReal, PetscReal);
969: PETSC_EXTERN PetscBool PetscEqualReal(PetscReal, PetscReal);
970: PETSC_EXTERN PetscBool PetscEqualScalar(PetscScalar, PetscScalar);

972: /*@C
973:   PetscIsCloseAtTolScalar - Like `PetscIsCloseAtTol()` but for `PetscScalar`

975:   Input Parameters:
976: + lhs  - The first number
977: . rhs  - The second number
978: . rtol - The relative tolerance
979: - atol - The absolute tolerance

981:   Level: beginner

983:   Note:
984:   This routine is equivalent to `PetscIsCloseAtTol()` when PETSc is configured without complex
985:   numbers.

987: .seealso: `PetscIsCloseAtTol()`
988: @*/
989: static inline PetscBool PetscIsCloseAtTolScalar(PetscScalar lhs, PetscScalar rhs, PetscReal rtol, PetscReal atol)
990: {
991:   PetscBool close = PetscIsCloseAtTol(PetscRealPart(lhs), PetscRealPart(rhs), rtol, atol);

993:   if (PetscDefined(USE_COMPLEX)) close = (PetscBool)(close && PetscIsCloseAtTol(PetscImaginaryPart(lhs), PetscImaginaryPart(rhs), rtol, atol));
994:   return close;
995: }

997: /*
998:     These macros are currently hardwired to match the regular data types, so there is no support for a different
999:     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
1000:  */
1001: #define MPIU_MATSCALAR MPIU_SCALAR
1002: typedef PetscScalar MatScalar;
1003: typedef PetscReal   MatReal;

1005: struct petsc_mpiu_2scalar {
1006:   PetscScalar a, b;
1007: };
1008: PETSC_EXTERN MPI_Datatype MPIU_2SCALAR PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2scalar);

1010: /* MPI Datatypes for composite reductions */
1011: struct petsc_mpiu_real_int {
1012:   PetscReal v;
1013:   PetscInt  i;
1014: };

1016: struct petsc_mpiu_scalar_int {
1017:   PetscScalar v;
1018:   PetscInt    i;
1019: };

1021: PETSC_EXTERN MPI_Datatype MPIU_REAL_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_real_int);
1022: PETSC_EXTERN MPI_Datatype MPIU_SCALAR_INT PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_scalar_int);

1024: #if defined(PETSC_USE_64BIT_INDICES)
1025: struct /* __attribute__((packed, aligned(alignof(PetscInt *)))) */ petsc_mpiu_2int {
1026:   PetscInt a;
1027:   PetscInt b;
1028: };
1029: struct __attribute__((packed)) petsc_mpiu_int_mpiint {
1030:   PetscInt    a;
1031:   PetscMPIInt b;
1032: };
1033: /*
1034:  static_assert(sizeof(struct petsc_mpiu_2int) == 2 * sizeof(PetscInt), "");
1035:  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt *), "");
1036:  static_assert(alignof(struct petsc_mpiu_2int) == alignof(PetscInt[2]), "");

1038:  clang generates warnings that petsc_mpiu_2int is not layout compatible with PetscInt[2] or
1039:  PetscInt *, even though (with everything else uncommented) both of the static_asserts above
1040:  pass! So we just comment it out...
1041: */
1042: PETSC_EXTERN MPI_Datatype MPIU_2INT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_2int) */;
1043: PETSC_EXTERN MPI_Datatype MPIU_INT_MPIINT /* PETSC_ATTRIBUTE_MPI_TYPE_TAG_LAYOUT_COMPATIBLE(struct petsc_mpiu_int_mpiint) */;
1044: #else
1045:   #define MPIU_2INT       MPI_2INT
1046:   #define MPIU_INT_MPIINT MPI_2INT
1047: #endif
1048: PETSC_EXTERN MPI_Datatype MPI_4INT;
1049: PETSC_EXTERN MPI_Datatype MPIU_4INT;

1051: static inline PetscInt PetscPowInt(PetscInt base, PetscInt power)
1052: {
1053:   PetscInt result = 1;
1054:   while (power) {
1055:     if (power & 1) result *= base;
1056:     power >>= 1;
1057:     if (power) base *= base;
1058:   }
1059:   return result;
1060: }

1062: static inline PetscInt64 PetscPowInt64(PetscInt base, PetscInt power)
1063: {
1064:   PetscInt64 result = 1;
1065:   while (power) {
1066:     if (power & 1) result *= base;
1067:     power >>= 1;
1068:     if (power) base *= base;
1069:   }
1070:   return result;
1071: }

1073: static inline PetscReal PetscPowRealInt(PetscReal base, PetscInt power)
1074: {
1075:   PetscReal result = 1;
1076:   if (power < 0) {
1077:     power = -power;
1078:     base  = ((PetscReal)1) / base;
1079:   }
1080:   while (power) {
1081:     if (power & 1) result *= base;
1082:     power >>= 1;
1083:     if (power) base *= base;
1084:   }
1085:   return result;
1086: }

1088: static inline PetscScalar PetscPowScalarInt(PetscScalar base, PetscInt power)
1089: {
1090:   PetscScalar result = (PetscReal)1;
1091:   if (power < 0) {
1092:     power = -power;
1093:     base  = ((PetscReal)1) / base;
1094:   }
1095:   while (power) {
1096:     if (power & 1) result *= base;
1097:     power >>= 1;
1098:     if (power) base *= base;
1099:   }
1100:   return result;
1101: }

1103: static inline PetscScalar PetscPowScalarReal(PetscScalar base, PetscReal power)
1104: {
1105:   PetscScalar cpower = power;
1106:   return PetscPowScalar(base, cpower);
1107: }

1109: /*MC
1110:    PetscApproximateLTE - Performs a less than or equal to on a given constant with a fudge for floating point numbers

1112:    Synopsis:
1113: #include <petscmath.h>
1114:    bool PetscApproximateLTE(PetscReal x,constant float)

1116:    Not Collective

1118:    Input Parameters:
1119: +   x - the variable
1120: -   b - the constant float it is checking if `x` is less than or equal to

1122:    Level: advanced

1124:    Notes:
1125:    The fudge factor is the value `PETSC_SMALL`

1127:    The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2

1129:    This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1130:    floating point results.

1132:    Example\:
1133: .vb
1134:   PetscReal x;
1135:   if (PetscApproximateLTE(x, 3.2)) { // replaces if (x <= 3.2) {
1136: .ve

1138: .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateGTE()`
1139: M*/
1140: #define PetscApproximateLTE(x, b) ((x) <= (PetscRealConstant(b) + PETSC_SMALL))

1142: /*MC
1143:    PetscApproximateGTE - Performs a greater than or equal to on a given constant with a fudge for floating point numbers

1145:    Synopsis:
1146: #include <petscmath.h>
1147:    bool PetscApproximateGTE(PetscReal x,constant float)

1149:    Not Collective

1151:    Input Parameters:
1152: +   x - the variable
1153: -   b - the constant float it is checking if `x` is greater than or equal to

1155:    Level: advanced

1157:    Notes:
1158:    The fudge factor is the value `PETSC_SMALL`

1160:    The constant numerical value is automatically set to the appropriate precision of PETSc so can just be provided as, for example, 3.2

1162:    This is used in several examples for setting initial conditions based on coordinate values that are computed with i*h that produces inexact
1163:    floating point results.

1165:    Example\:
1166: .vb
1167:   PetscReal x;
1168:   if (PetscApproximateGTE(x, 3.2)) {  // replaces if (x >= 3.2) {
1169: .ve

1171: .seealso: `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1172: M*/
1173: #define PetscApproximateGTE(x, b) ((x) >= (PetscRealConstant(b) - PETSC_SMALL))

1175: /*@C
1176:    PetscCeilInt - Returns the ceiling of the quotation of two positive integers

1178:    Not Collective

1180:    Input Parameters:
1181: +   x - the numerator
1182: -   y - the denominator

1184:    Level: advanced

1186:   Example\:
1187: .vb
1188:   PetscInt n = PetscCeilInt(10, 3); // n has the value of 4
1189: .ve

1191: .seealso: `PetscCeilInt64()`, `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1192: @*/
1193: static inline PetscInt PetscCeilInt(PetscInt x, PetscInt y)
1194: {
1195:   return x / y + (x % y ? 1 : 0);
1196: }

1198: /*@C
1199:    PetscCeilInt64 - Returns the ceiling of the quotation of two positive integers

1201:    Not Collective

1203:    Input Parameters:
1204: +   x - the numerator
1205: -   y - the denominator

1207:    Level: advanced

1209:   Example\:
1210: .vb
1211:   PetscInt64 n = PetscCeilInt64(10, 3); // n has the value of 4
1212: .ve

1214: .seealso: `PetscCeilInt()`, `PetscMax()`, `PetscMin()`, `PetscAbsInt()`, `PetscAbsReal()`, `PetscApproximateLTE()`
1215: @*/
1216: static inline PetscInt64 PetscCeilInt64(PetscInt64 x, PetscInt64 y)
1217: {
1218:   return x / y + (x % y ? 1 : 0);
1219: }

1221: PETSC_EXTERN PetscErrorCode PetscLinearRegression(PetscInt, const PetscReal[], const PetscReal[], PetscReal *, PetscReal *);