The R Project SVN R-packages

/* ========================================================================== */
/* === MatrixOps/cholmod_norm =============================================== */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/MatrixOps Module.  Copyright (C) 2005-2006, Timothy A. Davis
 * The CHOLMOD/MatrixOps Module is licensed under Version 2.0 of the GNU
 * General Public License.  See gpl.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
 * -------------------------------------------------------------------------- */

/* r = norm (A), compute the infinity-norm, 1-norm, or 2-norm of a sparse or
 * dense matrix.  Can compute the 2-norm only for a dense column vector.
 * Returns -1 if an error occurs.
 *
 * Pattern, real, complex, and zomplex sparse matrices are supported.
 */

#ifndef NMATRIXOPS

#include "cholmod_internal.h"
#include "cholmod_matrixops.h"


/* ========================================================================== */
/* === abs_value ============================================================ */
/* ========================================================================== */

/* Compute the absolute value of a real, complex, or zomplex value */

static double abs_value
(
    int xtype,
    double *Ax,
    double *Az,
    Int p,
    cholmod_common *Common
)
{
    double s = 0 ;
    switch (xtype)
    {
    case CHOLMOD_PATTERN:
        s = 1 ;
        break ;

    case CHOLMOD_REAL:
        s = fabs (Ax [p]) ;
        break ;

    case CHOLMOD_COMPLEX:
        s = Common->hypotenuse (Ax [2*p], Ax [2*p+1]) ;
        break ;

    case CHOLMOD_ZOMPLEX:
        s = Common->hypotenuse (Ax [p], Az [p]) ;
        break ;
    }
    return (s) ;
}


/* ========================================================================== */
/* === cholmod_norm_dense =================================================== */
/* ========================================================================== */

double CHOLMOD(norm_dense)
(
    /* ---- input ---- */
    cholmod_dense *X,   /* matrix to compute the norm of */
    int norm,       /* type of norm: 0: inf. norm, 1: 1-norm, 2: 2-norm */
    /* --------------- */
    cholmod_common *Common
)
{
    double xnorm, s, x, z ;
    double *Xx, *Xz, *W ;
    Int nrow, ncol, d, i, j, use_workspace, xtype ;

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (EMPTY) ;
    RETURN_IF_NULL (X, EMPTY) ;
    RETURN_IF_XTYPE_INVALID (X, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, EMPTY) ;
    Common->status = CHOLMOD_OK ;
    ncol = X->ncol ;
    if (norm < 0 || norm > 2 || (norm == 2 && ncol > 1))
    {
    ERROR (CHOLMOD_INVALID, "invalid norm") ;
    return (EMPTY) ;
    }

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    nrow = X->nrow ;
    d = X->d ;
    Xx = X->x ;
    Xz = X->z ;
    xtype = X->xtype ;

    /* ---------------------------------------------------------------------- */
    /* allocate workspace, if needed */
    /* ---------------------------------------------------------------------- */

    W = NULL ;
    use_workspace = (norm == 0 && ncol > 4) ;
    if (use_workspace)
    {
    CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
    W = Common->Xwork ;
    if (Common->status < CHOLMOD_OK)
    {
        /* oops, no workspace */
        use_workspace = FALSE ;
    }
    }


    /* ---------------------------------------------------------------------- */
    /* compute the norm */
    /* ---------------------------------------------------------------------- */

    xnorm = 0 ;

    if (use_workspace)
    {

    /* ------------------------------------------------------------------ */
    /* infinity-norm = max row sum, using stride-1 access of X */
    /* ------------------------------------------------------------------ */

    DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;

    /* this is faster than stride-d, but requires O(nrow) workspace */
    for (j = 0 ; j < ncol ; j++)
    {
        for (i = 0 ; i < nrow ; i++)
        {
        W [i] += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
        }
    }
    for (i = 0 ; i < nrow ; i++)
    {
        s = W [i] ;
        if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
        {
        xnorm = s ;
        }
        W [i] = 0 ;
    }

    }
    else if (norm == 0)
    {

    /* ------------------------------------------------------------------ */
    /* infinity-norm = max row sum, using stride-d access of X */
    /* ------------------------------------------------------------------ */

    for (i = 0 ; i < nrow ; i++)
    {
        s = 0 ;
        for (j = 0 ; j < ncol ; j++)
        {
        s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
        }
        if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
        {
        xnorm = s ;
        }
    }

    }
    else if (norm == 1)
    {

    /* ------------------------------------------------------------------ */
    /* 1-norm = max column sum */
    /* ------------------------------------------------------------------ */

    for (j = 0 ; j < ncol ; j++)
    {
        s = 0 ;
        for (i = 0 ; i < nrow ; i++)
        {
        s += abs_value (xtype, Xx, Xz, i+j*d, Common) ;
        }
        if ((IS_NAN (s) || s > xnorm) && !IS_NAN (xnorm))
        {
        xnorm = s ;
        }
    }
    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* 2-norm = sqrt (sum (X.^2)) */
    /* ------------------------------------------------------------------ */

    switch (xtype)
    {

        case CHOLMOD_REAL:
        for (i = 0 ; i < nrow ; i++)
        {
            x = Xx [i] ;
            xnorm += x*x ;
        }
        break ; 

        case CHOLMOD_COMPLEX:
        for (i = 0 ; i < nrow ; i++)
        {
            x = Xx [2*i  ] ;
            z = Xx [2*i+1] ;
            xnorm += x*x + z*z ;
        }
        break ; 

        case CHOLMOD_ZOMPLEX:
        for (i = 0 ; i < nrow ; i++)
        {
            x = Xx [i] ;
            z = Xz [i] ;
            xnorm += x*x + z*z ;
        }
        break ; 
    }

    xnorm = sqrt (xnorm) ;
    }

    /* ---------------------------------------------------------------------- */
    /* return result */
    /* ---------------------------------------------------------------------- */

    return (xnorm) ;
}


