The nonzero entries A(i,j) = 10*i + j are chosen to make it easier to see where the entry is located; the first digit gives the row number and the second one gives the column number. For the problems and iterative methods we will use, the values are actually double precision floating point numbers. Links back to this example are in each data structure section, to make it easier to return to it. I recommend that you open two windows (or at least two tabs) to keep the example in view while reading about each data structure.

The first three of the data structures below require two 1D arrays of integers and one 1D array of doubles. In the tables displaying those arrays, to make it easier to keep track of the correct indices an extra table row proceeds them, in a greenish color. It is not part of the data structures, just a convenience similar to having column headers in a spreadsheet. The tables displaying the arrays of integers and doubles are split in COO and CSR, in case you have a narrow browser window. The names of the arrays used below are mostly standard, although the three (rowptr, col, val) in CSR are more commonly named (ia, ja, a). In the text empty parentheses are appended to "col" and "val" to emphasize that they are arrays. When those terms are embedded in text with a specific index (e.g., col(k)) , or are written as a triplet (row, col, val) the parentheses are omitted. So "col( )" and "val( )" just means they are arrays, not functions.

The example matrix above has nnz = 21 nonzeros, but some arrays are shown with an unused nnz+1 = 22^nd entry. In practice the arrays are allocated with only nnz, not nnz+1, components. For debugging purposes only, you might allocate nnz+1 entries instead and then set row(nnz+1) = col(nnz+1) = -5 and val(nnz+1) = NaN; that will cause your code to stop with an error condition if you accidentally run off the end of the array. If you do this, be sure to change the codes back to using only nnz entries before handing anything in.

Coordinate-wise (COO)

COO features:

• Easy to use, understand, debug.
• Easy to add nonzeros to data structure (for fill in)
• Rows/columns need not be ordered in any way.
• Storage: nnz floating point numbers and 2nnz integers.
• Row/column access equally easy.
• The matrix-vector product y = A*x becomes:

Initialize y(1:n) = 0

for k = 1 to nnz

   y(row(k)) = y(row(k)) + val(k)*x(col(k))

end for

For the example matrix listed above, the matrix would be stored as
 

Don't forget: in the table, the row labeled "Index" is for convenience, and is not part of the data structure. Also, for the example matrix the three arrays (row, col, val) in COO storage would be of length 21, not 22. Although the arrays show the nonzeros in row order, in practice the nonzeros could be placed in any order - row after row, column after column, or even randomly.

Compressed Sparse Row (CSR)

CSR stores the column indices and values just like COO, but entries must be listed one row after another. Doing that reduces the required information about row indices to just knowing where a row begins in the column index and value arrays col( ) and val( ). This is done using an integer array of length n+1, such that the k-th entry gives the index in col( ) and val( ) where the k-th row begins. That accounts for the first n entries in the row index information array (which is usually called "row pointers" or in the tables below, rowptr( ) ). The last entry, rowptr(n+1), has the value nnz + 1. That allows the loops implementing a matrix*vector product without having a test for the special case where the last row is being used, as shown below.

Beware that the term "row pointer" is common, but what the integer array row( ) really contains are the indices (usually 1-based) of where each row starts in the arrays col( ) and val( ). It does not contain a pointer in the usual C/C++/Fortran sense.

Features:

• Requires entries from a row to be consecutively located in data structure.
• Columns need not be ordered in any way.
• Storage: nnz floating point, nnz + n + 1 integer.
• Row access is easy, column access is difficult.
• Matrix-vector products y = A*x pseudo-code:

Initialize y(1:n) = 0
for i = 1 to n
   for k = row(i) to row(i+1)-1
       y(i) = y(i) + val(k)*x(col(k))
   end for
end for

For the example matrix at the top of the page, the CSR arrays are:

The arrays val( ) and col( ) above have been split into two separate tables with the rowptr( ) array between them to make it easier to see how the rowptr( ) array of indices matches up. As a simple example, the fourth row of the matrix begins with the index row(4) = 10 and ends with the index row(4+1) - 1 = row(5) - 1 = 13 - 1 = 12. What would the entries in row( ) be if the fourth row was all zeros and none of them are explicitly stored? Will the matrix-vector product pseudo-code still work correctly?

