Matrix Storage Schemes for LAPACK Routines

Developer Reference for Intel® oneAPI Math Kernel Library for C

Download PDF

ID 766684

Date 11/07/2023

Version

Public

A newer version of this document is available. Customers should click here to go to the newest version.

Visible to Intel only — GUID: GUID-BEC328BC-73C7-453F-B375-DB017AD24971

View Details

Matrix Storage Schemes for LAPACK Routines

LAPACK routines use the following matrix storage schemes:

Full Storage
Packed Storage
Band Storage
Rectangular Full Packed (RFP) Storage

Full Storage

Consider an m-by-n matrix A :

It is stored in a one-dimensional array a of length at least lda*n for column major layout or m*lda for row major layout. Element a_i,j is stored as array element a[k] where the mapping of k(i, j) is defined as

column major layout: k(i, j) = i - 1 + (j - 1)*lda
row major layout: k(i, j) = (i - 1)*lda + j - 1

NOTE:

Although LAPACK accepts parameter values of zero for matrix size, in general the size of the array used to store an m-by-n matrix A with leading dimension lda should be greater than or equal to max(1, n*lda) for column major layout and max (1, m*lda) for row major layout.

NOTE:

Even though the array used to store a matrix is one-dimensional, for simplicity the documentation sometimes refers parts of the array such as rows, columns, upper and lower triangular part, and diagonals. These refer to the parts of the matrix stored within the array. For example, the lower triangle of array a is defined as the subset of elements a[k(i,j)] with i≥j.

Packed Storage

The packed storage format compactly stores matrix elements when only one part of the matrix, the upper or lower triangle, is necessary to determine all of the elements of the matrix. This is the case when the matrix is upper triangular, lower triangular, symmetric, or Hermitian. For an n-by-n matrix of one of these types, a linear array ap of length n*(n + 1)/2 is adequate. Two parameters define the storage scheme: matrix_layout, which specifies column major (with the value LAPACK_COL_MAJOR) or row major (with the value LAPACK_ROW_MAJOR) matrix layout, and uplo, which specifies that the upper triangle (with the value 'U') or the lower triangle (with the value 'L') is stored.

Element a_i,j is stored as array element a[k] where the mapping of k(i, j) is defined as

		matrix_layout = LAPACK_COL_MAJOR	matrix_layout = LAPACK_ROW_MAJOR
uplo = `'U'`	1 ≤`i`≤`j`≤`n`	`k`(`i`, `j`) = `i` - 1 + `j`*(`j` - 1)/2	`k`(`i`, `j`) = `j` - 1 + (`i` - 1)(2`n` - `i`)/2
uplo = `'L'`	1 ≤`j`≤`i`≤`n`	`k`(`i`, `j`) = `i` - 1 + (`j` - 1)(2`n` - `j`)/2	`k`(`i`, `j`) = `j` - 1 + `i`*(`i` - 1)/2

NOTE:

Although LAPACK accepts parameter values of zero for matrix size, in general the size of the array should be greater than or equal to max(1, nx*(n + 1)/2).

Band Storage

When the non-zero elements of a matrix are confined to diagonal bands, it is possible to store the elements more efficiently using band storage. For example, consider an m-by-n band matrix A with kl subdiagonals and ku superdiagonals:

This matrix can be stored compactly in a one dimensional array ab. There are two operations involved in storing the matrix: packing the band matrix into matrix AB, and converting the packed matrix to a one-dimensional array.

Packing the Band Matrix: How the band matrix is packed depends on the matrix layout.
- Column major layout: matrix A is packed in an ldab-by-n matrix AB column-wise so that the diagonals of A become rows of array AB.
  
  The number of rows of ABldab≥kl + ku + 1, and the number of columns of AB is n.
- Row major layout: matrix A is packed in an m-by-ldab matrix AB row-wise so that the diagonals of A become columns of AB.
  
  The number of columns of ABldab≥kl + ku + 1, and the number of rows of AB is m.
NOTE:

