Sparse Matrix Storage Formats

Storage Formats for the Direct Sparse Solvers

The storing the non-zero elements of a sparse matrix into a linear array is done by walking down each column (column-major format) or across each row (row-major format) in order, and writing the non-zero elements to a linear array in the order they appear in the walk.

For symmetric matrices, it is necessary to store only the upper triangular half of the matrix (upper triangular format) or the lower triangular half of the matrix (lower triangular format).

The Intel MKL direct sparse solvers use a row-major upper triangular storage format: the matrix is compressed row-by-row and for symmetric matrices only non-zero elements in the upper triangular half of the matrix are stored.

The Intel MKL sparse matrix storage format for direct sparse solvers is specified by three arrays: values, columns, and rowIndex. The following table describes the arrays in terms of the values, row, and column positions of the non-zero elements in a sparse matrix.

values

A real or complex array that contains the non-zero elements of a sparse matrix. The non-zero elements are mapped into the values array using the row-major upper triangular storage mapping described above.

columns

Element i of the integer array columns is the number of the column that contains the i-th element in the values array.

rowIndex

Element j of the integer array rowIndex gives the index of the element in the values array that is first non-zero element in a row j.

The length of the values and columns arrays is equal to the number of non-zero elements in the matrix.

As the rowIndex array gives the location of the first non-zero element within a row, and the non-zero elements are stored consecutively, the number of non-zero elements in the i-th row is equal to the difference of rowIndex(i) and rowIndex(i+1).

To have this relationship hold for the last row of the matrix, an additional entry (dummy entry) is added to the end of rowIndex. Its value is equal to the number of non-zero elements plus one. This makes the total length of the rowIndex array one larger than the number of rows in the matrix.

Note

The Intel MKL sparse storage scheme for the direct sparse solvers supports both with one-based indexing and zero-based indexing.

Consider the symmetric matrix A:

Only elements from the upper triangle are stored. The actual arrays for the matrix A are as follows:

For a non-symmetric or non-Hermitian matrix, all non-zero elements need to be stored. Consider the non-symmetric matrix B:

The matrix B has 13 non-zero elements, and all of them are stored as follows:

Direct sparse solvers can also solve symmetrically structured systems of equations. A symmetrically structured system of equations is one where the pattern of non-zero elements is symmetric. That is, a matrix has a symmetric structure if a(j,i) is not zero if and only if a(j, i) is not zero. From the point of view of the solver software, a "non-zero" element of a matrix is any element stored in the values array, even if its value is equal to 0. In that sense, any non-symmetric matrix can be turned into a symmetrically structured matrix by carefully adding zeros to the values array. For example, the above matrix B can be turned into a symmetrically structured matrix by adding two non-zero entries:

The matrix B can be considered to be symmetrically structured with 15 non-zero elements and represented as:

Storage Format Restrictions

The storage format for the sparse solver must conform to two important restrictions:

- the non-zero values in a given row must be placed into the values array in the order in which they occur in the row (from left to right);

- no diagonal element can be omitted from the values array for any symmetric or structurally symmetric matrix.

The second restriction implies that if symmetric or structurally symmetric matrices have zero diagonal elements, then they must be explicitly represented in the values array.

Sparse Matrix Storage Formats for Sparse BLAS Levels 2 and Level 3

This section describes in detail the sparse matrix storage formats supported in the current version of the Intel MKL Sparse BLAS Level 2 and Level 3.

CSR Format

The Intel MKL compressed sparse row (CSR) format is specified by four arrays: the values, columns, pointerB, and pointerE. The following table describes the arrays in terms of the values, row, and column positions of the non-zero elements in a sparse matrix A.

values

A real or complex array that contains the non-zero elements of A. Values of the non-zero elements of A are mapped into the values array using the row-major storage mapping described above.

columns

Element i of the integer array columns is the number of the column in A that contains the i-th value in the values array.

pointerB

