The R Project SVN R-packages

/* ========================================================================== */
/* === Core/cholmod_complex ================================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Core Module.  Copyright (C) 2005-2006,
 * Univ. of Florida.  Author: Timothy A. Davis
 * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
 * Lesser General Public License.  See lesser.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* If you convert a matrix that contains uninitialized data, valgrind will
 * complain.  This can occur in a factor L which has gaps (a partial
 * factorization, or after updates that change the nonzero pattern), an
 * unpacked sparse matrix, a dense matrix with leading dimension d > # of rows,
 * or any matrix (dense, sparse, triplet, or factor) with more space allocated
 * than is used.  You can safely ignore any of these complaints by valgrind. */

#include "cholmod_internal.h"
#include "cholmod_core.h"

/* ========================================================================== */
/* === cholmod_hypot ======================================================== */
/* ========================================================================== */

/* There is an equivalent routine called hypot in <math.h>, which conforms
 * to ANSI C99.  However, CHOLMOD does not assume that ANSI C99 is available.
 * You can use the ANSI C99 hypot routine with:
 *
 *  #include <math.h>
 *  Common->hypotenuse = hypot ;
 *
 * Default value of the Common->hypotenuse pointer is cholmod_hypot.
 *
 * s = hypot (x,y) computes s = sqrt (x*x + y*y) but does so more accurately.
 * The NaN cases for the double relops x >= y and x+y == x are safely ignored.
 */

double CHOLMOD(hypot) (double x, double y)
{
    double s, r ;
    x = fabs (x) ;
    y = fabs (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) ;
}


/* ========================================================================== */
/* === cholmod_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 1 if the denominator is zero,
 * 0 otherwise.  Note that this return value is the single exception to the
 * rule that all CHOLMOD routines that return int return TRUE if successful
 * or FALSE otherise.
 *
 * 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.
 *
 * Default value of the Common->complex_divide pointer is cholmod_divcomplex.
 */

