The R Project SVN R-packages

/* Based on UMF_triplet.c from UMFPACK which carries the following copyright  */

/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A.      */
/* Davis.  All Rights Reserved.  See ../README for License.                   */
/* email: davis@cise.ufl.edu    CISE Department, Univ. of Florida.            */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
/* -------------------------------------------------------------------------- */

#include "triplet_to_col.h"
void triplet_to_col
(
    int n_row,
    int n_col,
    int nz,
    const int Ti [ ],       /* size nz */
    const int Tj [ ],       /* size nz */
    const double Tx [ ],    /* size nz */
    int Ap [ ],         /* size n_col + 1 */
    int Ai [ ],         /* size nz */
    double Ax [ ]       /* size nz */
)
{
    int i, j, k, p, cp, p1, p2, pdest, pj;
    int *Rp = Calloc((n_row + 1), int), 
    *Rj = Calloc(nz, int),
    *W = Calloc((n_row > n_col) ? n_row : n_col, int),
    *RowCount = Calloc(n_row, int);
    double *Rx = (double *) NULL;

    if (Tx) Rx = Calloc(nz, double);

    /* count the entries in each row (including duplicates) */
    /* use W as workspace for row counts (including duplicates) */
    memset(W, 0, sizeof(int) * n_row);
    for (k = 0 ; k < nz ; k++) {
    i = Ti [k] ;
    j = Tj [k] ;
    if (i < 0 || i >= n_row || j < 0 || j >= n_col)
        error(_("entry %d in matrix[%d,%d] has row %d and column %d"),
          k, n_row, n_col, i, j);
    W [i]++ ;
    }
    /* compute the row pointers */
    Rp [0] = 0 ;
    for (i = 0 ; i < n_row ; i++) {
    Rp [i+1] = Rp [i] + W [i] ;
    W [i] = Rp [i] ;
    }
    /* W is now equal to the row pointers */

    /* ---------------------------------------------------------------------- */
    /* construct the row form */
    /* ---------------------------------------------------------------------- */
    for (k = 0 ; k < nz ; k++) {
    p = W [Ti [k]]++ ;
    Rj [p] = Tj [k] ;
    if (Tx) Rx [p] = Tx [k] ;
    }
    /* Rp stays the same, but W [i] is advanced to the start of row i+1 */
    /* ---------------------------------------------------------------------- */
    /* sum up duplicates */
    /* ---------------------------------------------------------------------- */
    /* use W [j] to hold position in Ri/Rx/Rz of a_ij, for row i [ */
    for (j = 0 ; j < n_col ; j++) {
    W [j] = -1;
    }
    for (i = 0 ; i < n_row ; i++) {
    p1 = Rp [i] ;
    p2 = Rp [i+1] ;
    pdest = p1 ;
    /* At this point, W [j] < p1 holds true for all columns j, */
    /* because Ri/Rx/Rz is stored in row oriented order. */
    for (p = p1; p < p2; p++) {
        j = Rj [p] ;
        pj = W [j] ;
        if (pj >= p1) {
        /* this column index, j, is already in row i, at position pj */
        /* sum the entry */
        if (Tx) Rx [pj] += Rx [p] ;
        } else {
        /* keep the entry */
        /* also keep track in W[j] of position of a_ij for case above */
        W [j] = pdest ;
        /* no need to move the entry if pdest is equal to p */
        if (pdest != p)
        {
            Rj [pdest] = j ;
            if (Tx) Rx [pdest] = Rx [p] ;
        }
        pdest++ ;
        }
    }
    RowCount [i] = pdest - p1 ;
    }
    /* done using W for position of a_ij ] */
    /* ---------------------------------------------------------------------- */
    /* count the entries in each column */
    /* ---------------------------------------------------------------------- */
    /* [ use W as work space for column counts of A */
    memset(W, 0, sizeof(int) * n_col);
    for (i = 0 ; i < n_row ; i++) {
    for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++) {
        j = Rj [p] ;
        W [j]++ ;
    }
    }
    /* ---------------------------------------------------------------------- */
    /* create the column pointers */
    /* ---------------------------------------------------------------------- */
    Ap [0] = 0 ;
    for (j = 0 ; j < n_col ; j++) {
    Ap [j+1] = Ap [j] + W [j] ;
    }
    /* done using W as workspace for column counts of A ] */
    for (j = 0 ; j < n_col ; j++) {
    W [j] = Ap [j] ;
    }
    /* ---------------------------------------------------------------------- */
    /* construct the column form */
    /* ---------------------------------------------------------------------- */

    for (i = 0 ; i < n_row ; i++) {
    for (p = Rp [i] ; p < Rp [i] + RowCount [i] ; p++) {
        cp = W [Rj [p]]++ ;
        Ai [cp] = i ;
        if (Tx) Ax [cp] = Rx [p] ;
    }
    }
    Free(Rp); Free(Rj); Free(W); Free(RowCount);
    if (Tx) Free(Rx);
}