Element j of this integer array gives the index of the element in the values array that is first non-zero element in a row j of A. Note that this index is equal to pointerB(j) - pointerB(1)+1 .

pointerE

An integer array that contains row indices, such that pointerE(j)-pointerB(1) is the index of the element in the values array that is last non-zero element in a row j of A.

The length of the values and columns arrays is equal to the number of non-zero elements in A.The length of the pointerB and pointerE arrays is equal to the number of rows in A.

Note

Note that the Intel MKL Sparse BLAS routines support the CSR format both with one-based indexing and zero-based indexing.

The matrix B

can be represented in the CSR format as:

This storage format is used in the NIST Sparse BLAS library [Rem05].

Note that the storage format accepted for the direct sparse solvers and described above (see Storage Formats for the Direct Sparse Solvers) is a variation of the CSR format. It also is used in the Intel MKL Sparse BLAS Level 2 both with one-based indexing and zero-based indexing. The above matrix B can be represented in this format (referred to as the 3-array variation of the CSR format) as:

The 3-array variation of the CSR format has a restriction: all non-zero elements are stored continuously, that is the set of non-zero elements in the row J goes just after the set of non-zero elements in the row J-1.

There are no such restrictions in the general (NIST) CSR format. This may be useful, for example, if there is a need to operate with different submatrices of the matrix at the same time. In this case, it is enough to define the arrays pointerB and pointerE for each needed submatrix so that all these arrays are pointers to the same array values.

Comparing the array rowIndex from the Table "Storage Arrays for a Non-Symmetric Example Matrix" with the arrays pointerB and pointerE from the Table "Storage Arrays for an Example Matrix in CSR Format" it is easy to see that

pointerB(i) = rowIndex(i) for i=1, ..5;

pointerE(i) = rowIndex(i+1) for i=1, ..5.

This enables calling a routine that has values, columns, pointerB and pointerE as input parameters for a sparse matrix stored in the format accepted for the direct sparse solvers. For example, a routine with the interface:

   Subroutine name_routine(.... ,  values, columns, pointerB, pointerE, ...)

can be called with parameters values, columns, rowIndex as follows:

   call name_routine(.... ,  values, columns, rowIndex, rowindex(2), ...).

CSC Format

The compressed sparse column format (CSC) is similar to the CSR format, but the columns are used instead the rows. In other words, the CSC format is identical to the CSR format for the transposed matrix. The CSR format is specified by four arrays: values, columns, pointerB, and pointerE. The following table describes the arrays in terms of the values, row, and column positions of the non-zero elements in a sparse matrix A.

values

A real or complex array that contains the non-zero elements of A. Values of the non-zero elements of A are mapped into the values array using the column-major storage mapping.

rows

Element i of the integer array rows is the number of the row in A that contains the i-th value in the values array.

pointerB

Element j of this integer array gives the index of the element in the values array that is first non-zero element in a column j of A. Note that this index is equal to pointerB(j) - pointerB(1)+1 .

pointerE

An integer array that contains column indices, such that pointerE(j)-pointerB(1) is the index of the element in the values array that is last non-zero element in a column j of A.

The length of the values and columns arrays is equal to the number of non-zero elements in A.The length of the pointerB and pointerE arrays is equal to the number of columns in A.

Note

Note that the Intel MKL Sparse BLAS routines support the CSC format both with one-based indexing and zero-based indexing.

The above matrix B can be represented in the CSC format as:

Coordinate Format

The coordinate format is the most flexible and simplest format for the sparse matrix representation. Only non-zero elements are stored, and the coordinates of each non-zero element are given explicitly. Many commercial libraries support the matrix-vector multiplication for the sparse matrices in the coordinate format.

The Intel MKL coordinate format is specified by three arrays: values, rows, and column, and a parameter nnz which is number of non-zero elements in A. All three arrays have dimension nnz. The following table describes the arrays in terms of the values, row, and column positions of the non-zero elements in a sparse matrix A.

values