int CHOLMOD(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 (fabs (br) >= fabs (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 (IS_ZERO (den)) ;
}


/* ========================================================================== */
/* === change_complexity ==================================================== */
/* ========================================================================== */

/* X and Z represent an array of size nz, with numeric xtype given by xtype_in.
 *
 * If xtype_in is:
 * CHOLMOD_PATTERN: X and Z must be NULL.
 * CHOLMOD_REAL:    X is of size nz, Z must be NULL.
 * CHOLMOD_COMPLEX: X is of size 2*nz, Z must be NULL.
 * CHOLMOD_ZOMPLEX: X is of size nz, Z is of size nz.
 *
 * The array is changed into the numeric xtype given by xtype_out, with the
 * same definitions of X and Z above.  Note that the input conditions, above,
 * are not checked.  These are checked in the caller routine.
 *
 * Returns TRUE if successful, FALSE otherwise.  X and Z are not modified if
 * not successful.
 */

static int change_complexity
(
    /* ---- input ---- */
    Int nz,     /* size of X and/or Z */
    int xtype_in,   /* xtype of X and Z on input */
    int xtype_out,  /* requested xtype of X and Z on output */
    int xtype1,     /* xtype_out must be in the range [xtype1 .. xtype2] */
    int xtype2,
    /* ---- in/out --- */
    void **XX,      /* old X on input, new X on output */
    void **ZZ,      /* old Z on input, new Z on output */
    /* --------------- */
    cholmod_common *Common
)
{
    double *Xold, *Zold, *Xnew, *Znew ;
    Int k ;
    size_t nz2 ;

    if (xtype_out < xtype1 || xtype_out > xtype2)
    {
    ERROR (CHOLMOD_INVALID, "invalid xtype") ;
    return (FALSE) ;
    }

    Common->status = CHOLMOD_OK ;
    Xold = *XX ;
    Zold = *ZZ ;

    switch (xtype_in)
    {

    /* ------------------------------------------------------------------ */
    /* converting from pattern */
    /* ------------------------------------------------------------------ */

    case CHOLMOD_PATTERN:

        switch (xtype_out)
        {

        /* ---------------------------------------------------------- */
        /* pattern -> real */
        /* ---------------------------------------------------------- */

        case CHOLMOD_REAL:
            /* allocate X and set to all ones */
            Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [k] = 1 ;
            }
            *XX = Xnew ;
            break ;

        /* ---------------------------------------------------------- */
        /* pattern -> complex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_COMPLEX:
            /* allocate X and set to all ones */
            Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [2*k  ] = 1 ;
            Xnew [2*k+1] = 0 ;
            }
            *XX = Xnew ;
            break ;

        /* ---------------------------------------------------------- */
        /* pattern -> zomplex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_ZOMPLEX:
            /* allocate X and Z and set to all ones */
            Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
            CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [k] = 1 ;
            Znew [k] = 0 ;
            }
            *XX = Xnew ;
            *ZZ = Znew ;
            break ;
        }
        break ;

    /* ------------------------------------------------------------------ */
    /* converting from real */
    /* ------------------------------------------------------------------ */

    case CHOLMOD_REAL:

        switch (xtype_out)
        {

        /* ---------------------------------------------------------- */
        /* real -> pattern */
        /* ---------------------------------------------------------- */

        case CHOLMOD_PATTERN:
            /* free X */
            *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
            break ;

        /* ---------------------------------------------------------- */
        /* real -> complex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_COMPLEX:
            /* allocate a new X and copy the old X */
            Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [2*k  ] = Xold [k] ;
            Xnew [2*k+1] = 0 ;
            }
            CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
            *XX = Xnew ;
            break ;

        /* ---------------------------------------------------------- */
        /* real -> zomplex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_ZOMPLEX:
            /* allocate a new Z and set it to zero */
            Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Znew [k] = 0 ;
            }
            *ZZ = Znew ;
            break ;
        }
        break ;

    /* ------------------------------------------------------------------ */
    /* converting from complex */
    /* ------------------------------------------------------------------ */

    case CHOLMOD_COMPLEX:

        switch (xtype_out)
        {

        /* ---------------------------------------------------------- */
        /* complex -> pattern */
        /* ---------------------------------------------------------- */

        case CHOLMOD_PATTERN:
            /* free X */
            *XX = CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
            break ;

        /* ---------------------------------------------------------- */
        /* complex -> real */
        /* ---------------------------------------------------------- */

        case CHOLMOD_REAL:
            /* pack the real part of X, discarding the imaginary part */
            for (k = 0 ; k < nz ; k++)
            {
            Xold [k] = Xold [2*k] ;
            }
            /* shrink X in half (this cannot fail) */
            nz2 = 2*nz ;
            *XX = CHOLMOD(realloc) (nz, sizeof (double), *XX, &nz2,
                Common) ;
            break ;

        /* ---------------------------------------------------------- */
        /* complex -> zomplex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_ZOMPLEX:
            /* allocate X and Z and copy the old X into them */
            Xnew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            Znew = CHOLMOD(malloc) (nz, sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            CHOLMOD(free) (nz, sizeof (double), Xnew, Common) ;
            CHOLMOD(free) (nz, sizeof (double), Znew, Common) ;
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [k] = Xold [2*k  ] ;
            Znew [k] = Xold [2*k+1] ;
            }
            CHOLMOD(free) (nz, 2*sizeof (double), *XX, Common) ;
            *XX = Xnew ;
            *ZZ = Znew ;
            break ;
        }
        break ;

    /* ------------------------------------------------------------------ */
    /* converting from zomplex */
    /* ------------------------------------------------------------------ */

    case CHOLMOD_ZOMPLEX:

        switch (xtype_out)
        {

        /* ---------------------------------------------------------- */
        /* zomplex -> pattern */
        /* ---------------------------------------------------------- */

        case CHOLMOD_PATTERN:
            /* free X and Z */
            *XX = CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
            *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
            break ;

        /* ---------------------------------------------------------- */
        /* zomplex -> real */
        /* ---------------------------------------------------------- */

        case CHOLMOD_REAL:
            /* free the imaginary part */
            *ZZ = CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
            break ;

        /* ---------------------------------------------------------- */
        /* zomplex -> complex */
        /* ---------------------------------------------------------- */

        case CHOLMOD_COMPLEX:
            Xnew = CHOLMOD(malloc) (nz, 2*sizeof (double), Common) ;
            if (Common->status < CHOLMOD_OK)
            {
            return (FALSE) ;
            }
            for (k = 0 ; k < nz ; k++)
            {
            Xnew [2*k  ] = Xold [k] ;
            Xnew [2*k+1] = Zold [k] ;
            }
            CHOLMOD(free) (nz, sizeof (double), *XX, Common) ;
            CHOLMOD(free) (nz, sizeof (double), *ZZ, Common) ;
            *XX = Xnew ;
            *ZZ = NULL ;
            break ;

        }
        break ;
    }

    return (TRUE) ;
}


