Revision f584db857045d794fb5a587b18f1ca3ed40b8679 authored by Douglas Bates on 15 January 2006, 00:00:00 UTC, committed by Gabor Csardi on 15 January 2006, 00:00:00 UTC
1 parent a2d2fd0
triplet_to_col.c
/* 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);
}
Computing file changes ...