The R Project SVN R-packages

/* ========================================================================== */
/* === UMF_global =========================================================== */
/* ========================================================================== */

/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.5, Copyright (c) 2005 by Timothy A. Davis.  CISE Dept,   */
/* 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/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)) ;
}