The R Project SVN R-packages

/* ========================================================================== */
/* === MatrixOps/cholmod_symmetry =========================================== */
/* ========================================================================== */

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

/* Determines if a sparse matrix is rectangular, unsymmetric, symmetric,
 * skew-symmetric, or Hermitian.  It does so by looking at its numerical values
 * of both upper and lower triangular parts of a CHOLMOD "unsymmetric"
 * matrix, where A->stype == 0.  The transpose of A is NOT constructed.
 *
 * If not unsymmetric, it also determines if the matrix has a diagonal whose
 * entries are all real and positive (and thus a candidate for sparse Cholesky
 * if A->stype is changed to a nonzero value).
 *
 * Note that a Matrix Market "general" matrix is either rectangular or
 * unsymmetric.
 *
 * The row indices in the column of each matrix MUST be sorted for this function
 * to work properly (A->sorted must be TRUE).  This routine returns EMPTY if
 * A->stype is not zero, or if A->sorted is FALSE.  The exception to this rule
 * is if A is rectangular.
 *
 * If option == 0, then this routine returns immediately when it finds a
 * non-positive diagonal entry (or one with nonzero imaginary part).   If the
 * matrix is not a candidate for sparse Cholesky, it returns the value
 * CHOLMOD_MM_UNSYMMETRIC, even if the matrix might in fact be symmetric or
 * Hermitian.
 *
 * This routine is useful inside the MATLAB backslash, which must look at an
 * arbitrary matrix (A->stype == 0) and determine if it is a candidate for
 * sparse Cholesky.  In that case, option should be 0.
 *
 * This routine is also useful when writing a MATLAB matrix to a file in
 * Rutherford/Boeing or Matrix Market format.  Those formats require a
 * determination as to the symmetry of the matrix, and thus this routine should
 * not return upon encountering the first non-positive diagonal.  In this case,
 * option should be 1.
 *
 * If option is 2, this function can be used to compute the numerical and
 * pattern symmetry, where 0 is a completely unsymmetric matrix, and 1 is a
 * perfectly symmetric matrix.  This option is used when computing the following
 * statistics for the matrices in the UF Sparse Matrix Collection.
 *
 *  numerical symmetry: number of matched offdiagonal nonzeros over
 *  the total number of offdiagonal entries.  A real entry A(i,j), i ~= j,
 *  is matched if A (j,i) == A (i,j), but this is only counted if both
 *  A(j,i) and A(i,j) are nonzero.  This does not depend on Z.
 *  (If A is complex, then the above test is modified; A (i,j) is matched
 *  if conj (A (j,i)) == A (i,j)).
 *
 *  Then numeric symmetry = xmatched / nzoffdiag, or 1 if nzoffdiag = 0.
 *  
 *  pattern symmetry: number of matched offdiagonal entries over the
 *  total number of offdiagonal entries.  An entry A(i,j), i ~= j, is
 *  matched if A (j,i) is also an entry.
 *
 *  Then pattern symmetry = pmatched / nzoffdiag, or 1 if nzoffdiag = 0.
 *  
 * The symmetry of a matrix with no offdiagonal entries is equal to 1.
 *
 * A workspace of size ncol integers is allocated; EMPTY is returned if this
 * allocation fails.
 *
 * Summary of return values:
 *
 *  EMPTY (-1)              out of memory, stype not zero, A not sorted
 *  CHOLMOD_MM_RECTANGULAR 1        A is rectangular
 *  CHOLMOD_MM_UNSYMMETRIC 2        A is unsymmetric
 *  CHOLMOD_MM_SYMMETRIC 3      A is symmetric, but with non-pos. diagonal
 *  CHOLMOD_MM_HERMITIAN 4      A is Hermitian, but with non-pos. diagonal
 *  CHOLMOD_MM_SKEW_SYMMETRIC 5     A is skew symmetric
 *  CHOLMOD_MM_SYMMETRIC_POSDIAG 6  A is symmetric with positive diagonal
 *  CHOLMOD_MM_HERMITIAN_POSDIAG 7  A is Hermitian with positive diagonal
 *
 * See also the spsym mexFunction, which is a MATLAB interface for this code.
 *
 * If the matrix is a candidate for sparse Cholesky, it will return a result
 * CHOLMOD_MM_SYMMETRIC_POSDIAG if real, or CHOLMOD_MM_HERMITIAN_POSDIAG if
 * complex.  Otherwise, it will return a value less than this.  This is true
 * regardless of the value of the option parameter.
 */