/* ========================================================================== */
/* === cholmod_norm_sparse ================================================== */
/* ========================================================================== */

double CHOLMOD(norm_sparse)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to compute the norm of */
    int norm,       /* type of norm: 0: inf. norm, 1: 1-norm */
    /* --------------- */
    cholmod_common *Common
)
{
    double anorm, s ;
    double *Ax, *Az, *W ;
    Int *Ap, *Ai, *Anz ;
    Int i, j, p, pend, nrow, ncol, packed, xtype ;

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (EMPTY) ;
    RETURN_IF_NULL (A, EMPTY) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
    Common->status = CHOLMOD_OK ;
    ncol = A->ncol ;
    nrow = A->nrow ;
    if (norm < 0 || norm > 1)
    {
    ERROR (CHOLMOD_INVALID, "invalid norm") ;
    return (EMPTY) ;
    }
    if (A->stype && nrow != ncol)
    {
    ERROR (CHOLMOD_INVALID, "matrix invalid") ;
    return (EMPTY) ;
    }

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    Ap = A->p ;
    Ai = A->i ;
    Ax = A->x ;
    Az = A->z ;
    Anz = A->nz ;
    packed = A->packed ;
    xtype = A->xtype ;

    /* ---------------------------------------------------------------------- */
    /* allocate workspace, if needed */
    /* ---------------------------------------------------------------------- */

    W = NULL ;
    if (A->stype || norm == 0)
    {
    CHOLMOD(allocate_work) (0, 0, nrow, Common) ;
    W = Common->Xwork ;
    if (Common->status < CHOLMOD_OK)
    {
        /* out of memory */
        return (EMPTY) ;
    }
    DEBUG (for (i = 0 ; i < nrow ; i++) ASSERT (W [i] == 0)) ;
    }

    /* ---------------------------------------------------------------------- */
    /* compute the norm */
    /* ---------------------------------------------------------------------- */

    anorm = 0 ;

    if (A->stype > 0)
    {

    /* ------------------------------------------------------------------ */
    /* A is symmetric with upper triangular part stored */
    /* ------------------------------------------------------------------ */

    /* infinity-norm = 1-norm = max row/col sum */
    for (j = 0 ; j < ncol ; j++)
    {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
        i = Ai [p] ;
        s = abs_value (xtype, Ax, Az, p, Common) ;
        if (i == j)
        {
            W [i] += s ;
        }
        else if (i < j)
        {
            W [i] += s ;
            W [j] += s ;
        }
        }
    }

    }
    else if (A->stype < 0)
    {

    /* ------------------------------------------------------------------ */
    /* A is symmetric with lower triangular part stored */
    /* ------------------------------------------------------------------ */

    /* infinity-norm = 1-norm = max row/col sum */
    for (j = 0 ; j < ncol ; j++)
    {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
        i = Ai [p] ;
        s = abs_value (xtype, Ax, Az, p, Common) ;
        if (i == j)
        {
            W [i] += s ;
        }
        else if (i > j)
        {
            W [i] += s ;
            W [j] += s ;
        }
        }
    }

    }
    else if (norm == 0)
    {

    /* ------------------------------------------------------------------ */
    /* A is unsymmetric, compute the infinity-norm */
    /* ------------------------------------------------------------------ */

    /* infinity-norm = max row sum */
    for (j = 0 ; j < ncol ; j++)
    {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
        W [Ai [p]] += abs_value (xtype, Ax, Az, p, Common) ;
        }
    }

    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* A is unsymmetric, compute the 1-norm */
    /* ------------------------------------------------------------------ */

    /* 1-norm = max column sum */
    for (j = 0 ; j < ncol ; j++)
    {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        if (xtype == CHOLMOD_PATTERN)
        {
        s = pend - p ;
        }
        else
        {
        s = 0 ;
        for ( ; p < pend ; p++)
        {
            s += abs_value (xtype, Ax, Az, p, Common) ;
        }
        }
        if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
        {
        anorm = s ;
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* compute the max row sum */
    /* ---------------------------------------------------------------------- */

    if (A->stype || norm == 0)
    {
    for (i = 0 ; i < nrow ; i++)
    {
        s = W [i] ;
        if ((IS_NAN (s) || s > anorm) && !IS_NAN (anorm))
        {
        anorm = s ;
        }
        W [i] = 0 ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* return result */
    /* ---------------------------------------------------------------------- */

    return (anorm) ;
}
#endif