Matrix Storage Modes
In this section, the word matrix will be used to refer to a mathematical object, and the word array will be used to refer to its representation as a FORTRAN data structure.
General Mode
A general matrix is an N × N matrix A. It is stored in a FORTRAN array that is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as N. IMSL general matrix subprograms only refer to values Aij for i = 1, …, N and j = 1, …, N. The data type of a general array can be one of REAL, DOUBLE PRECISION, or COMPLEX. If your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX can also be declared.
Rectangular Mode
A rectangular matrix is an M × N matrix A. It is stored in a FORTRAN array that is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as M. IMSL rectangular matrix subprograms only refer to values Aij for i = 1, …, M and j = 1, …, N. The data type of a rectangular array can be REAL, DOUBLE PRECISION, or COMPLEX. If your FORTRAN compiler allows, you can declare the nonstandard data type DOUBLE COMPLEX.
Symmetric Mode
A symmetric matrix is a square N × N matrix A, such that AT = A. (AT is the transpose of A.) It is stored in a FORTRAN array that is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as N. IMSL symmetric matrix subprograms only refer to the upper or to the lower half of A (i.e., to values Aij for i = 1, …, N and j = i, …, N, or Aij for j = 1, …, N and i = j, …, N). The data type of a symmetric array can be one of REAL or DOUBLE PRECISION. Use of the upper half of the array is denoted in the BLAS that compute with symmetric matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library using the CHARACTER*1 flag UPLO = ’U’. Otherwise, UPLO = ’L’ denotes that the lower half of the array is used.
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 matrix is stored in a FORTRAN array that is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as N. IMSL Hermitian matrix subprograms only refer to the upper or to the lower half of A (i.e., to values Aij for i = 1, …, N and j = i, …, N, or Aij for j = 1, …, N and i = j, …, N). Use of the upper half of the array is denoted in the BLAS that compute with Hermitian matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library, using the CHARACTER*1 flag UPLO = ’U’. Otherwise, UPLO = ’L’ denotes that the lower half of the array is used. The data type of a Hermitian array can be COMPLEX or, if your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX.
Triangular Mode
A triangular matrix is a square N × N matrix A such that values Aij = 0 for i < j or Aij = 0 for i > j. The first condition defines a lower triangular matrix while the second condition defines an upper triangular matrix. A lower triangular matrix A is stored in the lower triangular part of a FORTRAN array A. An upper triangular matrix is stored in the upper triangular part of a FORTRAN array. Triangular matrices are called unit triangular whenever Ajj = 1, j = 1, …, N. For unit triangular matrices, only the strictly lower or upper parts of the array are referenced. This is denoted in the BLAS that compute with triangular matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library, using the CHARACTER*1 flag DIAG = ’U’. Otherwise, DIAG = ’N’ denotes that the diagonal array terms should be used. For unit triangular matrices, the diagonal terms are each used with the mathematical value 1. The array diagonal term does not need to be 1.0 in this usage. Use of the upper half of the array is denoted in the BLAS that compute with triangular matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library, sing the CHARACTER*1 flag UPLO = ’U’. Otherwise, UPLO = ’L’ denotes that the lower half of the array is used. The data type of an array that contains a triangular matrix can be one of REAL, DOUBLE PRECISION, or COMPLEX. If your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX can also be declared.
Band Storage Mode
A band matrix is an M × N matrix A with all of its nonzero elements “close” to the main diagonal. Specifically, values Aij = 0 if i ‑ j > NLCA or j ‑ i > NUCA. The integers NLCA and NUCA are the lower and upper band widths. 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, the band matrix mode is most useful only when the number of nonzero codiagonals is much less than m.
In the band storage mode, the NLCA lower codiagonals and NUCA upper codiagonals are stored in the rows of a FORTRAN array of dimension 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 array positions (i – j + NUCA + 1, j). This array is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as m. The data type of a band matrix array can be one of REAL, DOUBLE PRECISION, COMPLEX or, if your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX. Use of the CHARACTER*1 flag TRANS=’N’ in the BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library, specifies that the matrix A is used. The flag value
while
For example, consider a real 5 × 5 band matrix with 1 lower and 2 upper codiagonals, stored in the FORTRAN array declared by the following statements:
PARAMETER (N=5, NLCA=1, NUCA=2)
REAL A(NLCA+NUCA+1, N)
The matrix A has the form
As a FORTRAN array, it is
The entries marked with an x in the above array are not referenced by the IMSL band subprograms.
Band Symmetric Storage Mode
A band symmetric matrix is a band matrix that is also symmetric. The band symmetric storage mode is similar to the band mode except only the lower or upper codiagonals are stored.
In the band symmetric storage mode, the NCODA upper codiagonals are stored in the rows of a FORTRAN array of dimension (NCODA + 1) × N. The elements are stored in the same column of the array as they are in the matrix. Specifically, values Aij, j ≤ i inside the band are stored in array positions i ‑ j + NCODA + 1, j). This is the storage mode designated by using the CHARACTER*1 flag UPLO = ’U’ in Level 2 BLAS that compute with band symmetric matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library. Alternatively, Aij, j ≤ i, inside the band, are stored in array positions (i ‑ j + 1, j). This is the storage mode designated by using the CHARACTER*1 flag UPLO = ’L’ in these Level 2 BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library. The array is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as NCODA + 1. The data type of a band symmetric array can be REAL or DOUBLE PRECISION.
For example, consider a real 5 × 5 band matrix with 2 codiagonals. Its FORTRAN declaration is
PARAMETER (N=5, NCODA=2)
REAL A(NCODA+1, N)
The matrix A has the form
Since A is symmetric, the values Aij = Aji. In the FORTRAN array, it is
The entries marked with an × in the above array are not referenced by the IMSL band symmetric subprograms.
An alternate storage mode for band symmetric matrices is designated using the CHARACTER*1 flag UPLO = ’L’ in Level 2 BLAS that compute with band symmetric matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library. In that case, the example matrix is represented as
Band Hermitian Storage Mode
A band Hermitian matrix is a band matrix that is also Hermitian. The band Hermitian mode is a complex analogue of the band symmetric mode.
In the band Hermitian storage mode, the NCODA upper codiagonals are stored in the rows of a FORTRAN array of dimension (NCODA + 1) × N. The elements are stored in the same column of the array as they are in the matrix. In the Level 2 BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library, this is denoted by using the CHARACTER*1 flag UPLO =’U’. The array is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as
(NCODA + 1). The data type of a band Hermitian array can be COMPLEX or, if your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX.
For example, consider a complex 5 × 5 band matrix with 2 codiagonals. Its FORTRAN declaration is
PARAMETER (N=5, NCODA = 2)
COMPLEX A(NCODA + 1, N)
The matrix A has the form
where the value
is the complex conjugate of Aij. This matrix represented as a FORTRAN array is
The entries marked with an × in the above array are not referenced by the IMSL band Hermitian subprograms.
An alternate storage mode for band Hermitian matrices is designated using the CHARACTER*1 flag UPLO = ’L’ in Level 2 BLAS that compute with band Hermitian matrices, see Chapter 9, “Programming Notes for BLAS” in the Math Library. In that case, the example matrix is represented as
Band Triangular Storage Mode
A band triangular matrix is a band matrix that is also triangular. In the band triangular storage mode, the NCODA codiagonals are stored in the rows of a FORTRAN array of dimension (NCODA + 1) × N. The elements are stored in the same column of the array as they are in the matrix. For usage in the Level 2 BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library, the CHARACTER*1 flag DIAG has the same meaning as used in section “Triangular Storage Mode”. The flag UPLO has the meaning analogous with its usage in the section “Banded Symmetric Storage Mode”. This array is declared by the following statement:
DIMENSION A(LDA,N)
The parameter LDA is called the leading dimension of A. It must be at least as large as (NCODA + 1).
For example, consider a 5 × 5 band upper triangular matrix with 2 codiagonals. Its FORTRAN declaration is
PARAMETER (N = 5, NCODA = 2)
COMPLEX A(NCODA + 1, N)
The matrix A has the form
This matrix represented as a FORTRAN array is
This corresponds to the CHARACTER*1 flags DIAG = ’N’ and UPLO = ’U’. The matrix AT is represented as the FORTRAN array
This corresponds to the CHARACTER*1 flags DIAG = ’N’ and UPLO = ’L’. In both examples, the entries indicated with an × are not referenced by IMSL subprograms.
Codiagonal Band Symmetric Storage Mode
This is an alternate storage mode for band symmetric matrices. It is not used by any of the BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library. Storing data in a form transposed from the Band Symmetric Storage Mode maintains unit spacing between consecutive referenced array elements. This data structure is used to get good performance in the Cholesky decomposition algorithm that solves positive definite symmetric systems of linear equations Ax = b. The data type can be REAL or DOUBLE PRECISION. In the codiagonal band symmetric storage mode, the NCODA upper codiagonals and right‑hand‑side are stored in columns of this FORTRAN array. This array is declared by the following statement:
DIMENSION A(LDA, NCODA + 2)
The parameter LDA is the leading positive dimension of A. It must be at least as large as N + NCODA.
Consider a real symmetric 5 × 5 matrix with 2 codiagonals
and a right‑hand‑side vector
A FORTRAN declaration for the array to hold this matrix and right‑hand‑side vector is
PARAMETER (N = 5, NCODA = 2, LDA = N + NCODA)
REAL A(LDA, NCODA + 2)
The matrix and right‑hand‑side entries are placed in the FORTRAN array A as follows:
Entries marked with an × do not need to be defined. Certain of the IMSL band symmetric subprograms will initialize and use these values during the solution process. When a solution is computed, the bi, i = 1, …, 5, are replaced by xi, i = 1, …, 5.
The nonzero Aij, j ≥ i, are stored in array locations A(j + NCODA, (j ‑ i) + 1) . The right‑hand‑side entries bj are stored in locations A(j + NCODA, NCODA + 2). The solution entries xj are returned in A(j + NCODA, NCODA + 2).
Codiagonal Band Hermitian Storage Mode
This is an alternate storage mode for band Hermitian matrices. It is not used by any of the BLAS, see Chapter 9, “Programming Notes for BLAS” in the Math Library. In the codiagonal band Hermitian storage mode, the real and imaginary parts of the 2 * NCODA + 1 upper codiagonals and right‑hand‑side are stored in columns of a FORTRAN array. Note that there is no explicit use of the COMPLEX or the nonstandard data type DOUBLE COMPLEX data type in this storage mode.
For Hermitian complex matrices,
where U and V are real matrices. They satisfy the conditions U = UT and V = ‑VT. The right‑hand‑side
where c and d are real vectors. The solution vector is denoted as
where u and v are real. The storage is declared with the following statement
DIMENSION A(LDA, 2*NCODA + 3)
The parameter LDA is the leading positive dimension of A. It must be at least as large as N + NCODA.
The diagonal terms Ujj are stored in array locations A (j + NCODA, 1). The diagonal Vjj are zero and are not stored. The nonzero Uij, j > i, are stored in locations A(j + NCODA, 2 * (j ‑ i)).
The nonzero Vij are stored in locations A(j + NCODA, 2*(j ‑ i) + 1). The right side vector b is stored with cj and dj in locations A(j + NCODA, 2*NCODA + 2) and A(j + NCODA, 2*NCODA + 3) respectively. The real and imaginary parts of the solution, uj and vj, respectively overwrite cj and dj.
Consider a complex hermitian 5 × 5 matrix with 2 codiagonals
and a right‑hand‑side vector
A FORTRAN declaration for the array to hold this matrix and right‑hand‑side vector is
PARAMETER (N = 5, NCODA = 2, LDA = N + NCODA)
REAL A(LDA,2*NCODA + 3)
The matrix and right‑hand‑side entries are placed in the FORTRAN array A as follows:
Entries marked with an × do not need to be defined.
Sparse Matrix Storage Mode
The sparse linear algebraic equation solvers in Chapter 1, “Basic Statistics” accept the input matrix in sparse storage mode. This structure consists of INTEGER values N and NZ, the matrix dimension and the total number of nonzero entries in the matrix. In addition, there are two INTEGER arrays IROW(*) and JCOL(*) that contain unique matrix row and column coordinates where values are given. There is also an array A(*) of values. All other entries of the matrix are zero. Each of the arrays IROW(*), JCOL(*), A(*) must be of size NZ. The correspondence between matrix and array entries is given by
AIROW(i),JCOL(i) = A(i), i = 1, …, NZ
The data type for A(*) can be one of REAL, DOUBLE PRECISION, or COMPLEX. If your FORTRAN compiler allows, the nonstandard data type DOUBLE COMPLEX can also be declared.
For example, consider a real 5 × 5 sparse matrix with 11 nonzero entries. The matrix A has the form
Declarations of arrays and definitions of the values for this sparse matrix are
PARAMETER (NZ = 11, N = 5)
DIMENSION IROW(NZ), JCOL(NZ), A(NZ)
DATA IROW /1,1,1,2,2,3,3,3,4,5,5/
DATA JCOL /1,3,4,1,2,2,3,4,3,4,5/
DATA A /A11,A13,A14,A21,A22,A32,A33,A34, &
A43,A54,A55/