For both column major and row major layout, elements of the upper left triangle of AB are not used. Depending on the relationship of the dimensions m, n, kl, and ku, the lower right triangle might not be used.
Converting the Packed Matrix to a One-Dimensional Array: The packed matrix AB is stored in a linear array ab as described in Full Storage . The size of ab should be greater than or equal to the total number of elements of matrix AB: ldab*n for column major layout or ldab*m for row major layout. The leading dimension of ab, ldab, must be greater than or equal to kl + ku + 1 (and some routines require it to be even larger).

Element a_i,j is stored as array element a[k(i, j)] where the mapping of k(i, j) is defined as
- column major layout: k(i, j) = i + ku - j + (j - 1)*ldab; 1 ≤j≤n, max(1, j - ku) ≤i≤ min(m, j + kl)
- row major layout: k(i,j) = j-i+k_l+(i-1)(k_l+k_u+1), 1 ≤ i ≤ m, max(1, i - k_l) ≤ j ≤ min(n, i + k_u)

NOTE:

Although LAPACK accepts parameter values of zero for matrix size, in general the size of the array should be greater than or equal to max(1, n*ldab) for column major layout and max (1, m*ldab) for row major layout.

Rectangular Full Packed Storage

A combination of full and packed storage, rectangular full packed storage can be used to store the upper or lower triangle of a matrix which is upper triangular, lower triangular, symmetric, or Hermitian. It offers the storage savings of packed storage plus the efficiency of using full storage Level 3 BLAS and LAPACK routines. Three parameters define the storage scheme: matrix_layout, which specifies column major (with the value LAPACK_COL_MAJOR) or row major (with the value LAPACK_ROW_MAJOR) matrix layout; uplo, which specifies that the upper triangle (with the value 'U') or the lower triangle (with the value 'L') is stored;and transr, which specifies normal (with the value 'N'), transpose (with the value 'T'), or conjugate transpose (with the value 'C') operation on the matrix.

Consider an N-by-N matrix A:

The upper or lower triangle of A can be stored in the array ap of length N*(N + 1)/2.

Additionally, define k as the integer part of N/2, such that N=2*k if N is even, and N=2*k + 1 if N is odd.

Storing the matrix involves packing the matrix into a rectangular matrix, and then storing the matrix in a one-dimensional array. The size of rectangular matrix AP required for the N-by-N matrix A is N + 1 by N/2 for even N, and N by (N + 1)/2 for odd N.

These examples illustrate the rectangular full packed storage method.

Upper triangular - uplo = 'U'

Consider a matrix A with N = 6:
- Not transposed - transr = 'N'
  
  The elements of the upper triangle of A can be packed in a matrix with the dimensions (N + 1)-by-(N/2) = 7 by 3:
- Transposed or conjugate transposed - transr = 'T' or transr = 'C'
  
  The elements of the upper triangle of A can be packed in a matrix with the dimensions (N/2) by (N + 1) = 3 by 7:
Consider a matrix A with N = 5:
- Not transposed - transr = 'N'
  
  The elements of the upper triangle of A can be packed in a matrix with the dimensions (N)-by-((N+1)/2) = 5 by 3:
- Transposed or conjugate transposed - transr = 'T' or transr = 'C'
  
  The elements of the upper triangle of A can be packed in a matrix with the dimensions ((N+1)/2) by (N ) = 5 by 3:
Lower triangular - uplo = 'L'

Consider a matrix A with N = 6:
- Not transposed - transr = 'N'
  
  The elements of the lower triangle of A can be packed in a matrix with the dimensions (N + 1)-by-(N/2) = 7 by 3:
- Transposed or conjugate transposed - transr = 'T' or transr = 'C'
  
  The elements of the lower triangle of A can be packed in a matrix with the dimensions (N/2) by (N + 1) = 3 by 7:
Consider a matrix A with N = 5:
- Not transposed - transr = 'N'
  
  The elements of the lower triangle of A can be packed in a matrix with the dimensions (N)-by-((N+1)/2) = 5 by 3:
- Transposed or conjugate transposed - transr = 'T' or transr = 'C'
  
  The elements of the lower triangle of A can be packed in a matrix with the dimensions ((N+1)/2) by (N ) = 5 by 3:

The packed matrix AP can be stored using column major layout or row major layout.