/* ========================================================================== */
/* === cholmod_sparse_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a sparse matrix.  Supports any type on input
 * and output (pattern, real, complex, or zomplex). */

int CHOLMOD(sparse_xtype)
(
    /* ---- input ---- */
    int to_xtype,   /* requested xtype */
    /* ---- in/out --- */
    cholmod_sparse *A,  /* sparse matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;

    ok = change_complexity (A->nzmax, A->xtype, to_xtype,
        CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(A->x), &(A->z), Common) ;
    if (ok)
    {
    A->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_triplet_xtype ================================================ */
/* ========================================================================== */

/* Change the numeric xtype of a triplet matrix.  Supports any type on input
 * and output (pattern, real, complex, or zomplex). */

int CHOLMOD(triplet_xtype)
(
    /* ---- input ---- */
    int to_xtype,   /* requested xtype */
    /* ---- in/out --- */
    cholmod_triplet *T, /* triplet matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (T, FALSE) ;
    RETURN_IF_XTYPE_INVALID (T, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    ok = change_complexity (T->nzmax, T->xtype, to_xtype,
        CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, &(T->x), &(T->z), Common) ;
    if (ok)
    {
    T->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_dense_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a dense matrix.  Supports real, complex or
 * zomplex on input and output */

int CHOLMOD(dense_xtype)
(
    /* ---- input ---- */
    int to_xtype,   /* requested xtype */
    /* ---- in/out --- */
    cholmod_dense *X,   /* dense matrix to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (X, FALSE) ;
    RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    ok = change_complexity (X->nzmax, X->xtype, to_xtype,
        CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(X->x), &(X->z), Common) ;
    if (ok)
    {
    X->xtype = to_xtype ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_factor_xtype ================================================= */
/* ========================================================================== */

/* Change the numeric xtype of a factor.  Supports real, complex or zomplex on
 * input and output.   Supernodal zomplex factors are not supported. */

int CHOLMOD(factor_xtype)
(
    /* ---- input ---- */
    int to_xtype,   /* requested xtype */
    /* ---- in/out --- */
    cholmod_factor *L,  /* factor to change */
    /* --------------- */
    cholmod_common *Common
)
{
    Int ok ;
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (L, FALSE) ;
    RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    if (L->is_super &&
        (L->xtype == CHOLMOD_ZOMPLEX || to_xtype == CHOLMOD_ZOMPLEX))
    {
    ERROR (CHOLMOD_INVALID, "invalid xtype for supernodal L") ;
    return (FALSE) ;
    }
    ok = change_complexity ((L->is_super ? L->xsize : L->nzmax), L->xtype,
        to_xtype, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, &(L->x), &(L->z), Common) ;
    if (ok)
    {
    L->xtype = to_xtype ;
    }
    return (ok) ;
}