Compressed Sparse column (CSC)

Modified Sparse Row (MSR)

The first n entries in the integer array bindx stores the indices where each row's off-diagonal entries begin. The n+1 entry contains nnz + 1 so that (as with CSR) the matrix-vector product can be written with a single nested loop and no special cases handled for the last row of the matrix. The array of values stores first the n diagonal entries, then skips a single unused entry in position n+1, followed by the off-diagonal entries, row-by-row.

For the example matrix at the top of the page, MSR storage is:
 

You should try to write what the sparse matrix-vector multiplication looks like in MSR.

Something to beware of with COO, CSR, and MSR: all three of these were originally used with Fortran, and hence the indices into the arrays are sometimes 1-based. However, many modern packages are written in C/C++, so sometimes 0-based indexing is used. To make it even more confusing, to maintain compatibility with Fortran drivers, some C/C++ packages use 1-based indexing. You need to carefully check the documentation, and sometimes even the source code, to determine which is used.

Sometimes a C++ package will assume the indexing is 1-based, and shift the integer arrays accordingly to use the more convenient 0-based C indexing internally. Then the package shifts indices back to 1-based before returning to the calling program unit. Others have the user pass in a flag that indicates whether or not the data structure is using 0- or 1-based indexing. Most of the recent ones allow the data to be input with either 0 or 1 based indexing, and the user provides a flag indicating if the index is "C style" or "Fortran style". RTFM (which I believe is the acronym for "read the fine manual").

Compressed Diagonal Storage (CDS)

• Need only store array of values and two integers: p = lower and q = upper bandwidth of matrix. The array containing the nonzero values is indexed as val(-p:q,1:n).
• Easy to add entries to data structure inside of bandwidth.
• Storage: p+q+1 floating point numbers and 2 integers (for p and q).
• Row and column access are equally difficult.
• This storage scheme is also called "banded storage", and is used by LAPACK and LINPACK in routines that have GB as the second and third letters in their name. For example, DGBSV solves a general linear system Ax = b, where A is stored in CDS format.
• Matrix-vector products y = A*x are performed as:

Initialize y(1:n) = 0
for diag = -p to q
   for k = max(1,1-diag) to min(n,n-diag)
       y(k) = y(k) + val(k,diag)*x(k+diag)
   end for
end for

The first row and column (labeled "Diag") in the table are not part of the storage scheme, and are just for convenience in reading the table.  Note that the numbering of the diagonals,  with the main diagonal as zero, subdiagonals with negative numbers, and superdiagonals positive numbers, matches what MATLAB does in its tril( ) and triu( ) operators.

Entries marked with an asterisk (*) are nonexistent ones. For example, the (1,1) entry in the array above would correspond to the (1,-1) entry in A, which does not exist. Similarly the (6,8) entry in the array would be the A(8,11) but A is only 8x8. Any algorithm should not reference those entries, or require them to be set to zero. That is the reason for the max( ) and min( ) operators in the loop controls in the matrix-vector product. It also provides an easy check for some indexing errors: set those entries to NaN, and set the compiler flags so that encountering NaN will trigger a stopping error.

Envelope or skyline storage

Three arrays are required: val(1:nenv) for nonzero values, rowptr(1:n+1) for row pointers, and fstcol(1:n) for column number of first nonzero in each row (fstcol is shorthand for "first column"). Matrix-vector products become

for k = 1:n
  y(k) = 0
  xcol = fstcol(k)
  istart = rowptr(i)
  iend = rowptr(i+1)-1
  for j = istart:iend
    y(k) = y(k) + val(j)*x(xcol+j-istart)
  end for
end for

For the example matrix, skyline storage is :
 

• Started: Sun 05 Feb 2012, 04:34 PM
• Last Modified: Sun 05 Feb 2012, 08:13 PM