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