A real or complex array that contains the non-zero elements of A in any order.

rows

Element i of the integer array rows is the number of the row in A that contains the i-th value in the values array.

columns

Element i of the integer array columns is the number of the column in A that contains the i-th value in the values array.

Note

Note that the Intel MKL Sparse BLAS routines support the coordinate format both with one-based indexing and zero-based indexing.

For example, the sparse matrix C

can be represented in the coordinate format as follows:

Diagonal Storage Format

If the sparse matrix has diagonals containing only zero elements, then the diagonal storage format can be used to reduce the amount of information needed to locate the non-zero elements. This storage format is particularly useful in many applications where the matrix arises from a finite element or finite difference discretization. The Intel MKL diagonal storage format is specified by two arrays: values and distance, and two parameters: ndiag, which is the number of non-empty diagonals, and lval, which is the declared leading dimension in the calling (sub)programs. The following table describes the arrays values and distance:

values

A real or complex two-dimensional array is dimensioned as lval by ndiag. Each column of it contains the non-zero elements of certain diagonal of A. The key point of the storage is that each element in values retains the row number of the original matrix. To achieve this diagonals in the lower triangular part of the matrix are padded from the top, and those in the upper triangular part are padded from the bottom. Note that the value of distance(i) is the number of elements to be padded for diagonal i.

distance

An integer array with dimension ndiag. Element i of the array distance is the distance between i-diagonal and the main diagonal. The distance is positive if the diagonal is above the main diagonal, and negative if the diagonal is below the main diagonal. The main diagonal has a distance equal to zero.

The above matrix C can be represented in the diagonal storage format as follows:

where the asterisks denote padded elements.

When storing symmetric, Hermitian, or skew-symmetric matrices, it is necessary to store only the upper or the lower triangular part of the matrix.

For the Intel MKL triangular solver routines elements of the array distance must be sorted in increasing order. In all other cases the diagonals and distances can be stored in arbitrary order.

Skyline Storage Format

The skyline storage format is important for the direct sparse solvers, and it is well suited for Cholesky or LU decomposition when no pivoting is required.

The skyline storage format accepted in Intel MKL can store only triangular matrix or triangular part of a matrix. This format is specified by two arrays: values and pointers. The following table describes these arrays:

values

A scalar array. For a lower triangular matrix it contains the set of elements from each row of the matrix starting from the first non-zero element to and including the diagonal element. For an upper triangular matrix it contains the set of elements from each column of the matrix starting with the first non-zero element down to and including the diagonal element. Encountered zero elements are included in the sets.

pointers

An integer array with dimension (m+1), where m is the number of rows for lower triangle (columns for the upper triangle). pointers(i) - pointers(1)+1 gives the index of element in values that is first non-zero element in row (column) i. The value of pointers(m+1) is set to nz+pointers(1), where nnz is the number of elements in the array values.

For example, the low triangle of the matrix C given above can be stored as follows:

values  =  ( 1  -2   5   4  -4   0   2   7   8   0   0   -5 )

pointers = ( 1   2   4   5   9   13 )

and the upper triangle of this matrix C can be stored as follows:

values   = ( 1  -1   5  -3   0   4   6   7  4   0   -5 )

pointers = ( 1   2   4   7   9   12 )

This storage format is supported by the NIST Sparse BLAS library [Rem05].

Note that the Intel MKL Sparse BLAS routines operating with the skyline storage format does not support general matrices.

BSR Format

The Intel MKL block compressed sparse row (BSR) format for sparse matrices is specified by four arrays: values, columns, pointerB, and pointerE. The following table describes these arrays.

values

A real array that contains the elements of the non-zero blocks of a sparse matrix. The elements are stored block-by-block in row-major order. A non-zero block is the block that contains at least one non-zero element. All elements of non-zero blocks are stored, even if some of them is equal to zero. Within each non-zero block elements are stored in column-major order in the case of one-based indexing, and in row-major order in the case of the zero-based indexing.

