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