The R Project SVN R-packages

/* ========================================================================== */
/* === MatrixOps/cholmod_ssmult ============================================= */
/* ========================================================================== */

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

/* C = A*B.  Multiply two sparse matrices.
 *
 * A and B can be packed or unpacked, sorted or unsorted, and of any stype.
 * If A or B are symmetric, an internal unsymmetric copy is made first, however.
 * C is computed as if A and B are unsymmetric, and then if the stype input
 * parameter requests a symmetric form (upper or lower) the matrix is converted
 * into that form.
 *
 * C is returned as packed, and either unsorted or sorted, depending on the
 * "sorted" input parameter.  If C is returned sorted, then either C = (B'*A')'
 * or C = (A*B)'' is computed, depending on the number of nonzeros in A, B, and
 * C.
 *
 * workspace:
 *  if C unsorted: Flag (A->nrow), W (A->nrow) if values
 *  if C sorted:   Flag (B->ncol), W (B->ncol) if values
 *  Iwork (max (A->ncol, A->nrow, B->nrow, B->ncol))
 *  allocates temporary copies for A, B, and C, if required.
 *
 * Only pattern and real matrices are supported.  Complex and zomplex matrices
 * are supported only when the numerical values are not computed ("values"
 * is FALSE).
 */

#ifndef NMATRIXOPS

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


/* ========================================================================== */
/* === cholmod_ssmult ======================================================= */
/* ========================================================================== */

