1 /* ========================================================================== */
2 /* === UMFPACK_global ======================================================= */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* -------------------------------------------------------------------------- */
10
11 /*
12 Global variables. UMFPACK uses these function pointers for several
13 user-redefinable functions. The amd_* functions are defined in
14 AMD/Source/amd_global.c.
15
16 Function pointer default for mexFunction
17 (see MATLAB/umfpackmex.c)
18 ---------------- ------- ---------------
19 amd_malloc malloc mxMalloc
20 amd_free free mxFree
21 amd_realloc realloc mxRealloc
22 amd_calloc calloc mxCalloc
23 amd_printf printf mexPrintf
24
25 umfpack_hypot umf_hypot umf_hypot
26 umfpack_divcomplex umf_divcomplex umf_divcomplex
27
28 This routine is compiled just once for all four versions of UMFPACK
29 (int/UF_long, double/complex).
30 */
31
32 #include "umf_internal.h"
33
34 double (*umfpack_hypot) (double, double) = umf_hypot ;
35 int (*umfpack_divcomplex) (double, double, double, double, double *, double *)
36 = umf_divcomplex ;
37
38
39 /* ========================================================================== */
40 /* === umf_hypot ============================================================ */
41 /* ========================================================================== */
42
43 /* There is an equivalent routine called hypot in <math.h>, which conforms
44 * to ANSI C99. However, UMFPACK does not assume that ANSI C99 is available.
45 * You can use the ANSI C99 hypot routine with:
46 *
47 * #include <math.h>
48 * umfpack_hypot = hypot ;
49 *
50 * prior to calling any UMFPACK routine.
51 *
52 * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
53 *
54 * The NaN case for the double relops x >= y and x+y == x is safely ignored.
55 */
56
umf_hypot(double x,double y)57 double umf_hypot (double x, double y)
58 {
59 double s, r ;
60 x = SCALAR_ABS (x) ;
61 y = SCALAR_ABS (y) ;
62 if (x >= y)
63 {
64 if (x + y == x)
65 {
66 s = x ;
67 }
68 else
69 {
70 r = y / x ;
71 s = x * sqrt (1.0 + r*r) ;
72 }
73 }
74 else
75 {
76 if (y + x == y)
77 {
78 s = y ;
79 }
80 else
81 {
82 r = x / y ;
83 s = y * sqrt (1.0 + r*r) ;
84 }
85 }
86 return (s) ;
87 }
88
89
90 /* ========================================================================== */
91 /* === umf_divcomplex ======================================================= */
92 /* ========================================================================== */
93
94 /* c = a/b where c, a, and b are complex. The real and imaginary parts are
95 * passed as separate arguments to this routine. The NaN case is ignored
96 * for the double relop br >= bi. Returns TRUE (1) if the denominator is
97 * zero, FALSE (0) otherwise.
98 *
99 * This uses ACM Algo 116, by R. L. Smith, 1962, which tries to avoid
100 * underflow and overflow.
101 *
102 * c can be the same variable as a or b.
103 */
104
umf_divcomplex(double ar,double ai,double br,double bi,double * cr,double * ci)105 int umf_divcomplex
106 (
107 double ar, double ai, /* real and imaginary parts of a */
108 double br, double bi, /* real and imaginary parts of b */
109 double *cr, double *ci /* real and imaginary parts of c */
110 )
111 {
112 double tr, ti, r, den ;
113 if (SCALAR_ABS (br) >= SCALAR_ABS (bi))
114 {
115 r = bi / br ;
116 den = br + r * bi ;
117 tr = (ar + ai * r) / den ;
118 ti = (ai - ar * r) / den ;
119 }
120 else
121 {
122 r = br / bi ;
123 den = r * br + bi ;
124 tr = (ar * r + ai) / den ;
125 ti = (ai * r - ar) / den ;
126 }
127 *cr = tr ;
128 *ci = ti ;
129 return (SCALAR_IS_ZERO (den)) ;
130 }
131