Blame view
fvn_sparse/UMFPACK/Source/umfpack_global.c
3.77 KB
422234dc3
git-svn-id: https...
|
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)) ;
}
|