Sparse matrices are typically represented by storing only their
non-zero entries. Along with the actual matrix values, their
row and column information needs to be explicitly stored.
For example, the sparse matrix
- -
| 1 1 0 0 0 0 0 |
| 0 0 0 0 0 0 0 |
| 0 0 3 3 3 0 0 |
| 0 6 0 0 0 4 0 |
| 2 0 0 0 0 5 5 |
- -
may be stored column-by-column by the three arrays:
Column pointers (colptrs) : 0 2 4 5 6 7 9 10
Row indices (rowindices) : 0 4 0 3 2 2 2 3 4 4
Non-zero values (val) : 1 2 1 6 3 3 3 4 5 5
(In the above, note that all indices start at 0 (C-style) )
Assume the sparse matrix has m rows and n columns.
"colptrs" is an integer array of size n+1. This array
indexes into the "rowindices" and "val" arrays. In particular,
the i-th column starts at colptrs(i) and ends at colptrs(i+1)-1.
For programming convenience, colptrs(n+1) is set to nz, where
nz is the number of nonzeros in the matrix.
"rowindices" is an integer array of size nz (the number of
nonzeros in the matrix), and it contains the row indices
of the matrix elements. For example, rowindices[colptrs[i]]
gives the row index of the 1st non-zero element in the i-th column.
"val" is an array of size nz (the number of nonzeros in the
matrix), and it contains the matrix values. Whether val is an
integer, float or double (or whatever) array depends on the type
of your matrix elements. For example, val[colptrs[i]] gives
the value of the 1st non-zero element in the i-th column.
Since many of you will be exchanging sparse matrices, I would
like you to input and output the m x n sparse matrices (to files) in
the following form:
5 7 10 % m, n, nz
% Note that m is number of rows
% and n is the number of columns
0 2 4 5 6 7 9 10 % colptrs
0 4 0 3 2 2 2 3 4 4 % rowindices
1.0 2.0 1.0 6.0 3.0 3.0 3.0 4.0 5.0 5.0 % val
The following is an example of using this sparse matrix format:
void at_times_x(int n, int nz, int *colptrs, int *rowindices, double *val,
double *x, double *y)
/* Subroutine at_times_x computes y = A^T * x, A is mxn and stored as
a sparse matrix, and y and x are dense vectors */
{ double yi;
int i,j;
for(i = 0;i < n;i++){
yi = 0.0;
for(j = colptrs[i];j < colptrs[i+1];j++){
yi += val[j] * x[rowindices[j]]; }
y[i] = yi; }
}