The R Project SVN R-packages

/* ========================================================================== */
/* === Core/cholmod_transpose =============================================== */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Core Module.  Version 0.6.  Copyright (C) 2005, 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
 * -------------------------------------------------------------------------- */

/* Core utility routines for the cholmod_sparse object to
 * compute the transpose or permuted transpose of a matrix:
 *
 * Primary routines:
 * -----------------
 * cholmod_transpose        transpose sparse matrix
 * cholmod_ptranspose       transpose and permute sparse matrix
 * cholmod_sort         sort row indices in each column of sparse matrix
 *
 * Secondary routines:
 * -------------------
 * cholmod_transpose_unsym  transpose unsymmetric sparse matrix
 * cholmod_transpose_sym    transpose symmetric sparse matrix
 *
 * All xtypes (pattern, real, complex, and zomplex) are supported.
 *
 * ---------------------------------------
 * Unsymmetric case: A->stype is zero.
 * ---------------------------------------
 *
 * Computes F = A', F = A (:,f)' or F = A (p,f)', except that the indexing by
 * f does not work the same as the MATLAB notation (see below).  A->stype
 * is zero, which denotes that both the upper and lower triangular parts of
 * A are present (and used).  A may in fact be symmetric in pattern and/or
 * value; A->stype just denotes which part of A are stored.  A may be
 * rectangular.
 *
 * p is a permutation of 0:m-1, and f is a subset of 0:n-1, where A is m-by-n.
 * There can be no duplicate entries in p or f.
 *
 * The set f is held in fset and fsize.
 *  fset = NULL means ":" in MATLAB. fset is ignored.
 *  fset != NULL means f = fset [0..fset-1].
 *  fset != NULL and fsize = 0 means f is the empty set.
 *
 * Columns not in the set f are considered to be zero.  That is,
 * if A is 5-by-10 then F = A (:,[3 4])' is not 2-by-5, but 10-by-5, and rows
 * 3 and 4 of F are equal to columns 3 and 4 of A (the other rows of F are
 * zero).  More precisely, in MATLAB notation:
 *
 *  [m n] = size (A) ;
 *  F = A ;
 *  notf = ones (1,n) ;
 *  notf (f) = 0 ;
 *  F (:, find (notf)) = 0
 *  F = F'
 *
 * If you want the MATLAB equivalent F=A(p,f) operation, use cholmod_submatrix
 * instead (which does not compute the transpose).
 *
 * F->nzmax must be large enough to hold the matrix F.  It is not modified.
 * If F->nz is present then F->nz [j] = # of entries in column j of F.
 *
 * A can be sorted or unsorted, with packed or unpacked columns.
 *
 * If f is present and not sorted in ascending order, then F is unsorted
 * (that is, it may contain columns whose row indices do not appear in
 * ascending order).  Otherwise, F is sorted (the row indices in each
 * column of F appear in strictly ascending order).
 *
 * F is returned in packed or unpacked form, depending on F->packed on input.
 * If F->packed is false, then F is returned in unpacked form (F->nz must be
 * present).  Each row i of F is large enough to hold all the entries in row i
 * of A, even if f is provided.  That is, F->i and
 * F->x [F->p [i] .. F->p [i] + F->nz [i] - 1] contain all entries in A (i,f),
 * but F->p [i+1] - F->p [i] is equal to the number of nonzeros in A (i,:),
 * not just A (i,f).
 *
 * The cholmod_transpose_unsym routine is the only operation in CHOLMOD that
 * can produce an unpacked matrix.
 *
 * ---------------------------------------
 * Symmetric case: A->stype is nonzero.
 * ---------------------------------------
 *
 * Computes F = A' or F = A(p,p)', the transpose or permuted transpose, where
 * A->stype is nonzero.
 *
 * If A->stype > 0, then A is a symmetric matrix where just the upper part
 * of the matrix is stored.  Entries in the lower triangular part may be
 * present, but are ignored.  A must be square.  If F=A', then F is returned
 * sorted; otherwise F is unsorted for the F=A(p,p)' case.
 *
 * There can be no duplicate entries in p.
 * The fset and fsize parameters are not used.
 *
 * Three kinds of transposes are available, depending on the "values" parameter:
 * 0: do not transpose the numerical values; create a CHOLMOD_PATTERN matrix
 * 1: array transpose
 * 2: complex conjugate transpose (same as 2 if input is real or pattern)
 *
 * -----------------------------------------------------------------------------
 *
 * For cholmod_transpose_unsym and cholmod_transpose_sym, the output matrix
 * F must already be pre-allocated by the caller, with the correct dimensions.
 * If F is not valid or has the wrong dimensions, it is not modified.
 * Otherwise, if F is too small, the transpose is not computed; the contents
 * of F->p contain the column pointers of the resulting matrix, where
 * F->p [F->ncol] > F->nzmax.  In this case, the remaining contents of F are
 * not modified.  F can still be properly free'd with cholmod_free_sparse.
 */

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


