xref: /petsc/include/petscmath.h (revision 3923b477fd0dced8a2d147b4fb4519fe3af97d3f)
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 
11 #if !defined(__PETSCMATH_H)
12 #define __PETSCMATH_H
13 #include <math.h>
14 
15 PETSC_EXTERN MPI_Datatype MPIU_2SCALAR;
16 PETSC_EXTERN MPI_Datatype MPIU_2INT;
17 
18 /*
19 
20      Defines operations that are different for complex and real numbers;
21    note that one cannot mix the use of complex and real in the same
22    PETSc program. All PETSc objects in one program are built around the object
23    PetscScalar which is either always a real or a complex.
24 
25 */
26 
27 #define PetscExpPassiveScalar(a) PetscExpScalar()
28 #if defined(PETSC_USE_REAL_SINGLE)
29 #define MPIU_REAL   MPI_FLOAT
30 typedef float PetscReal;
31 #define PetscSqrtReal(a)    sqrt(a)
32 #elif defined(PETSC_USE_REAL_DOUBLE)
33 #define MPIU_REAL   MPI_DOUBLE
34 typedef double PetscReal;
35 #define PetscSqrtReal(a)    sqrt(a)
36 #elif defined(PETSC_USE_REAL___FLOAT128)
37 #if defined(__cplusplus)
38 extern "C" {
39 #endif
40 #include <quadmath.h>
41 #if defined(__cplusplus)
42 }
43 #endif
44 #define MPIU_REAL MPIU___FLOAT128
45 typedef __float128 PetscReal;
46 #define PetscSqrtReal(a)    sqrtq(a)
47 #endif /* PETSC_USE_REAL_* */
48 
49 /*
50     Complex number definitions
51  */
52 #if defined(PETSC_CLANGUAGE_CXX) && defined(PETSC_HAVE_CXX_COMPLEX)
53 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_DESIRE_COMPLEX)
54 #define PETSC_HAVE_COMPLEX 1
55 /* C++ support of complex number */
56 #if defined(PETSC_HAVE_CUSP)
57 #define complexlib cusp
58 #include <cusp/complex.h>
59 #else
60 #define complexlib std
61 #include <complex>
62 #endif
63 
64 #define PetscRealPartComplex(a)      (a).real()
65 #define PetscImaginaryPartComplex(a) (a).imag()
66 #define PetscAbsComplex(a)           complexlib::abs(a)
67 #define PetscConjComplex(a)          complexlib::conj(a)
68 #define PetscSqrtComplex(a)          complexlib::sqrt(a)
69 #define PetscPowComplex(a,b)         complexlib::pow(a,b)
70 #define PetscExpComplex(a)           complexlib::exp(a)
71 #define PetscLogComplex(a)           complexlib::log(a)
72 #define PetscSinComplex(a)           complexlib::sin(a)
73 #define PetscCosComplex(a)           complexlib::cos(a)
74 
75 #if defined(PETSC_USE_REAL_SINGLE)
76 typedef complexlib::complex<float> PetscComplex;
77 #elif defined(PETSC_USE_REAL_DOUBLE)
78 typedef complexlib::complex<double> PetscComplex;
79 #elif defined(PETSC_USE_REAL___FLOAT128)
80 typedef complexlib::complex<__float128> PetscComplex; /* Notstandard and not expected to work, use __complex128 */
81 #endif  /* PETSC_USE_REAL_ */
82 #endif  /* PETSC_USE_COMPLEX && PETSC_DESIRE_COMPLEX */
83 
84 #elif defined(PETSC_CLANGUAGE_C) && defined(PETSC_HAVE_C99_COMPLEX)
85 /* Use C99 _Complex for the type. Do not include complex.h by default to define "complex" because of symbol conflicts in Hypre. */
86 /* Compilation units that can safely use complex should define PETSC_DESIRE_COMPLEX before including any headers */
87 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_DESIRE_COMPLEX)
88 #define PETSC_HAVE_COMPLEX 1
89 #include <complex.h>
90 
91 #if defined(PETSC_USE_REAL_SINGLE)
92 typedef float _Complex PetscComplex;
93 
94 #define PetscRealPartComplex(a)      crealf(a)
95 #define PetscImaginaryPartComplex(a) cimagf(a)
96 #define PetscAbsComplex(a)           cabsf(a)
97 #define PetscConjComplex(a)          conjf(a)
98 #define PetscSqrtComplex(a)          csqrtf(a)
99 #define PetscPowComplex(a,b)         cpowf(a,b)
100 #define PetscExpComplex(a)           cexpf(a)
101 #define PetscLogComplex(a)           clogf(a)
102 #define PetscSinComplex(a)           csinf(a)
103 #define PetscCosComplex(a)           ccosf(a)
104 
105 #elif defined(PETSC_USE_REAL_DOUBLE)
106 typedef double _Complex PetscComplex;
107 
108 #define PetscRealPartComplex(a)      creal(a)
109 #define PetscImaginaryPartComplex(a) cimag(a)
110 #define PetscAbsComplex(a)           cabs(a)
111 #define PetscConjComplex(a)          conj(a)
112 #define PetscSqrtComplex(a)          csqrt(a)
113 #define PetscPowComplex(a,b)         cpow(a,b)
114 #define PetscExpComplex(a)           cexp(a)
115 #define PetscLogComplex(a)           clog(a)
116 #define PetscSinComplex(a)           csin(a)
117 #define PetscCosComplex(a)           ccos(a)
118 
119 #elif defined(PETSC_USE_REAL___FLOAT128)
120 typedef __complex128 PetscComplex;
121 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128;
122 PETSC_EXTERN MPI_Datatype MPIU___COMPLEX128;
123 
124 #define PetscRealPartComplex(a)      crealq(a)
125 #define PetscImaginaryPartComplex(a) cimagq(a)
126 #define PetscAbsComplex(a)           cabsq(a)
127 #define PetscConjComplex(a)          conjq(a)
128 #define PetscSqrtComplex(a)          csqrtq(a)
129 #define PetscPowComplex(a,b)         cpowq(a,b)
130 #define PetscExpComplex(a)           cexpq(a)
131 #define PetscLogComplex(a)           clogq(a)
132 #define PetscSinComplex(a)           csinq(a)
133 #define PetscCosComplex(a)           ccosq(a)
134 #endif /* PETSC_USE_REAL_* */
135 #elif defined(PETSC_USE_COMPLEX)
136 #error "PETSc was configured --with-scalar-type=complex, but a language-appropriate complex library is not available"
137 #endif /* PETSC_USE_COMPLEX || PETSC_DESIRE_COMPLEX */
138 #endif /* (PETSC_CLANGUAGE_CXX && PETSC_HAVE_CXX_COMPLEX) else-if (PETSC_CLANGUAGE_C && PETSC_HAVE_C99_COMPLEX) */
139 
140 #if defined(PETSC_HAVE_MPI_C_DOUBLE_COMPLEX)
141 #define MPIU_C_DOUBLE_COMPLEX MPI_C_DOUBLE_COMPLEX
142 #define MPIU_C_COMPLEX MPI_C_COMPLEX
143 #else
144 PETSC_EXTERN MPI_Datatype MPIU_C_DOUBLE_COMPLEX;
145 PETSC_EXTERN MPI_Datatype MPIU_C_COMPLEX;
146 #endif /* PETSC_HAVE_MPI_C_DOUBLE_COMPLEX */
147 
148 #if defined(PETSC_USE_COMPLEX)
149 typedef PetscComplex PetscScalar;
150 #define PetscRealPart(a)      PetscRealPartComplex(a)
151 #define PetscImaginaryPart(a) PetscImaginaryPartComplex(a)
152 #define PetscAbsScalar(a)     PetscAbsComplex(a)
153 #define PetscConj(a)          PetscConjComplex(a)
154 #define PetscSqrtScalar(a)    PetscSqrtComplex(a)
155 #define PetscPowScalar(a,b)   PetscPowComplex(a,b)
156 #define PetscExpScalar(a)     PetscExpComplex(a)
157 #define PetscLogScalar(a)     PetscLogComplex(a)
158 #define PetscSinScalar(a)     PetscSinComplex(a)
159 #define PetscCosScalar(a)     PetscCosComplex(a)
160 
161 #if defined(PETSC_USE_REAL_SINGLE)
162 #define MPIU_SCALAR MPIU_C_COMPLEX
163 #elif defined(PETSC_USE_REAL_DOUBLE)
164 #define MPIU_SCALAR MPIU_C_DOUBLE_COMPLEX
165 #elif defined(PETSC_USE_REAL___FLOAT128)
166 #define MPIU_SCALAR MPIU___COMPLEX128
167 #endif /* PETSC_USE_REAL_* */
168 
169 /*
170     real number definitions
171  */
172 #else /* PETSC_USE_COMPLEX */
173 #if defined(PETSC_USE_REAL_SINGLE)
174 #define MPIU_SCALAR           MPI_FLOAT
175 typedef float PetscScalar;
176 #elif defined(PETSC_USE_REAL_DOUBLE)
177 #define MPIU_SCALAR           MPI_DOUBLE
178 typedef double PetscScalar;
179 #elif defined(PETSC_USE_REAL___FLOAT128)
180 PETSC_EXTERN MPI_Datatype MPIU___FLOAT128;
181 #define MPIU_SCALAR MPIU___FLOAT128
182 typedef __float128 PetscScalar;
183 #endif /* PETSC_USE_REAL_* */
184 #define PetscRealPart(a)      (a)
185 #define PetscImaginaryPart(a) ((PetscReal)0.)
186 PETSC_STATIC_INLINE PetscReal PetscAbsScalar(PetscScalar a) {return a < 0.0 ? -a : a;}
187 #define PetscConj(a)          (a)
188 #if !defined(PETSC_USE_REAL___FLOAT128)
189 #define PetscSqrtScalar(a)    sqrt(a)
190 #define PetscPowScalar(a,b)   pow(a,b)
191 #define PetscExpScalar(a)     exp(a)
192 #define PetscLogScalar(a)     log(a)
193 #define PetscSinScalar(a)     sin(a)
194 #define PetscCosScalar(a)     cos(a)
195 #else /* PETSC_USE_REAL___FLOAT128 */
196 #define PetscSqrtScalar(a)    sqrtq(a)
197 #define PetscPowScalar(a,b)   powq(a,b)
198 #define PetscExpScalar(a)     expq(a)
199 #define PetscLogScalar(a)     logq(a)
200 #define PetscSinScalar(a)     sinq(a)
201 #define PetscCosScalar(a)     cosq(a)
202 #endif /* PETSC_USE_REAL___FLOAT128 */
203 
204 #endif /* PETSC_USE_COMPLEX */
205 
206 #define PetscSign(a) (((a) >= 0) ? ((a) == 0 ? 0 : 1) : -1)
207 #define PetscAbs(a)  (((a) >= 0) ? (a) : -(a))
208 
209 /* --------------------------------------------------------------------------*/
210 
211 /*
212    Certain objects may be created using either single or double precision.
213    This is currently not used.
214 */
215 typedef enum { PETSC_SCALAR_DOUBLE,PETSC_SCALAR_SINGLE, PETSC_SCALAR_LONG_DOUBLE } PetscScalarPrecision;
216 
217 #if defined(PETSC_HAVE_COMPLEX)
218 /* PETSC_i is the imaginary number, i */
219 PETSC_EXTERN PetscComplex PETSC_i;
220 #endif
221 
222 /*MC
223    PetscMin - Returns minimum of two numbers
224 
225    Synopsis:
226    type PetscMin(type v1,type v2)
227 
228    Not Collective
229 
230    Input Parameter:
231 +  v1 - first value to find minimum of
232 -  v2 - second value to find minimum of
233 
234 
235    Notes: type can be integer or floating point value
236 
237    Level: beginner
238 
239 
240 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
241 
242 M*/
243 #define PetscMin(a,b)   (((a)<(b)) ?  (a) : (b))
244 
245 /*MC
246    PetscMax - Returns maxium of two numbers
247 
248    Synopsis:
249    type max PetscMax(type v1,type v2)
250 
251    Not Collective
252 
253    Input Parameter:
254 +  v1 - first value to find maximum of
255 -  v2 - second value to find maximum of
256 
257    Notes: type can be integer or floating point value
258 
259    Level: beginner
260 
261 .seealso: PetscMin(), PetscClipInterval(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
262 
263 M*/
264 #define PetscMax(a,b)   (((a)<(b)) ?  (b) : (a))
265 
266 /*MC
267    PetscClipInterval - Returns a number clipped to be within an interval
268 
269    Synopsis:
270    type clip PetscClipInterval(type x,type a,type b)
271 
272    Not Collective
273 
274    Input Parameter:
275 +  x - value to use if within interval (a,b)
276 .  a - lower end of interval
277 -  b - upper end of interval
278 
279    Notes: type can be integer or floating point value
280 
281    Level: beginner
282 
283 .seealso: PetscMin(), PetscMax(), PetscAbsInt(), PetscAbsReal(), PetscSqr()
284 
285 M*/
286 #define PetscClipInterval(x,a,b)   (PetscMax((a),PetscMin((x),(b))))
287 
288 /*MC
289    PetscAbsInt - Returns the absolute value of an integer
290 
291    Synopsis:
292    int abs PetscAbsInt(int v1)
293 
294    Not Collective
295 
296    Input Parameter:
297 .   v1 - the integer
298 
299    Level: beginner
300 
301 .seealso: PetscMax(), PetscMin(), PetscAbsReal(), PetscSqr()
302 
303 M*/
304 #define PetscAbsInt(a)  (((a)<0)   ? -(a) : (a))
305 
306 /*MC
307    PetscAbsReal - Returns the absolute value of an real number
308 
309    Synopsis:
310    Real abs PetscAbsReal(PetscReal v1)
311 
312    Not Collective
313 
314    Input Parameter:
315 .   v1 - the double
316 
317 
318    Level: beginner
319 
320 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscSqr()
321 
322 M*/
323 #define PetscAbsReal(a) (((a)<0)   ? -(a) : (a))
324 
325 /*MC
326    PetscSqr - Returns the square of a number
327 
328    Synopsis:
329    type sqr PetscSqr(type v1)
330 
331    Not Collective
332 
333    Input Parameter:
334 .   v1 - the value
335 
336    Notes: type can be integer or floating point value
337 
338    Level: beginner
339 
340 .seealso: PetscMax(), PetscMin(), PetscAbsInt(), PetscAbsReal()
341 
342 M*/
343 #define PetscSqr(a)     ((a)*(a))
344 
345 /* ----------------------------------------------------------------------------*/
346 /*
347      Basic constants
348 */
349 #if defined(PETSC_USE_REAL___FLOAT128)
350 #define PETSC_PI                 M_PIq
351 #elif defined(M_PI)
352 #define PETSC_PI                 M_PI
353 #else
354 #define PETSC_PI                 3.14159265358979323846264338327950288419716939937510582
355 #endif
356 
357 #if !defined(PETSC_USE_64BIT_INDICES)
358 #define PETSC_MAX_INT            2147483647
359 #define PETSC_MIN_INT            (-PETSC_MAX_INT - 1)
360 #else
361 #define PETSC_MAX_INT            9223372036854775807L
362 #define PETSC_MIN_INT            (-PETSC_MAX_INT - 1)
363 #endif
364 
365 #if defined(PETSC_USE_REAL_SINGLE)
366 #  define PETSC_MAX_REAL                3.40282346638528860e+38F
367 #  define PETSC_MIN_REAL                -PETSC_MAX_REAL
368 #  define PETSC_MACHINE_EPSILON         1.19209290e-07F
369 #  define PETSC_SQRT_MACHINE_EPSILON    3.45266983e-04F
370 #  define PETSC_SMALL                   1.e-5
371 #elif defined(PETSC_USE_REAL_DOUBLE)
372 #  define PETSC_MAX_REAL                1.7976931348623157e+308
373 #  define PETSC_MIN_REAL                -PETSC_MAX_REAL
374 #  define PETSC_MACHINE_EPSILON         2.2204460492503131e-16
375 #  define PETSC_SQRT_MACHINE_EPSILON    1.490116119384766e-08
376 #  define PETSC_SMALL                   1.e-10
377 #elif defined(PETSC_USE_REAL___FLOAT128)
378 #  define PETSC_MAX_REAL                FLT128_MAX
379 #  define PETSC_MIN_REAL                -FLT128_MAX
380 #  define PETSC_MACHINE_EPSILON         FLT128_EPSILON
381 #  define PETSC_SQRT_MACHINE_EPSILON    1.38777878078e-17
382 #  define PETSC_SMALL                   1.e-20
383 #endif
384 
385 #if defined PETSC_HAVE_ADIC
386 /* Use MPI_Allreduce when ADIC is not available. */
387 PETSC_EXTERN PetscErrorCode PetscGlobalMax(MPI_Comm, const PetscReal*,PetscReal*);
388 PETSC_EXTERN PetscErrorCode PetscGlobalMin(MPI_Comm, const PetscReal*,PetscReal*);
389 PETSC_EXTERN PetscErrorCode PetscGlobalSum(MPI_Comm, const PetscScalar*,PetscScalar*);
390 #endif
391 
392 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanScalar(PetscScalar);
393 PETSC_EXTERN PetscErrorCode PetscIsInfOrNanReal(PetscReal);
394 
395 /* ----------------------------------------------------------------------------*/
396 #define PassiveReal   PetscReal
397 #define PassiveScalar PetscScalar
398 
399 /*
400     These macros are currently hardwired to match the regular data types, so there is no support for a different
401     MatScalar from PetscScalar. We left the MatScalar in the source just in case we use it again.
402  */
403 #define MPIU_MATSCALAR MPIU_SCALAR
404 typedef PetscScalar MatScalar;
405 typedef PetscReal MatReal;
406 
407 
408 #endif
409