Matrix Storage Modes
In this section, the word matrix is used to refer to a mathematical object and the word array is used to refer to its representation as a C data structure. In the following list of array types, the IMSL C Math Library functions require input consisting of matrix dimension values and all values for the matrix entries. These values are stored in row-major order in the arrays.
Each function processes the input array and typically returns a pointer to a “result.” For example, in solving linear algebraic systems, the pointer is to the solution. For general, real eigenvalue problems, the pointer is to the eigenvalues. Normally, the input array values are not changed by the functions.
In the IMSL C Math Library, an array is a pointer to a contiguous block of data. They are not pointers to pointers to the rows of the matrix. Typical declarations are:
float *a = {1, 2, 3, 4};
float b[2][2] = {1, 2, 3, 4};
float c[] = {1, 2, 3, 4};
General Mode
A general matrix is a square n × n matrix. The data type of a general array can be float, double, f_complex, or d_complex.
Rectangular Mode
A rectangular matrix is an m × n matrix. The data type of a rectangular array can be float, double, f_complex, or d_complex.
Symmetric Mode
A symmetric matrix is a square n × n matrix A, such that AT = A. (The matrix AT is the transpose of A.) The data type of a symmetric array can be float or double.
Hermitian Mode
A Hermitian matrix is a square n × n matrix A, such that
The matrix Ā is the complex conjugate of A, and
is the conjugate transpose of A. For Hermitian matrices AH = A. The data type of a Hermitian array can be f_complex or d_complex.
Sparse Coordinate Storage Format
Only the nonzero elements of a sparse matrix need to be communicated to a function. Sparse coordinate storage format stores the value of each matrix entry along with that entry’s row and column index. The following four non-homogeneous data structures are defined to support this concept:
typedef struct {
int row;
int col;
float val;
} Imsl_f_sparse_elem;
typedef struct {
int row;
int col;
double val;
} Imsl_d_sparse_elem;
typedef struct {
int row;
int col;
f_complex val;
} Imsl_c_sparse_elem;
typedef struct {
int row;
int col;
d_complex val;
} Imsl_z_sparse_elem;
See the Complex Data Types and Functions in the Reference Material at the end of this manual for a discussion of the complex data types f_complex and d_complex. Note that the only difference in these structures involves changes in underlying data types. A sparse matrix is passed to functions that accept sparse coordinate format by forming an array of one of these data types. The number of elements in that array will be equal to the number of nonzeros in the sparse matrix.
As an example consider the 6 × 6 matrix:
The matrix A has 15 nonzero elements, and the sparse coordinate representation would be
|
row |
0 |
1 |
1 |
1 |
2 |
3 |
3 |
3 |
4 |
4 |
4 |
4 |
5 |
5 |
5 |
|
col |
0 |
1 |
2 |
3 |
2 |
0 |
3 |
4 |
0 |
3 |
4 |
5 |
0 |
1 |
5 |
|
val |
2 |
9 |
-3 |
-1 |
5 |
-2 |
-7 |
-1 |
-1 |
-5 |
1 |
-3 |
-1 |
-2 |
6 |
Since this representation does not rely on order, an equivalent form would be
|
row |
5 |
4 |
3 |
0 |
5 |
1 |
2 |
1 |
4 |
3 |
1 |
4 |
3 |
5 |
4 |
|
col |
0 |
0 |
0 |
0 |
1 |
1 |
2 |
2 |
3 |
3 |
3 |
4 |
4 |
5 |
5 |
|
val |
-1 |
-1 |
-2 |
2 |
-2 |
9 |
5 |
-3 |
-5 |
-7 |
-1 |
1 |
-1 |
6 |
-3 |
There are different ways this data could be used to initialize an array of type, for example, Imsl_f_sparse_elem. Consider the following program fragment:
#include <imsl.h>
int main()
{
Imsl_f_sparse_elem a[] = {
{0, 0, 2.0},
{1, 1, 9.0},
{1, 2, -3.0},
{1, 3, -1.0},
{2, 2, 5.0},
{3, 0, -2.0},
{3, 3, -7.0},
{3, 4, -1.0},
{4, 0, -1.0},
{4, 3, -5.0},
{4, 4, 1.0},
{4, 5, -3.0},
{5, 0, -1.0},
{5, 1, -2.0},
{5, 5, 6.0} };
Imsl_f_sparse_elem b[15];
b[0].row = b[0].col = 0; b[0].val = 2.0;
b[1].row = b[1].col = 1; b[1].val = 9.0;
b[2].row = 1; b[2].col = 2; b[2].val = -3.0;
b[3].row = 1; b[3].col = 3; b[3].val = -1.0;
b[4].row = b[4].col = 2; b[4].val = 5.0;
b[5].row = 3; b[5].col = 0; b[5].val = -2.0;
b[6].row = b[6].col = 3; b[6].val = -7.0;
b[7].row = 3; b[7].col = 4; b[7].val = -1;
b[8].row = 4; b[8].col = 0; b[8].val = -1.0;
b[9].row = 4; b[9].col = 3; b[9].val = -5.0;
b[10].row = b[10].col = 4; b[10].val = 1.0;
b[11].row = 4; b[11].col = 5; b[11].val = -3.0;
b[12].row = 5; b[12].col = 0; b[12].val = -1.0;
b[13].row = 5; b[13] = 1; b[13].val = -2.0;
b[14].row = b[14].col = 5; b[14].val = 6.0;
}
Both a and b represent the sparse matrix A, and the functions in this module would produce identical results regardless of which identifier was sent through the argument list.
A sparse symmetric or Hermitian matrix is a special case, since it is only necessary to store the diagonal and either the upper or lower triangle. As an example, consider the 5 × 5 linear system:
The Hermitian and symmetric positive definite system solvers in this library expect the diagonal and lower triangle to be specified. The sparse coordinate form for the lower triangle is given by
|
row |
0 |
1 |
2 |
3 |
1 |
2 |
3 |
|
col |
0 |
1 |
2 |
3 |
0 |
1 |
2 |
|
val |
(4,0) |
(4,0) |
(4,0) |
(4,0) |
(1,1) |
(1,1) |
(1,1) |
As before, an equivalent form would be
|
row |
0 |
1 |
1 |
2 |
2 |
3 |
(4,0) |
|
col |
0 |
0 |
1 |
1 |
2 |
2 |
3 |
|
val |
(4,0) |
(1,1) |
(4,0) |
(1,1) |
(4,0) |
(1,1) |
(4,0) |
The following program fragment will initialize both a and b to H.
#include <imsl.h>
int main()
{
Imsl_c_sparse_elem a[] = {
{0, 0, {4.0, 0.0}},
{1, 1, {4.0, 0.0}},
{2, 2, {4.0, 0.0}},
{3, 3, {4.0, 0.0}},
{1, 0, {1.0, 1.0}},
{2, 1, {1.0, 1.0}},
{3, 2, {1.0, 1.0}}
}
Imsl_c_sparse_elem b[7];
b[0].row = b[0].col = 0;
b[0].val = imsl_cf_convert (4.0, 0.0);
b[1].row = 1; b[1].col = 0;
b[1].val = imsl_cf_convert (1.0, 1.0);
b[2].row = b[2].col = 1;
b[2].val = imsl_cf_convert (4.0, 0.0);
b[3].row = 2; b[3].col = 1;
b[3].val = imsl_cf_convert (1.0, 1.0);
b[4].row = b[4].col = 2;
b[4].val = imsl_cf_convert (4.0, 0.0);
b[5].row = 3; b[5].col = 2;
b[5].val = imsl_cf_convert (1.0, 1.0);
b[6].row = b[6].col = 3;
b[6].val = imsl_cf_convert (4.0, 0.0);
}
There are some important points to note here. H is not symmetric, but rather Hermitian. The functions that accept Hermitian data understand this and operate assuming that
The IMSL C Math Library cannot take advantage of the symmetry in matrices that are not positive definite. The implication here is that a symmetric matrix that happens to be indefinite cannot be stored in this compact symmetric form. Rather, both upper and lower triangles must be specified and the sparse general solver called.
Band Storage Format
A band matrix is an M × N matrix with all of its nonzero elements “close” to the main diagonal. Specifically, values Aij = 0 if i - j > nlca or j – i > nuca. The integer m = nlca + nuca + 1 is the total band width. The diagonals, other than the main diagonal, are called codiagonals. While any M × N matrix is a band matrix, band storage format is only useful when the number of nonzero codiagonals is much less than N.
In band storage format, the nlca lower codiagonals and the nuca upper codiagonals are stored in the rows of an array of size M × N. The elements are stored in the same column of the array as they are in the matrix. The values Aij inside the band width are stored in the linear array in positions [(i ‑ j + nuca + 1) * n + j]. This results in a row-major, one-dimensional mapping from the two-dimensional notion of the matrix.
For example, consider the 5 × 5 matrix A with 1 lower and 2 upper codiagonals:
In band storage format, the data would be arranged as
This data would then be stored contiguously, row-major order, in an array of length 20.
As an example, consider the following tridiagonal matrix:
The following declaration will store this matrix in band storage format:
float a[] = {
0.0, 1.0, 2.0, 3.0, 4.0,
10.0, 20.0, 30.0, 40.0, 50.0,
5.0, 6.0, 7.0, 8.0, 0.0
};
As in the sparse coordinate representation, there is a space saving symmetric version of band storage. As an example, look at the following 5 × 5 symmetric problem:
In band symmetric storage format, the data would be arranged as
The following Hermitian example illustrates the procedure:
The following program fragments would store H in h, using band symmetric storage format.
f_complex h[] = {
{0.0, 0.0}, {0.0, 0.0}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0},
{0.0, 0.0}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0}, {1.0, 1.0},
{8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}, {8.0, 0.0}};
or equivalently
f_complex h[15];
h[0] = h[1] = h[5] = imsl_cf_convert (0.0, 0.0);
h[2] = h[3] = h[4] = h[6] = h[7] = h[8] = h[9] =
imsl_cf_convert (1.0, 1.0);
h[10] = h[11] = h[12] = h[13] = h[14] =
imsl_cf_convert (8.0, 0.0);
Choosing Between Banded and Coordinate Forms
It is clear that any matrix can be stored in either sparse coordinate or band format. The choice depends on the sparsity pattern of the matrix. A matrix with all nonzero data stored in bands close to the main diagonal would probably be a good candidate for band format. If nonzero information is scattered more or less uniformly through the matrix, sparse coordinate format is the best choice. As extreme examples, consider the following two cases: (1) an n × n matrix with all elements on the main diagonal and the (0, n - 1) and (n - 1, 0) entries nonzero. The sparse coordinate vector would be n + 2 units long. An array of length n(2n - 1) would be required to store the band representation, nearly twice as much storage as a dense solver might require. (2) a tridiagonal matrix with all diagonal, superdiagonal and subdiagonal entries nonzero. In band format, an array of length 3n is needed. In sparse coordinate, format a vector of length 3n - 2 is required. But the problem is that, for example, for float precision, each of those 3n - 2 units in coordinate format requires three times as much storage as any of the 3n units needed for band representation. This is due to carrying the row and column indices in coordinate form. Band storage evades this requirement by being essentially an ordered list, and defining location in the original matrix by position in the list.
Compressed Sparse Column (CSC) Format
Functions that accept data in coordinate format can also accept data stored in the format described in the Users’ Guide for the Harwell-Boeing Sparse Matrix Collection (via optional argument IMSL_CSC_FORMAT_none). The scheme is column oriented, with each column held as a sparse vector, represented by a list of the row indices of the entries in an integer array (“rowind” below) and a list of the corresponding values in a separate float (double, f_complex, d_complex) array (“values” below). Data for each column are stored consecutively and the columns are stored in order. A third array (“colptr” below) indicates the location in array “values” in which to place the first nonzero value of each succeeding column of the original sparse matrix. So colptr[i] contains the index of the first free location in array “values” in which to place the values from the ith column of the original sparse matrix. In other words, values[colptr[i]] holds the first nonzero value of the i-th column of the original sparse matrix. Only entries in the lower triangle and diagonal are stored for symmetric and Hermitian matrices. All arrays are based at zero, which is in contrast to the Harwell-Boeing test suite’s one-based arrays.
As in the Harwell-Boeing user guide (link above), the storage scheme is illustrated with the following example: The 5 × 5 matrix
would be stored in the arrays colptr (location of first entry), rowind (row indices), and values (nonzero entries) as follows:
|
Subscripts |
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
|
Colptr |
0 |
3 |
5 |
7 |
9 |
11 |
|
|
|
|
|
|
Rowind |
0 |
2 |
4 |
0 |
3 |
1 |
4 |
0 |
3 |
1 |
4 |
|
Values |
1 |
2 |
5 |
-3 |
4 |
-2 |
-5 |
-1 |
-4 |
3 |
6 |
The following program fragment shows the relation between CSC storage format and coordinate representation:
int main() {
int i, j, k, n = 5, nz, start, stop;
int colptr[] = { 0, 3, 5, 7, 9, 11 };
int rowind[] = { 0, 2, 4, 0, 3, 1, 4, 0, 3, 1, 4 };
float values[] = { 1.0, 2.0, 5.0, -3.0, 4.0, -2.0,
-5.0, -1.0, -4.0, 3.0, 6.0 };
Imsl_f_sparse_elem a[11];
k = 0;
for (i = 0; i < n; i++) {
start = colptr[i];
stop = colptr[i + 1];
for (j = start; j < stop; j++) {
a[k].row = rowind[j];
a[k].col = i;
a[k++].val = values[j];
}
}
nz = k;
}
