A.2 Sparse Matrices

Sparse matrices are stored in compressed column storage (CCS) format. For a general nrows by ncols sparse matrix with nnz nonzero entries this means the following. The sparsity pattern and the nonzero values are stored in three fields:

values:
A ’d’ or ’z’ matrix of size (nnz,1) with the nonzero entries of the matrix stored columnwise.
rowind:
An array of integers of length nnz containing the row indices of the nonzero entries, stored in the same order as values.
colptr:
An array of integers of length ncols+1 with for each column of the matrix the index of the first element in values from that column. More precisely, colptr[0] is 0, and for k = 0, 1, …, ncols-1, colptr[k+1] is equal to colptr[k] plus the number of nonzeros in column k of the matrix. Thus, colptr[ncols] is equal to nnz, the number of nonzero entries.

For example, for the matrix

    ⌊ 1  0  0 5 ⌋
    | 2  0  4 0 |
A = |⌈ 0  0  0 6 |⌉
      3  0  0 0

the elements of values, rowind and colptr are:

values: 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, rowind: 0, 1,3, 1, 0, 2, colptr: 0, 3, 3, 4, 6.

It is crucial that for each column the row indices in rowind are sorted; the equivalent representation

values: 3.0, 2.0, 1.0, 4.0, 5.0, 6.0, rowind: 3, 1, 0, 1, 0, 2, colptr: 0, 3, 3, 4, 6

is not allowed (and will likely cause the program to crash).

The nzmax field specifies the number of non-zero elements the matrix can store. It is equal to the length of rowind and values; this number can be larger that colptr[nrows], but never less. This field makes it possible to preallocate a certain amount of memory to avoid reallocations if the matrix is constructed sequentially by filling in elements. In general the nzmax field can safely be ignored, however, since it will always be adjusted automatically as the number of non-zero elements grows.

The id field controls the type of the matrix and can have values DOUBLE and COMPLEX.

Sparse matrices are created using the following functions from the API.

spmatrix *SpMatrix_New(int nrows, int ncols, int nzmax, int id)

Returns a sparse zero matrix with nrows rows and ncols columns. nzmax is the number of elements that will be allocated (the length of the values and rowind fields).

spmatrix *SpMatrix_NewFromMatrix(spmatrix *src, int id)

Returns a copy the sparse matrix src.

spmatrix *SpMatrix_NewFromIJV(matrix *I, matrix *J, matrix *V, int nrows, int ncols, int nzmax, int id)

Creates a sparse matrix with nrows rows and ncols columns from a triplet description. I and J must be integer matrices and V either a double or complex matrix, or NULL. If V is NULL the values of the entries in the matrix are undefined, otherwise they are specified by V. Repeated entries in V are summed. The number of allocated elements is given by nzmax, which is adjusted if it is smaller than the required amount.

We illustrate use of the sparse matrix class by listing the source code for the real() method, which returns the real part of a sparse matrix:

static PyObject * spmatrix_real(spmatrix *self) {  
 
  if (SP_ID(self) != COMPLEX)  
    return (PyObject *)SpMatrix_NewFromMatrix(self, 0, SP_ID(self));  
 
  spmatrix *ret = SpMatrix_New(SP_NROWS(self), SP_NCOLS(self),  
      SP_NNZ(self), DOUBLE);  
  if (!ret) return PyErr_NoMemory();  
 
  int i;  
  for (i=0; i < SP_NNZ(self); i++)  
    SP_VALD(ret)[i] = creal(SP_VALZ(self)[i]);  
 
  memcpy(SP_COL(ret), SP_COL(self), (SP_NCOLS(self)+1)*sizeof(int_t));  
  memcpy(SP_ROW(ret), SP_ROW(self), SP_NNZ(self)*sizeof(int_t));  
  return (PyObject *)ret;  
}