#ifndef NMATRIXOPS

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


/* ========================================================================== */
/* === get_value ============================================================ */
/* ========================================================================== */

/* Get the pth value in the matrix. */

static void get_value
(
    double *Ax,     /* real values, or real/imag. for CHOLMOD_COMPLEX type */
    double *Az,     /* imaginary values for CHOLMOD_ZOMPLEX type */
    Int p,      /* get the pth entry */
    Int xtype,      /* A->xtype: pattern, real, complex, or zomplex */
    double *x,      /* the real part */
    double *z       /* the imaginary part */
)
{
    switch (xtype)
    {
    case CHOLMOD_PATTERN:
        *x = 1 ;
        *z = 0 ;
        break ;

    case CHOLMOD_REAL:
        *x = Ax [p] ;
        *z = 0 ;
        break ;

    case CHOLMOD_COMPLEX:
        *x = Ax [2*p] ;
        *z = Ax [2*p+1] ;
        break ;

    case CHOLMOD_ZOMPLEX:
        *x = Ax [p] ;
        *z = Az [p] ;
        break ;
    }
}


/* ========================================================================== */
/* === cholmod_symmetry ===================================================== */
/* ========================================================================== */

/* Determine the symmetry of a matrix, and check its diagonal.
 *
 * option 0:  Do not count # of matched pairs.  Quick return if the
 *        the matrix has a zero, negative, or imaginary diagonal entry.
 *
 * option 1:  Do not count # of matched pairs.  Do not return quickly if
 *        the matrix has a zero, negative, or imaginary diagonal entry.
 *  The result 1 to 7 is accurately computed:
 *
 *  EMPTY (-1)      out of memory, stype not zero, A not sorted
 *  CHOLMOD_MM_RECTANGULAR 1    A is rectangular
 *  CHOLMOD_MM_UNSYMMETRIC 2    A is unsymmetric
 *  CHOLMOD_MM_SYMMETRIC 3      A is symmetric, with non-pos. diagonal
 *  CHOLMOD_MM_HERMITIAN 4      A is Hermitian, with non-pos. diagonal
 *  CHOLMOD_MM_SKEW_SYMMETRIC 5 A is skew symmetric
 *  CHOLMOD_MM_SYMMETRIC_POSDIAG 6  is symmetric with positive diagonal
 *  CHOLMOD_MM_HERMITIAN_POSDIAG 7  A is Hermitian with positive diagonal
 *
 *  The routine returns as soon as the above is determined (that is, it
 *  can return as soon as it determines the matrix is unsymmetric).
 *
 * option 2:  All of the above, but also compute the number of matched off-
 *  diagonal entries (of two types).  xmatched is the number of 
 *  nonzero entries for which A(i,j) = conj(A(j,i)).  pmatched is
 *  the number of entries (i,j) for which A(i,j) and A(j,i) are both in
 *  the pattern of A (the value doesn't matter).  nzoffdiag is the total
 *  number of off-diagonal entries in the pattern.  nzdiag is the number of
 *  diagonal entries in the pattern.
 *
 * With option 0 or 1, or if the matrix is rectangular, xmatched, pmatched,
 * nzoffdiag, and nzdiag are not computed.
 *
 * Note that a matched pair, A(i,j) and A(j,i) for i != j, is counted twice
 * (once per entry).
 */