/* ========================================================================== */
/* === TEMPLATE ============================================================= */
/* ========================================================================== */

#define PATTERN
#include "t_cholmod_transpose.c"
#define REAL
#include "t_cholmod_transpose.c"
#define COMPLEX
#include "t_cholmod_transpose.c"
#define COMPLEX
#define NCONJUGATE
#include "t_cholmod_transpose.c"
#define ZOMPLEX
#include "t_cholmod_transpose.c"
#define ZOMPLEX
#define NCONJUGATE
#include "t_cholmod_transpose.c"


/* ========================================================================== */
/* === cholmod_transpose_unsym ============================================== */
/* ========================================================================== */

/* Compute F = A', A (:,f)', or A (p,f)', where A is unsymmetric and F is
 * already allocated.  See cholmod_transpose for a simpler routine.
 *
 * workspace:
 * Iwork (MAX (nrow,ncol)) if fset is present
 * Iwork (nrow) if fset is NULL
 *
 * The xtype of A and F must match, unless values is zero or F->xtype is
 * CHOLMOD_PATTERN (in which case only the pattern of A is transpose into F).
 */

int CHOLMOD(transpose_unsym)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to transpose */
    int values,     /* 2: complex conj. transpose, 1: array transpose,
               0: do not transpose the numerical values */
    Int *Perm,      /* size nrow, if present (can be NULL) */
    Int *fset,      /* subset of 0:(A->ncol)-1 */
    size_t fsize,   /* size of fset */
    /* ---- output --- */
    cholmod_sparse *F,  /* F = A', A(:,f)', or A(p,f)' */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Fp, *Fnz, *Ap, *Ai, *Anz, *Wi ;
    Int nrow, ncol, ok = FALSE, permute, use_fset, Apacked, Fpacked, p, pend,
    i, j, k, Fsorted, nf, jj, jlast ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_NULL (F, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    if (A->nrow != F->ncol || A->ncol != F->nrow)
    {
    ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
    return (FALSE) ;
    }
    Common->status = CHOLMOD_OK ;

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

    nf = fsize ;
    use_fset = (fset != NULL) ;
    nrow = A->nrow ;
    ncol = A->ncol ;

    Ap = A->p ;     /* size A->ncol+1, column pointers of A */
    Ai = A->i ;     /* size nz = Ap [A->ncol], row indices of A */
    Anz = A->nz ;
    Apacked = A->packed ;
    ASSERT (IMPLIES (!Apacked, Anz != NULL)) ;

    permute = (Perm != NULL) ;

    Fp = F->p ;     /* size A->nrow+1, row pointers of F */
    Fnz = F->nz ;
    Fpacked = F->packed ;
    ASSERT (IMPLIES (!Fpacked, Fnz != NULL)) ;

    nf = (use_fset) ? nf : ncol ;

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

    CHOLMOD(allocate_work) (0, nrow + (fset != NULL) ? ncol : 0, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;    /* out of memory */
    }

    Wi = Common->Iwork ;    /* size nrow (i/l/l) */

    /* ---------------------------------------------------------------------- */
    /* check Perm and fset */
    /* ---------------------------------------------------------------------- */

    if (permute)
    {
    for (i = 0 ; i < nrow ; i++)
    {
        Wi [i] = 1 ;
    }
    for (k = 0 ; k < nrow ; k++)
    {
        i = Perm [k] ;
        if (i < 0 || i > nrow || Wi [i] == 0)
        {
        ERROR (CHOLMOD_INVALID, "invalid permutation") ;
        return (FALSE) ;
        }
        Wi [i] = 0 ;
    }
    }

    if (use_fset)
    {
    for (j = 0 ; j < ncol ; j++)
    {
        Wi [j] = 1 ;
    }
    for (k = 0 ; k < nf ; k++)
    {
        j = fset [k] ;
        if (j < 0 || j > ncol || Wi [j] == 0)
        {
        ERROR (CHOLMOD_INVALID, "invalid fset") ;
        return (FALSE) ;
        }
        Wi [j] = 0 ;
    }
    }

    /* Perm and fset are now valid */
    ASSERT (CHOLMOD(dump_perm) (Perm, nrow, nrow, "Perm", Common)) ;
    ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;

    /* ---------------------------------------------------------------------- */
    /* count the entries in each row of A or A(:,f) */
    /* ---------------------------------------------------------------------- */

    for (i = 0 ; i < nrow ; i++)
    {
    Wi [i] = 0 ;
    }

    jlast = EMPTY ;
    Fsorted = TRUE ;

    if (use_fset)
    {
    /* count entries in each row of A(:,f) */
    for (jj = 0 ; jj < nf ; jj++)
    {
        j = fset [jj] ;
        if (j <= jlast)
        {
        Fsorted = FALSE ;
        }
        p = Ap [j] ;
        pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
        Wi [Ai [p]]++ ;
        }
        jlast = j ;
    }

    /* save the nz counts if F is unpacked, and recount all of A */
    if (!Fpacked)
    {
        if (permute)
        {
        for (i = 0 ; i < nrow ; i++)
        {
            Fnz [i] = Wi [Perm [i]] ;
        }
        }
        else
        {
        for (i = 0 ; i < nrow ; i++)
        {
            Fnz [i] = Wi [i] ;
        }
        }
        for (i = 0 ; i < nrow ; i++)
        {
        Wi [i] = 0 ;
        }

        /* count entries in each row of A */
        for (j = 0 ; j < ncol ; j++)
        {
        p = Ap [j] ;
        pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
            Wi [Ai [p]]++ ;
        }
        }
    }

    }
    else
    {

    /* count entries in each row of A */
    for (j = 0 ; j < ncol ; j++)
    {
        p = Ap [j] ;
        pend = (Apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
        for ( ; p < pend ; p++)
        {
        Wi [Ai [p]]++ ;
        }
    }

    /* save the nz counts if F is unpacked */
    if (!Fpacked)
    {
        if (permute)
        {
        for (i = 0 ; i < nrow ; i++)
        {
            Fnz [i] = Wi [Perm [i]] ;
        }
        }
        else
        {
        for (i = 0 ; i < nrow ; i++)
        {
            Fnz [i] = Wi [i] ;
        }
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* compute the row pointers */
    /* ---------------------------------------------------------------------- */

    p = 0 ;
    if (permute)
    {
    for (i = 0 ; i < nrow ; i++)
    {
        Fp [i] = p ;
        p += Wi [Perm [i]] ;
    }
    for (i = 0 ; i < nrow ; i++)
    {
        Wi [Perm [i]] = Fp [i] ;
    }
    }
    else
    {
    for (i = 0 ; i < nrow ; i++)
    {
        Fp [i] = p ;
        p += Wi [i] ;
    }
    for (i = 0 ; i < nrow ; i++)
    {
        Wi [i] = Fp [i] ;
    }
    }
    Fp [nrow] = p ;

    if (p > (Int) (F->nzmax))
    {
    ERROR (CHOLMOD_INVALID, "F is too small") ;
    return (FALSE) ;
    }

    /* ---------------------------------------------------------------------- */
    /* transpose matrix, using template routine */
    /* ---------------------------------------------------------------------- */

    if (values == 0 || F->xtype == CHOLMOD_PATTERN)
    {
    ok = p_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    else if (F->xtype == CHOLMOD_REAL)
    {
    ok = r_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    else if (F->xtype == CHOLMOD_COMPLEX)
    {
    if (values == 1)
    {
        /* array transpose */
        ok = ct_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    else
    {
        /* complex conjugate transpose */
        ok = c_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    }
    else if (F->xtype == CHOLMOD_ZOMPLEX)
    {
    if (values == 1)
    {
        /* array transpose */
        ok = zt_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    else
    {
        /* complex conjugate transpose */
        ok = z_cholmod_transpose_unsym (A, Perm, fset, nf, F, Common) ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* finalize result F */
    /* ---------------------------------------------------------------------- */

    if (ok)
    {
    F->sorted = Fsorted ;
    }
    ASSERT (CHOLMOD(dump_sparse) (F, "output F unsym", Common) >= 0) ;
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_transpose_sym ================================================ */
/* ========================================================================== */

/* Compute F = A' or A (p,p)', where A is symmetric and F is already allocated.
 * See cholmod_transpose for a simpler routine.
 *
 * workspace:  Iwork (nrow) if Perm NULL, Iwork (2*nrow) if Perm non-NULL.
 */

int CHOLMOD(transpose_sym)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to transpose */
    int values,     /* 2: complex conj. transpose, 1: array transpose,
               0: do not transpose the numerical values */
    Int *Perm,      /* size nrow, if present (can be NULL) */
    /* ---- output --- */
    cholmod_sparse *F,  /* F = A' or A(p,p)' */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Ap, *Anz, *Ai, *Fp, *Wi, *Pinv, *Iwork ;
    Int p, pend, packed, upper, permute, jold, n, i, j, k, iold ;
    Int ok = FALSE ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_NULL (F, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    RETURN_IF_XTYPE_INVALID (F, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    if (A->nrow != A->ncol || A->stype == 0)
    {
    /* this routine handles square symmetric matrices only */
    ERROR (CHOLMOD_INVALID, "matrix must be symmetric") ;
    return (FALSE) ;
    }
    if (A->nrow != F->ncol || A->ncol != F->nrow)
    {
    ERROR (CHOLMOD_INVALID, "F has the wrong dimensions") ;
    return (FALSE) ;
    }
    Common->status = CHOLMOD_OK ;

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

    permute = (Perm != NULL) ;
    n = A->nrow ;
    Ap = A->p ;     /* size A->ncol+1, column pointers of A */
    Ai = A->i ;     /* size nz = Ap [A->ncol], row indices of A */
    Anz = A->nz ;
    packed = A->packed ;
    ASSERT (IMPLIES (!packed, Anz != NULL)) ;
    upper = (A->stype > 0) ;

    Fp = F->p ;     /* size A->nrow+1, row pointers of F */

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

    CHOLMOD(allocate_work) (0, (Perm != NULL) ? 2*n : n, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;    /* out of memory */
    }

    /* ---------------------------------------------------------------------- */
    /* get workspace */
    /* ---------------------------------------------------------------------- */

    Iwork = Common->Iwork ;
    Wi   = Iwork ;      /* size n (i/l/l) */
    Pinv = Iwork + n ;      /* size n (i/i/l) , unused if Perm NULL */

    /* ---------------------------------------------------------------------- */
    /* check Perm and construct inverse permutation */
    /* ---------------------------------------------------------------------- */

    if (permute)
    {
    for (i = 0 ; i < n ; i++)
    {
        Pinv [i] = EMPTY ;
    }
    for (k = 0 ; k < n ; k++)
    {
        i = Perm [k] ;
        if (i < 0 || i > n || Pinv [i] != EMPTY)
        {
        ERROR (CHOLMOD_INVALID, "invalid permutation") ;
        return (FALSE) ;
        }
        Pinv [i] = k ;
    }
    }

    /* Perm is now valid */
    ASSERT (CHOLMOD(dump_perm) (Perm, n, n, "Perm", Common)) ;

    /* ---------------------------------------------------------------------- */
    /* count the entries in each row of F */
    /* ---------------------------------------------------------------------- */

    for (i = 0 ; i < n ; i++)
    {
    Wi [i] = 0 ;
    }

    if (packed)
    {
    if (permute)
    {
        if (upper)
        {
        /* packed, permuted, upper */
        for (j = 0 ; j < n ; j++)
        {
            jold = Perm [j] ;
            pend = Ap [jold+1] ;
            for (p = Ap [jold] ; p < pend ; p++)
            {
            iold = Ai [p] ;
            if (iold <= jold)
            {
                i = Pinv [iold] ;
                Wi [MIN (i, j)]++ ;
            }
            }
        }
        }
        else
        {
        /* packed, permuted, lower */
        for (j = 0 ; j < n ; j++)
        {
            jold = Perm [j] ;
            pend = Ap [jold+1] ;
            for (p = Ap [jold] ; p < pend ; p++)
            {
            iold = Ai [p] ;
            if (iold >= jold)
            {
                i = Pinv [iold] ;
                Wi [MAX (i, j)]++ ;
            }
            }
        }
        }
    }
    else
    {
        if (upper)
        {
        /* packed, unpermuted, upper */
        for (j = 0 ; j < n ; j++)
        {
            pend = Ap [j+1] ;
            for (p = Ap [j] ; p < pend ; p++)
            {
            i = Ai [p] ;
            if (i <= j)
            {
                Wi [i]++ ;
            }
            }
        }
        }
        else
        {
        /* packed, unpermuted, lower */
        for (j = 0 ; j < n ; j++)
        {
            pend = Ap [j+1] ;
            for (p = Ap [j] ; p < pend ; p++)
            {
            i = Ai [p] ;
            if (i >= j)
            {
                Wi [i]++ ;
            }
            }
        }
        }
    }
    }
    else
    {
    if (permute)
    {
        if (upper)
        {
        /* unpacked, permuted, upper */
        for (j = 0 ; j < n ; j++)
        {
            jold = Perm [j] ;
            p = Ap [jold] ;
            pend = p + Anz [jold] ;
            for ( ; p < pend ; p++)
            {
            iold = Ai [p] ;
            if (iold <= jold)
            {
                i = Pinv [iold] ;
                Wi [MIN (i, j)]++ ;
            }
            }
        }
        }
        else
        {
        /* unpacked, permuted, lower */
        for (j = 0 ; j < n ; j++)
        {
            jold = Perm [j] ;
            p = Ap [jold] ;
            pend = p + Anz [jold] ;
            for ( ; p < pend ; p++)
            {
            iold = Ai [p] ;
            if (iold >= jold)
            {
                i = Pinv [iold] ;
                Wi [MAX (i, j)]++ ;
            }
            }
        }
        }
    }
    else
    {
        if (upper)
        {
        /* unpacked, unpermuted, upper */
        for (j = 0 ; j < n ; j++)
        {
            p = Ap [j] ;
            pend = p + Anz [j] ;
            for ( ; p < pend ; p++)
            {
            i = Ai [p] ;
            if (i <= j)
            {
                Wi [i]++ ;
            }
            }
        }
        }
        else
        {
        /* unpacked, unpermuted, lower */
        for (j = 0 ; j < n ; j++)
        {
            p = Ap [j] ;
            pend = p + Anz [j] ;
            for ( ; p < pend ; p++)
            {
            i = Ai [p] ;
            if (i >= j)
            {
                Wi [i]++ ;
            }
            }
        }
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* compute the row pointers */
    /* ---------------------------------------------------------------------- */

    p = 0 ;
    for (i = 0 ; i < n ; i++)
    {
    Fp [i] = p ;
    p += Wi [i] ;
    }
    Fp [n] = p ;
    for (i = 0 ; i < n ; i++)
    {
    Wi [i] = Fp [i] ;
    }

    if (p > (Int) (F->nzmax))
    {
    ERROR (CHOLMOD_INVALID, "F is too small") ;
    return (FALSE) ;
    }

    /* ---------------------------------------------------------------------- */
    /* transpose matrix, using template routine */
    /* ---------------------------------------------------------------------- */

    if (values == 0 || F->xtype == CHOLMOD_PATTERN)
    {
    PRINT2 (("\n:::: p_transpose_sym Perm %p\n", Perm)) ;
    ok = p_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    else if (F->xtype == CHOLMOD_REAL)
    {
    PRINT2 (("\n:::: r_transpose_sym Perm %p\n", Perm)) ;
    ok = r_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    else if (F->xtype == CHOLMOD_COMPLEX)
    {
    if (values == 1)
    {
        /* array transpose */
        PRINT2 (("\n:::: ct_transpose_sym Perm %p\n", Perm)) ;
        ok = ct_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    else
    {
        /* complex conjugate transpose */
        PRINT2 (("\n:::: c_transpose_sym Perm %p\n", Perm)) ;
        ok = c_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    }
    else if (F->xtype == CHOLMOD_ZOMPLEX)
    {
    if (values == 1)
    {
        /* array transpose */
        PRINT2 (("\n:::: zt_transpose_sym Perm %p\n", Perm)) ;
        ok = zt_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    else
    {
        /* complex conjugate transpose */
        PRINT2 (("\n:::: z_transpose_sym Perm %p\n", Perm)) ;
        ok = z_cholmod_transpose_sym (A, Perm, F, Common) ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* finalize result F */
    /* ---------------------------------------------------------------------- */

    /* F is sorted if there is no permutation vector */
    if (ok)
    {
    F->sorted = !permute ;
    F->packed = TRUE ;
    F->stype = - SIGN (A->stype) ;  /* flip the stype */
    ASSERT (CHOLMOD(dump_sparse) (F, "output F sym", Common) >= 0) ;
    }
    return (ok) ;
}


/* ========================================================================== */
/* === cholmod_transpose ==================================================== */
/* ========================================================================== */

/* Returns A'.  See also cholmod_ptranspose below. */

cholmod_sparse *CHOLMOD(transpose)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to transpose */
    int values,     /* 2: complex conj. transpose, 1: array transpose,
               0: do not transpose the numerical values
               (returns its result as CHOLMOD_PATTERN) */
    /* --------------- */
    cholmod_common *Common
)
{
    return (CHOLMOD(ptranspose) (A, values, NULL, NULL, 0, Common)) ;
}


/* ========================================================================== */
/* === cholmod_ptranspose =================================================== */
/* ========================================================================== */

/* Return A' or A(p,p)' if A is symmetric.  Return A', A(:,f)', or A(p,f)' if
 * A is unsymmetric.
 *
 * workspace:
 * Iwork (MAX (nrow,ncol)) if unsymmetric and fset is non-NULL
 * Iwork (nrow) if unsymmetric and fset is NULL
 * Iwork (2*nrow) if symmetric and Perm is non-NULL.
 * Iwork (nrow) if symmetric and Perm is NULL.
 *
 * A simple worst-case upper bound on the workspace is nrow+ncol.
 */

cholmod_sparse *CHOLMOD(ptranspose)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to transpose */
    int values,     /* 2: complex conj. transpose, 1: array transpose,
               0: do not transpose the numerical values */
    Int *Perm,      /* if non-NULL, F = A(p,f) or A(p,p) */
    Int *fset,      /* subset of 0:(A->ncol)-1 */
    size_t fsize,   /* size of fset */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Ap, *Anz ;
    cholmod_sparse *F ;
    Int nrow, ncol, use_fset, j, jj, fnz, packed, stype, nf, ineed, ok,
    xtype ;

    nf = fsize ;

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

    RETURN_IF_NULL_COMMON (NULL) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
    stype = A->stype ;
    Common->status = CHOLMOD_OK ;

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

    nrow = A->nrow ;
    ncol = A->ncol ;

    if (stype != 0)
    {
    use_fset = FALSE ;
    if (Perm != NULL)
    {
        ineed = 2*nrow ;
    }
    else
    {
        ineed = nrow ;
    }
    }
    else
    {
    use_fset = (fset != NULL) ;
    if (use_fset)
    {
        ineed = MAX (nrow, ncol) ;
    }
    else
    {
        ineed = nrow ;
    }
    }

    CHOLMOD(allocate_work) (0, ineed, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (NULL) ;     /* out of memory */
    }

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

    Ap = A->p ;
    Anz = A->nz ;
    packed = A->packed ;
    ASSERT (IMPLIES (!packed, Anz != NULL)) ;
    xtype = values ? A->xtype : CHOLMOD_PATTERN ;

    /* ---------------------------------------------------------------------- */
    /* allocate F */
    /* ---------------------------------------------------------------------- */

    /* determine # of nonzeros in F */
    if (stype != 0)
    {
    /* F=A' or F=A(p,p)', fset is ignored */
    fnz = CHOLMOD(nnz) (A, Common) ;
    }
    else
    {
    nf = (use_fset) ? nf : ncol ;
    if (use_fset)
    {
        fnz = 0 ;
        /* F=A(:,f)' or F=A(p,f)' */
        for (jj = 0 ; jj < nf ; jj++)
        {
        /* The fset is not yet checked; it will be thoroughly checked
         * in cholmod_transpose_unsym.  For now, just make sure we don't
         * access Ap and Anz out of bounds. */
        j = fset [jj] ;
        if (j >= 0 && j < ncol)
        {
            fnz += packed ? (Ap [j+1] - Ap [j]) : MAX (0, Anz [j]) ;
        }
        }
    }
    else
    {
        /* F=A' or F=A(p,:)' */
        fnz = CHOLMOD(nnz) (A, Common) ;
    }
    }

    /* F is ncol-by-nrow, fnz nonzeros, sorted unless f is present and unsorted,
     * packed, of opposite stype as A, and with/without numerical values */
    F = CHOLMOD(allocate_sparse) (ncol, nrow, fnz, TRUE, TRUE, -SIGN(stype),
        xtype, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (NULL) ;     /* out of memory */
    }

    /* ---------------------------------------------------------------------- */
    /* transpose and optionally permute the matrix A */
    /* ---------------------------------------------------------------------- */

    if (stype != 0)
    {
    /* F = A (p,p)', using upper or lower triangular part of A only */
    ok = CHOLMOD(transpose_sym) (A, values, Perm, F, Common) ;
    }
    else
    {
    /* F = A (p,f)' */
    ok = CHOLMOD(transpose_unsym) (A, values, Perm, fset, nf, F, Common) ;
    }

    /* ---------------------------------------------------------------------- */
    /* return the matrix F, or NULL if an error occured */
    /* ---------------------------------------------------------------------- */

    if (!ok)
    {
    CHOLMOD(free_sparse) (&F, Common) ;
    }
    return (F) ;
}


/* ========================================================================== */
/* === cholmod_sort ========================================================= */
/* ========================================================================== */

/* Sort the columns of A, in place.  Returns A in packed form, even if it
 * starts as unpacked.  Removes entries in the ignored part of a symmetric
 * matrix.
 *
 * workspace: Iwork (max (nrow,ncol)).  Allocates additional workspace for a
 * temporary copy of A'.
 */

int CHOLMOD(sort)
(
    /* ---- in/out --- */
    cholmod_sparse *A,  /* matrix to sort */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Ap ;
    cholmod_sparse *F ;
    Int anz, ncol, nrow, stype ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    Common->status = CHOLMOD_OK ;
    nrow = A->nrow ;
    if (nrow <= 1)
    {
    /* a 1-by-n sparse matrix must be sorted */
    A->sorted = TRUE ;
    return (TRUE) ;
    }

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

    ncol = A->ncol ;
    CHOLMOD(allocate_work) (0, MAX (nrow, ncol), 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;    /* out of memory */
    }

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

    anz = CHOLMOD(nnz) (A, Common) ;
    stype = A->stype ;

    /* ---------------------------------------------------------------------- */
    /* sort the columns of the matrix */
    /* ---------------------------------------------------------------------- */

    /* allocate workspace for transpose: ncol-by-nrow, same # of nonzeros as A,
     * sorted, packed, same stype as A, and of the same numeric type as A. */
    F = CHOLMOD(allocate_sparse) (ncol, nrow, anz, TRUE, TRUE, stype,
        A->xtype, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;    /* out of memory */
    }

    if (stype != 0)
    {
    /* F = A', upper or lower triangular part only */
    CHOLMOD(transpose_sym) (A, 1, NULL, F, Common) ;
    A->packed = TRUE ;
    /* A = F' */
    CHOLMOD(transpose_sym) (F, 1, NULL, A, Common) ;
    }
    else
    {
    /* F = A' */
    CHOLMOD(transpose_unsym) (A, 1, NULL, NULL, 0, F, Common) ;
    A->packed = TRUE ;
    /* A = F' */
    CHOLMOD(transpose_unsym) (F, 1, NULL, NULL, 0, A, Common) ;
    }

    ASSERT (A->sorted && A->packed) ;
    ASSERT (CHOLMOD(dump_sparse) (A, "Asorted", Common) >= 0) ;

    /* ---------------------------------------------------------------------- */
    /* reduce A in size, if needed.  This must succeed. */
    /* ---------------------------------------------------------------------- */

    Ap = A->p ;
    anz = Ap [ncol] ;
    ASSERT ((size_t) anz <= A->nzmax) ;
    CHOLMOD(reallocate_sparse) (anz, A, Common) ;
    ASSERT (Common->status >= CHOLMOD_OK) ;

    /* ---------------------------------------------------------------------- */
    /* free workspace */
    /* ---------------------------------------------------------------------- */

    CHOLMOD(free_sparse) (&F, Common) ;
    return (TRUE) ;
}