build_csr.c: Build a Compressed Sparse Row matrix

To build a CSR matrix, first construct a Coordinate-form (COO) matrix, then use function coo_to_csr_matrix() to convert it into a CSR matrix. To build the COO matrix, first use create_coo_matrix() and then add individual elements using add_to_coo_matrix(). create_coo_matrix(n,r,c) Creates an empty COO matrix with enough space for up to n nonzero elements. All elements will be triples (i,j,x) where 0<=i<r and 0<=j<c.

struct coo_matrix *create_coo_matrix (unsigned maxn, unsigned rows, unsigned cols) {
  struct coo_matrix *p = surely_malloc (sizeof (*p));
  p->row_ind = (unsigned *)surely_malloc(maxn * sizeof(unsigned));
  p->col_ind = (unsigned *)surely_malloc(maxn * sizeof(unsigned));
  p->val = (double *)surely_malloc(maxn * sizeof(double));
  p->n=0;
  p->maxn=maxn;
  p->rows=rows; p->cols=cols;
  return p;
}

add_to_coo_matrix(p,i,j,x) Add the triple (i,j,x) to the COO matrix. Multiple elements at the same position are permitted, with the semantics that they are to be added together. That is, if another triple (i,j,y) is already present with the same i,j, the meaning is that the element at (i,j) has value (y+x).

void add_to_coo_matrix(struct coo_matrix *p, unsigned i, unsigned j, double x) {
  unsigned n = p->n;
  if (n>=p->maxn) exit(2);
  p->row_ind[n]=i;
  p->col_ind[n]=j;
  p->val[n]=x;
  p->n = n+1;
}

The functions swap, coo_less, coo_quicksort, coo_count are internal functions not part of the API.

/* adapted from qsort3 in cbench:
   https://github.com/cverified/cbench/blob/master/src/qsort/qsort3.c */

/* sort the coordinate elements of a coo_matrix */
void coo_quicksort(struct coo_matrix *p, unsigned base, unsigned n)
{
}

/* Count the number of distinct row/col entries in a sorted coordinate list */
unsigned coo_count (struct coo_matrix *p) {
  unsigned i,r,c,ri,ci,n,count;
  unsigned *row_ind = p->row_ind, *col_ind = p->col_ind;
  n = p->n;
  count=0;
  r=-1;
  c=0;
  for (i=0; i<n; i++) {
    ri=row_ind[i]; 
    ci=col_ind[i];
    if (ri!=r || ci!=c) {
      count++;
      r=ri; c=ci;
    }
  }
  return count;
}

coo_to_csr_matrix(p) Given a COO matrix p, convert to a CSR matrix. This takes NlogN time, where N is the number of triples in the COO matrix, because the first step is to sort the triples by lexicographic order of (i,j). This function has the side effect of rearranging (sorting) the triples.

struct csr_matrix *coo_to_csr_matrix(struct coo_matrix *p) {
  struct csr_matrix *q;
  unsigned count, i;
  unsigned r,c, ri, ci, cols, k, l, rows;
  unsigned *col_ind, *row_ptr, *prow_ind, *pcol_ind;
  double x, *val, *pval;
  unsigned n = p->n;
  coo_quicksort(p, 0, n);
  k = coo_count(p);
  rows = p->rows;
  prow_ind=p->row_ind;
  pcol_ind=p->col_ind;
  pval = p->val;
  q = surely_malloc(sizeof(struct csr_matrix));
  val = surely_malloc(k * sizeof(double));
  col_ind = surely_malloc(k * sizeof(unsigned));
  row_ptr = surely_malloc ((rows+1) * sizeof(unsigned));
  r=-1;
  c=0; /* this line is unnecessary for correctness but simplifies the proof */
  l=0;
  /* partial_csr_0 */
  for (i=0; i<n; i++) {
    ri = prow_ind[i];
    ci = pcol_ind[i];
    x = pval[i];
    if (ri==r)
      if (ci==c)
    val[l-1] += x; /* partial_CSR_duplicate */
      else {
    c=ci;
    col_ind[l] = ci;
    val[l] = x;
    l++;           /* partial_CSR_newcol */
      }
    else {
      while (r+1<=ri) row_ptr[++r]=l; /* partial_CSR_skiprow */
      c= ci;
      col_ind[l] = ci;
      val[l] = x;
      l++;            /* partial_CSR_newrow */
    }
  }
  cols = p->cols;
  while (r+1<=rows) row_ptr[++r]=l;  /* partial_CSR_lastrows */
  q->val = val;
  q->col_ind = col_ind;
  q->row_ptr = row_ptr;
  q->rows = rows;
  q->cols = cols;
  return q;          /* partial_CSR_properties */
}