int CHOLMOD(symmetry)
(
    /* ---- input ---- */
    cholmod_sparse *A,
    int option,         /* option 0, 1, or 2 (see above) */
    /* ---- output --- */   /* outputs ignored if any are NULL */
    Int *p_xmatched,        /* # of matched numerical entries */
    Int *p_pmatched,        /* # of matched entries in pattern */
    Int *p_nzoffdiag,       /* # of off diagonal entries */
    Int *p_nzdiag,      /* # of diagonal entries */
    /* --------------- */
    cholmod_common *Common
)
{
    double aij_real, aij_imag, aji_real, aji_imag ;
    double *Ax, *Az ;
    Int *Ap, *Ai, *Anz, *munch ;
    Int packed, nrow, ncol, xtype, is_symmetric, is_skew, is_hermitian, posdiag,
    j, p, pend, i, piend, result, xmatched, pmatched, nzdiag, i2, found ;

    /* ---------------------------------------------------------------------- */
    /* 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 ;
    ASSERT (CHOLMOD(dump_sparse) (A, "cholmod_symmetry", Common) >= 0) ;

    if (p_xmatched == NULL || p_pmatched == NULL
    || p_nzoffdiag == NULL || p_nzdiag == NULL)
    {
    /* option 2 is not performed if any output parameter is NULL */
    option = MAX (option, 1) ;
    }

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

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

    /* ---------------------------------------------------------------------- */
    /* check if rectangular, unsorted, or stype is not zero */
    /* ---------------------------------------------------------------------- */

    if (nrow != ncol)
    {
    /* matrix is rectangular */
    return (CHOLMOD_MM_RECTANGULAR) ;
    }

    if (!(A->sorted) || A->stype != 0)
    {
    /* this function cannot determine the type or symmetry */
    return (EMPTY) ;
    }

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    /* this function requires uninitialized Int workspace of size ncol */
    CHOLMOD(allocate_work) (0, ncol, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    /* out of memory */
    return (EMPTY) ;
    }

    munch = Common->Iwork ;     /* the munch array is size ncol */

    /* ---------------------------------------------------------------------- */
    /* determine symmetry of a square matrix */
    /* ---------------------------------------------------------------------- */

    /* a complex or zomplex matrix is Hermitian until proven otherwise */
    is_hermitian = (xtype >= CHOLMOD_COMPLEX) ;

    /* any matrix is symmetric until proven otherwise */
    is_symmetric = TRUE ;

    /* a non-pattern matrix is skew-symmetric until proven otherwise */
    is_skew = (xtype != CHOLMOD_PATTERN) ;

    /* a matrix has positive diagonal entries until proven otherwise */
    posdiag = TRUE ;

    /* munch pointers start at the top of each column */
    for (j = 0 ; j < ncol ; j++)
    {
    munch [j] = Ap [j] ;
    }

    xmatched = 0 ;
    pmatched = 0 ;
    nzdiag = 0 ;

    for (j = 0 ; j < ncol ; j++)    /* examine each column of A */
    {

    /* ------------------------------------------------------------------ */
    /* look at the entire munch column j */
    /* ------------------------------------------------------------------ */

    /* start at the munch point of column j, and go to end of the column */
    p = munch [j] ;
    pend = (packed) ? (Ap [j+1]) : (Ap [j] + Anz [j]) ;

    for ( ; p < pend ; p++)
    {
        /* get the row index of A(i,j) */
        i = Ai [p] ;

        if (i < j)
        {

        /* ---------------------------------------------------------- */
        /* A(i,j) in triu(A), but matching A(j,i) not in tril(A) */
        /* ---------------------------------------------------------- */

        /* entry A(i,j) is unmatched; it appears in the upper triangular
         * part, but not the lower triangular part.  The matrix is
         * unsymmetric. */
        is_hermitian = FALSE ;
        is_symmetric = FALSE ;
        is_skew = FALSE ;

        }
        else if (i == j)
        {

        /* ---------------------------------------------------------- */
        /* the diagonal A(j,j) is present; check its value */
        /* ---------------------------------------------------------- */

        get_value (Ax, Az, p, xtype, &aij_real, &aij_imag) ;
        if (aij_real != 0. || aij_imag != 0.)
        {
            /* diagonal is nonzero; matrix is not skew-symmetric */
            nzdiag++ ;
            is_skew = FALSE ;
        }
        if (aij_real <= 0. || aij_imag != 0.)
        {
            /* diagonal negative or imaginary; not chol candidate */
            posdiag = FALSE ;
        }

        }
        else /* i > j */
        {

        /* ---------------------------------------------------------- */
        /* consider column i, up to and including row j */
        /* ---------------------------------------------------------- */

        /* munch the entry at top of column i up to and incl row j */
        piend = (packed) ? (Ap [i+1]) : (Ap [i] + Anz [i]) ;

        found = FALSE ;

        for ( ; munch [i] < piend ; munch [i]++)
        {

            i2 = Ai [munch [i]] ;

            if (i2 < j)
            {

            /* -------------------------------------------------- */
            /* A(i2,i) in triu(A) but A(i,i2) not in tril(A) */
            /* -------------------------------------------------- */

            /* The matrix is unsymmetric. */
            is_hermitian = FALSE ;
            is_symmetric = FALSE ;
            is_skew = FALSE ;

            }
            else if (i2 == j)
            {

            /* -------------------------------------------------- */
            /* both A(i,j) and A(j,i) exist in the matrix */
            /* -------------------------------------------------- */

            /* this is one more matching entry in the pattern */
            pmatched += 2 ;
            found = TRUE ;

            /* get the value of A(i,j) */
            get_value (Ax, Az, p, xtype, &aij_real, &aij_imag) ;

            /* get the value of A(j,i) */
            get_value (Ax, Az, munch [i],
                xtype, &aji_real, &aji_imag) ;

            /* compare A(i,j) with A(j,i) */
            if (aij_real != aji_real || aij_imag != aji_imag)
            {
                /* the matrix cannot be symmetric */
                is_symmetric = FALSE ;
            }
            if (aij_real != -aji_real || aij_imag != aji_imag)
            {
                /* the matrix cannot be skew-symmetric */
                is_skew = FALSE ;
            }
            if (aij_real != aji_real || aij_imag != -aji_imag)
            {
                /* the matrix cannot be Hermitian */
                is_hermitian = FALSE ;
            }
            else
            {
                /* A(i,j) and A(j,i) are numerically matched */
                xmatched += 2 ;
            }

            }
            else /* i2 > j */
            {

            /* -------------------------------------------------- */
            /* entry A(i2,i) is not munched; consider it later */
            /* -------------------------------------------------- */

            break ;
            }
        }

        if (!found)
        {
            /* A(i,j) in tril(A) but A(j,i) not in triu(A).
             * The matrix is unsymmetric. */
            is_hermitian = FALSE ;
            is_symmetric = FALSE ;
            is_skew = FALSE ;
        }
        }

        if (option < 2 && !(is_symmetric || is_skew || is_hermitian))
        {
        /* matrix is unsymmetric; terminate the test */
        return (CHOLMOD_MM_UNSYMMETRIC) ;
        }
    }

    /* ------------------------------------------------------------------ */
    /* quick return if not Cholesky candidate */
    /* ------------------------------------------------------------------ */

    if (option < 1 && (!posdiag || nzdiag < ncol))
    {
        /* Diagonal entry not present, or present but negative or with
         * nonzero imaginary part.  Quick return for option 0. */
        return (CHOLMOD_MM_UNSYMMETRIC) ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* return the results */
    /* ---------------------------------------------------------------------- */

    if (option >= 2)
    {
    *p_xmatched = xmatched ;
    *p_pmatched = pmatched ;
    *p_nzoffdiag = CHOLMOD(nnz) (A, Common) - nzdiag ;
    *p_nzdiag = nzdiag ;
    }

    result = CHOLMOD_MM_UNSYMMETRIC ;
    if (is_hermitian)
    {
    /* complex Hermitian matrix, with either pos. or non-pos. diagonal */
    result = posdiag ? CHOLMOD_MM_HERMITIAN_POSDIAG : CHOLMOD_MM_HERMITIAN ;
    }
    else if (is_symmetric)
    {
    /* real or complex symmetric matrix, with pos. or non-pos. diagonal */
    result = posdiag ? CHOLMOD_MM_SYMMETRIC_POSDIAG : CHOLMOD_MM_SYMMETRIC ;
    }
    else if (is_skew)
    {
    /* real or complex skew-symmetric matrix */
    result = CHOLMOD_MM_SKEW_SYMMETRIC ;
    }
    return (result) ;
}
#endif