columns

Element i of the integer array columns is the number of the column in the block matrix that contains the i-th non-zero block.

pointerB

Element j of this integer array gives the index of the element in the columns array that is first non-zero block in a row j of the block matrix.

pointerE

Element j of this integer array gives the index of the element in the columns array that contains the last non-zero block in a row j of the block matrix plus 1.

The length of the values array is equal to the number of all elements in the non-zero blocks, the length of the columns array is equal to the number of non-zero blocks. The length of the pointerB and pointerE arrays is equal to the number of block rows in the block matrix.

Note

Note that the Intel MKL Sparse BLAS routines support BSR format both with one-based indexing and zero-based indexing.

For example, consider the sparse matrix D

If the size of the block equals 2, then the sparse matrix D can be represented as a 3x3 block matrix E with the following structure:

where

The matrix D can be represented in the BSR format as follows:

one-based indexing

values  =  (1 2 0 1 6 8 7 2 1 5 4 1 4 0 3 0 7 0 2 0)

columns  = (1   2   2   2   3)

pointerB = (1   3   4)

pointerE = (3   4   6)

zero-based indexing

values  =  (1 0 2 1 6 7 8 2 1 4 5 1 4 3 0 0 7 2 0 0)

columns  = (0   1   1   1   2)

pointerB = (0   2   3)

pointerE = (2   3   5)

This storage format is supported by the NIST Sparse BLAS library [Rem05].

Intel MKL supports the variation of the BSR format that is specified by three arrays: values, columns, and rowIndex. The following table describes these arrays.

values

A real array that contains the elements of the non-zero blocks of a sparse matrix. The elements are stored block by block in row-major order. A non-zero block is the block that contains at least one non-zero element. All elements of non-zero blocks are stored, even if some of them is equal to zero. Within each non-zero block the elements are stored in column major order in the case of the one-based indexing, and in row major order in the case of the zero-based indexing.

columns

Element i of the integer array columns is the number of the column in the block matrix that contains the i-th non-zero block.

rowIndex

Element j of this integer array gives the index of the element in the columns array that is first non-zero block in a row j of the block matrix.

The length of the values array is equal to the number of all elements in the non-zero blocks, the length of the columns array is equal to the number of non-zero blocks.

As the rowIndex array gives the location of the first non-zero block within a row, and the non-zero blocks are stored consecutively, the number of non-zero blocks in the i-th row is equal to the difference of rowIndex(i) and rowIndex(i+1).

To retain this relationship for the last row of the block matrix, an additional entry (dummy entry) is added to the end of rowIndex with value equal to the number of non-zeros blocks plus one. This makes the total length of the rowIndex array one larger than the number of rows of the block matrix.

The above matrix D can be represented in this 3-array variation of the BSR format as follows:

one-based indexing

values  =  (1 2 0 1 6 8 7 2 1 5 4 2 4 0 3 0 7 0 2 0)

columns  = (1   2   2   2   3)

rowIndex = (1   3   4   6)

zero-based indexing

values  =  (1 0 2 1 6 7 8 2 1 4 5 1 4 3 0 0 7 2 0 0)

columns  = (0   1   1   1   2)

rowIndex = (0   2   3 5)

When storing symmetric matrices, it is necessary to store only the upper or the lower triangular part of the matrix.

For example, consider the symmetric sparse matrix F:

If the size of the block equals 2, then the sparse matrix F can be represented as a 3x3 block matrix G with the following structure:

where

The symmetric matrix F can be represented in this 3-array variation of the BSR format (storing only upper triangular) as follows:

one-based indexing

values  =  (1 2 0 1 6 8 7 2 1 5 4 2 7 0 2 0)

columns  = (1   2   2   3)

rowIndex = (1   3   4 5)

zero-based indexing

values  =  (1 0 2 1 6 7 8 2 1 4 5 2 7 2 0 0)

columns  = (0   1   1   2)

rowIndex = (0   2   3 4)