cholmod_sparse *CHOLMOD(ssmult)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* left matrix to multiply */
    cholmod_sparse *B,  /* right matrix to multiply */
    int stype,      /* requested stype of C */
    int values,     /* TRUE: do numerical values, FALSE: pattern only */
    int sorted,     /* if TRUE then return C with sorted columns */
    /* --------------- */
    cholmod_common *Common
)
{
    double bjt ;
    double *Ax, *Bx, *Cx, *W ;
    Int *Ap, *Anz, *Ai, *Bp, *Bnz, *Bi, *Cp, *Ci, *Flag ;
    cholmod_sparse *C, *A2, *B2, *A3, *B3, *C2 ;
    Int apacked, bpacked, j, i, pa, paend, pb, pbend, ncol, mark, cnz, t, p,
    nrow, anz, bnz, do_swap_and_transpose, n1, n2 ;

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

    RETURN_IF_NULL_COMMON (NULL) ;
    RETURN_IF_NULL (A, NULL) ;
    RETURN_IF_NULL (B, NULL) ;
    values = values &&
    (A->xtype != CHOLMOD_PATTERN) && (B->xtype != CHOLMOD_PATTERN) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, 
        values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
    RETURN_IF_XTYPE_INVALID (B, CHOLMOD_PATTERN, 
        values ? CHOLMOD_REAL : CHOLMOD_ZOMPLEX, NULL) ;
    if (A->ncol != B->nrow)
    {
    /* inner dimensions must agree */
    ERROR (CHOLMOD_INVALID, "A and B inner dimensions must match") ;
    return (NULL) ;
    }
    /* A and B must have the same numerical type if values is TRUE (both must
     * be CHOLMOD_REAL, this is implicitly checked above) */
    Common->status = CHOLMOD_OK ;

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

    if (A->nrow <= 1)
    {
    /* C will be implicitly sorted, so no need to sort it here */
    sorted = FALSE ;
    }
    if (sorted)
    {
    n1 = MAX (A->nrow, B->ncol) ;
    }
    else
    {
    n1 = A->nrow ;
    }
    n2 = MAX4 (A->ncol, A->nrow, B->nrow, B->ncol) ;
    CHOLMOD(allocate_work) (n1, n2, values ? n1 : 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    /* out of memory */
    return (NULL) ;
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1 : 0, Common)) ;

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

    /* convert A to unsymmetric, if necessary */
    A2 = NULL ;
    B2 = NULL ;
    if (A->stype)
    {
    /* workspace: Iwork (max (A->nrow,A->ncol)) */
    A2 = CHOLMOD(copy) (A, 0, values, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        /* out of memory */
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;
        return (NULL) ;
    }
    A = A2 ;
    }

    /* convert B to unsymmetric, if necessary */
    if (B->stype)
    {
    /* workspace: Iwork (max (B->nrow,B->ncol)) */
    B2 = CHOLMOD(copy) (B, 0, values, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
        /* out of memory */
        CHOLMOD(free_sparse) (&A2, Common) ;
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;
        return (NULL) ;
    }
    B = B2 ;
    }

    ASSERT (CHOLMOD(dump_sparse) (A, "A", Common) >= 0) ;
    ASSERT (CHOLMOD(dump_sparse) (B, "B", Common) >= 0) ;

    /* get the A matrix */
    Ap  = A->p ;
    Anz = A->nz ;
    Ai  = A->i ;
    Ax  = A->x ;
    apacked = A->packed ;

    /* get the B matrix */
    Bp  = B->p ;
    Bnz = B->nz ;
    Bi  = B->i ;
    Bx  = B->x ;
    bpacked = B->packed ;

    /* get the size of C */
    nrow = A->nrow ;
    ncol = B->ncol ;

    /* get workspace */
    W = Common->Xwork ;     /* size nrow, unused if values is FALSE */
    Flag = Common->Flag ;   /* size nrow, Flag [0..nrow-1] < mark on input*/

    /* ---------------------------------------------------------------------- */
    /* count the number of entries in the result C */
    /* ---------------------------------------------------------------------- */

    cnz = 0 ;
    for (j = 0 ; j < ncol ; j++)
    {
    /* clear the Flag array */
    mark = CHOLMOD(clear_flag) (Common) ;

    /* for each nonzero B(t,j) in column j, do: */
    pb = Bp [j] ;
    pbend = (bpacked) ? (Bp [j+1]) : (pb + Bnz [j]) ;
    for ( ; pb < pbend ; pb++)
    {
        /* B(t,j) is nonzero */
        t = Bi [pb] ;

        /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
        pa = Ap [t] ;
        paend = (apacked) ? (Ap [t+1]) : (pa + Anz [t]) ;
        for ( ; pa < paend ; pa++)
        {
        i = Ai [pa] ;
        if (Flag [i] != mark)
        {
            Flag [i] = mark ;
            cnz++ ;
        }
        }
    }
    if (cnz < 0)
    {
        break ;     /* integer overflow case */
    }
    }

    mark = CHOLMOD(clear_flag) (Common) ;

    /* ---------------------------------------------------------------------- */
    /* check for integer overflow */
    /* ---------------------------------------------------------------------- */

    if (cnz < 0)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    CHOLMOD(free_sparse) (&A2, Common) ;
    CHOLMOD(free_sparse) (&B2, Common) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;
    return (NULL) ;
    }

    /* ---------------------------------------------------------------------- */
    /* Determine how to return C sorted (if requested) */
    /* ---------------------------------------------------------------------- */

    do_swap_and_transpose = FALSE ;

    if (sorted)
    {
    /* Determine the best way to return C with sorted columns.  Computing
     * C = (B'*A')' takes cnz + anz + bnz time (ignoring O(n) terms).
     * Sorting C when done, C = (A*B)'', takes 2*cnz time.  Pick the one
     * with the least amount of work. */

    anz = CHOLMOD(nnz) (A, Common) ;
    bnz = CHOLMOD(nnz) (B, Common) ;

    do_swap_and_transpose = (anz + bnz < cnz) ;

    if (do_swap_and_transpose)
    {

        /* -------------------------------------------------------------- */
        /* C = (B'*A')' */
        /* -------------------------------------------------------------- */

        /* workspace: Iwork (A->nrow) */
        A3 = CHOLMOD(ptranspose) (A, values, NULL, NULL, 0, Common) ;
        CHOLMOD(free_sparse) (&A2, Common) ;
        A2 = A3 ;
        if (Common->status < CHOLMOD_OK)
        {
        /* out of memory */
        CHOLMOD(free_sparse) (&A2, Common) ;
        CHOLMOD(free_sparse) (&B2, Common) ;
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common));
        return (NULL) ;
        }
        /* workspace: Iwork (B->nrow) */
        B3 = CHOLMOD(ptranspose) (B, values, NULL, NULL, 0, Common) ;
        CHOLMOD(free_sparse) (&B2, Common) ;
        B2 = B3 ;
        if (Common->status < CHOLMOD_OK)
        {
        /* out of memory */
        CHOLMOD(free_sparse) (&A2, Common) ;
        CHOLMOD(free_sparse) (&B2, Common) ;
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common));
        return (NULL) ;
        }
        A = B2 ;
        B = A2 ;

        /* get the new A matrix */
        Ap  = A->p ;
        Anz = A->nz ;
        Ai  = A->i ;
        Ax  = A->x ;
        apacked = A->packed ;

        /* get the new B matrix */
        Bp  = B->p ;
        Bnz = B->nz ;
        Bi  = B->i ;
        Bx  = B->x ;
        bpacked = B->packed ;

        /* get the size of C' */
        nrow = A->nrow ;
        ncol = B->ncol ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* allocate C */
    /* ---------------------------------------------------------------------- */

    C = CHOLMOD(allocate_sparse) (nrow, ncol, cnz, FALSE, TRUE, 0,
        values ? A->xtype : CHOLMOD_PATTERN, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    /* out of memory */
    CHOLMOD(free_sparse) (&A2, Common) ;
    CHOLMOD(free_sparse) (&B2, Common) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;
    return (NULL) ;
    }

    Cp = C->p ;
    Ci = C->i ;
    Cx = C->x ;

    /* ---------------------------------------------------------------------- */
    /* C = A*B */
    /* ---------------------------------------------------------------------- */

    cnz = 0 ;

    if (values)
    {

    /* pattern and values */
    for (j = 0 ; j < ncol ; j++)
    {
        /* clear the Flag array */
        mark = CHOLMOD(clear_flag (Common)) ;

        /* start column j of C */
        Cp [j] = cnz ;

        /* for each nonzero B(t,j) in column j, do: */
        pb = Bp [j] ;
        pbend = (bpacked) ? (Bp [j+1]) : (pb + Bnz [j]) ;
        for ( ; pb < pbend ; pb++)
        {
        /* B(t,j) is nonzero */
        t = Bi [pb] ;
        bjt = Bx [pb] ;

        /* add the nonzero pattern of A(:,t) to the pattern of C(:,j)
         * and scatter the values into W */
        pa = Ap [t] ;
        paend = (apacked) ? (Ap [t+1]) : (pa + Anz [t]) ;
        for ( ; pa < paend ; pa++)
        {
            i = Ai [pa] ;
            if (Flag [i] != mark)
            {
            Flag [i] = mark ;
            Ci [cnz++] = i ;
            }
            W [i] += Ax [pa] * bjt ;
        }
        }

        /* gather the values into C(:,j) */
        for (p = Cp [j] ; p < cnz ; p++)
        {
        i = Ci [p] ;
        Cx [p] = W [i] ;
        W [i] = 0 ;
        }
    }

    }
    else
    {

    /* pattern only */
    for (j = 0 ; j < ncol ; j++)
    {
        /* clear the Flag array */
        mark = CHOLMOD(clear_flag) (Common) ;

        /* start column j of C */
        Cp [j] = cnz ;

        /* for each nonzero B(t,j) in column j, do: */
        pb = Bp [j] ;
        pbend = (bpacked) ? (Bp [j+1]) : (pb + Bnz [j]) ;
        for ( ; pb < pbend ; pb++)
        {
        /* B(t,j) is nonzero */
        t = Bi [pb] ;

        /* add the nonzero pattern of A(:,t) to the pattern of C(:,j) */
        pa = Ap [t] ;
        paend = (apacked) ? (Ap [t+1]) : (pa + Anz [t]) ;
        for ( ; pa < paend ; pa++)
        {
            i = Ai [pa] ;
            if (Flag [i] != mark)
            {
            Flag [i] = mark ;
            Ci [cnz++] = i ;
            }
        }
        }
    }
    }

    Cp [ncol] = cnz ;
    ASSERT (MAX (1,cnz) == C->nzmax) ;

    /* ---------------------------------------------------------------------- */
    /* clear workspace and free temporary matrices */
    /* ---------------------------------------------------------------------- */

    CHOLMOD(free_sparse) (&A2, Common) ;
    CHOLMOD(free_sparse) (&B2, Common) ;
    CHOLMOD(clear_flag) (Common) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;

    /* ---------------------------------------------------------------------- */
    /* convert C to a symmetric upper/lower matrix if requested */
    /* ---------------------------------------------------------------------- */

    /* convert C in place, which cannot fail since no memory is allocated */
    if (stype > 0)
    {
    /* C = triu (C), in place */
    (void) CHOLMOD(band_inplace) (0, ncol, values, C, Common) ;
    C->stype = 1 ;
    }
    else if (stype < 0)
    {
    /* C = tril (C), in place */
    (void) CHOLMOD(band_inplace) (-nrow, 0, values, C, Common) ;
    C->stype = -1 ;
    }
    ASSERT (Common->status >= CHOLMOD_OK) ;

    /* ---------------------------------------------------------------------- */
    /* sort C, if requested */
    /* ---------------------------------------------------------------------- */

    if (sorted)
    {
    if (do_swap_and_transpose)
    {
        /* workspace: Iwork (C->ncol), which is A->nrow since C=(B'*A') */
        C2 = CHOLMOD(ptranspose) (C, values, NULL, NULL, 0, Common) ;
        CHOLMOD(free_sparse) (&C, Common) ;
        if (Common->status < CHOLMOD_OK)
        {
        /* out of memory */
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common));
        return (NULL) ;
        }
        C = C2 ;
    }
    else
    {
        /* workspace: Iwork (max (C->nrow,C->ncol)) */
        if (!CHOLMOD(sort) (C, Common))
        {
        /* out of memory */
        CHOLMOD(free_sparse) (&C, Common) ;
        ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common));
        return (NULL) ;
        }
    }
    }

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

    DEBUG (CHOLMOD(dump_sparse) (C, "ssmult", Common) >= 0) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, values ? n1:0, Common)) ;
    return (C) ;
}
#endif