NOTE:

The matrix_layout and transr parameters can specify the same storage scheme: for example, the storage scheme for matrix_layout = LAPACK_COL_MAJOR and transr = 'N' is the same as that for matrix_layout = LAPACK_ROW_MAJOR and transr = 'T'.

Element a_i,j is stored as array element ap[l] where the mapping of l(i, j) is defined in the following tables.

Column major layout: matrix_layout = LAPACK_COL_MAJOR a

transr	uplo	`N`	`l(i, j) =`	`i`	`j`
'N'	'U'	`2*k`	`(j - k)*(N + 1) + i`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k`	`i*(N + 1) + j + k + 1`	`0 ≤i < k`	`i≤j < k`
		`2*k + 1`	`(j - k)*N + i`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k + 1`	`i*N + j + k + 1`	`0 ≤i < k`	`i≤j < k`
	'L'	`2*k`	`j*(N + 1) + i + 1`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k`	`(i - k)*(N + 1) + j - k`	`k≤i < N`	`k≤j≤i`
		`2*k + 1`	`j*N + i`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k + 1`	`(i - k)*N + j - k - 1`	`k + 1 ≤i < N`	`k + 1 ≤j≤i`
'T' or 'C'	'U'	`2*k`	`i*k + j - k`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k`	`(j + k + 1)*k + i`	`0 ≤i < k`	`i≤j < k`
		`2*k + 1`	`i*(k + 1) + j - k`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k + 1`	`(j + k + 1)*(k + 1) + i`	`0 ≤i < k`	`i≤j < k`
	'L'	`2*k`	`(i + 1)*k + j`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k`	`(j - k)*k + i - k`	`k≤i < N`	`k≤j≤i`
		`2*k + 1`	`i*(k + 1) + j`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k + 1`	`(j - k - 1)*(k + 1) + i - k`	`k + 1 ≤i < N`	`k + 1 ≤j≤i`

Row major layout: matrix_layout = LAPACK_ROW_MAJOR

transr	uplo	`N`	`l(i, j) =`	`i`	`j`
'N'	'U'	`2*k`	`i*k + j - k`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k`	`(k + j + 1)*k + i`	`0 ≤i < k`	`i≤j < k`
		`2*k + 1`	`i*(k + 1) + j - k`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k + 1`	`(k + j + 1)*(k + 1) + i`	`0 ≤i < k`	`i≤j < k`
	'L'	`2*k`	`(i + 1)*k + j`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k`	`(j - k)*k + i - k`	`k≤i < N`	`k≤j≤i`
		`2*k + 1`	`i*(k + 1) + j`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k + 1`	`(j - k - 1)*(k + 1) + i - k`	`k + 1 ≤i < N`	`k + 1 ≤j≤i`
'T' or 'C'	'U'	`2*k`	`(j - k)*(N + 1) + i`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k`	`i*(N + 1) + k + j + 1`	`0 ≤i < k`	`i≤j < k`
		`2*k + 1`	`(j - k)*N + i`	`0 ≤i < N`	`max(i, k) ≤j < N`
		`2*k + 1`	`i*N + k + j + 1`	`0 ≤i < k`	`i≤j < k`
	'L'	`2*k`	`j*(N + 1) + i + 1`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k`	`(i - k)*(N + 1) + j - k`	`k≤i < N`	`k≤j≤i`
		`2*k + 1`	`j*N + i`	`0 ≤i < N`	`0 ≤j≤ min(i, k)`
		`2*k + 1`	`(i - k)*N + j - k - 1`	`k + 1 ≤i < N`	`k + 1 ≤j≤i`

NOTE:

Although LAPACK accepts parameter values of zero for matrix size, in general the size of the array should be greater than or equal to max(1, N*(N + 1)/2).

Parent topic: LAPACK Routines

Select Your Language

Using Intel.com Search

Quick Links

Recent Searches

Advanced Search

Only search in

Developer Reference for Intel® oneAPI Math Kernel Library for C

Matrix Storage Schemes for LAPACK Routines

Full Storage

Packed Storage

Band Storage

Rectangular Full Packed Storage