The R Project SVN R-packages

/* ========================================================================== */
/* === MatrixOps/cholmod_drop =============================================== */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * 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
 * -------------------------------------------------------------------------- */

/* Drop small entries from A, and entries in the ignored part of A if A
 * is symmetric.  None of the matrix operations drop small numerical entries
 * from a matrix, except for this one.  NaN's and Inf's are kept.
 *
 * workspace: none
 *
 * Supports pattern and real matrices, complex and zomplex not supported.
 */

#ifndef NMATRIXOPS

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


/* ========================================================================== */
/* === cholmod_drop ========================================================= */
/* ========================================================================== */

int CHOLMOD(drop)
(
    /* ---- input ---- */
    double tol,     /* keep entries with absolute value > tol */
    /* ---- in/out --- */
    cholmod_sparse *A,  /* matrix to drop entries from */
    /* --------------- */
    cholmod_common *Common
)
{
    double aij ;
    double *Ax ;
    Int *Ap, *Ai, *Anz ;
    Int packed, i, j, nrow, ncol, p, pend, nz, values ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_REAL, FALSE) ;
    Common->status = CHOLMOD_OK ;
    ASSERT (CHOLMOD(dump_sparse) (A, "A predrop", Common) >= 0) ;

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

    Ap = A->p ;
    Ai = A->i ;
    Ax = A->x ;
    Anz = A->nz ;
    packed = A->packed ;
    ncol = A->ncol ;
    nrow = A->nrow ;
    values = (A->xtype != CHOLMOD_PATTERN) ;
    nz = 0 ;

    if (values)
    {

    /* ------------------------------------------------------------------ */
    /* drop small numerical entries from A, and entries in ignored part */
    /* ------------------------------------------------------------------ */

    if (A->stype > 0)
    {

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

        for (j = 0 ; j < ncol ; j++)
        {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        Ap [j] = nz ;
        for ( ; p < pend ; p++)
        {
            i = Ai [p] ;
            aij = Ax [p] ;
            if (i <= j && (fabs (aij) > tol || IS_NAN (aij)))
            {
            Ai [nz] = i ;
            Ax [nz] = aij ;
            nz++ ;
            }
        }
        }

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

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

        for (j = 0 ; j < ncol ; j++)
        {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        Ap [j] = nz ;
        for ( ; p < pend ; p++)
        {
            i = Ai [p] ;
            aij = Ax [p] ;
            if (i >= j && (fabs (aij) > tol || IS_NAN (aij)))
            {
            Ai [nz] = i ;
            Ax [nz] = aij ;
            nz++ ;
            }
        }
        }
    }
    else
    {

        /* -------------------------------------------------------------- */
        /* both parts of A present, just drop small entries */
        /* -------------------------------------------------------------- */

        for (j = 0 ; j < ncol ; j++)
        {
        p = Ap [j] ;
        pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
        Ap [j] = nz ;
        for ( ; p < pend ; p++)
        {
            i = Ai [p] ;
            aij = Ax [p] ;
            if (fabs (aij) > tol || IS_NAN (aij))
            {
            Ai [nz] = i ;
            Ax [nz] = aij ;
            nz++ ;
            }
        }
        }
    }
    Ap [ncol] = nz ;

    /* reduce A->i and A->x in size */
    ASSERT (MAX (1,nz) <= A->nzmax) ;
    CHOLMOD(reallocate_sparse) (nz, A, Common) ;
    ASSERT (Common->status >= CHOLMOD_OK) ;

    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* consider only the pattern of A */
    /* ------------------------------------------------------------------ */

    /* Note that cholmod_band_inplace calls cholmod_reallocate_sparse */
    if (A->stype > 0)
    {
        CHOLMOD(band_inplace) (0, ncol, 0, A, Common) ;
    }
    else if (A->stype < 0)
    {
        CHOLMOD(band_inplace) (-nrow, 0, 0, A, Common) ;
    }
    }

    ASSERT (CHOLMOD(dump_sparse) (A, "A dropped", Common) >= 0) ;
    return (TRUE) ;
}
#endif