docs: add band matrix storage format documentation in the README of lapack/base

[aayush0325]

Description

What is the purpose of this pull request?

This pull request:

• adds band matrix storage format documentation in the README of lapack/base as discussed

Related Issues

Does this pull request have any related issues?

No.

Questions

Any questions for reviewers of this pull request?

No.

Other

Any other information relevant to this pull request? This may include screenshots, references, and/or implementation notes.

No.

Checklist

Please ensure the following tasks are completed before submitting this pull request.

• Read, understood, and followed the contributing guidelines.

---

@stdlib-js/reviewers

[stdlib-bot]

Coverage Report

Package	Statements	Branches	Functions	Lines
lapack/base	$\color{red}9786/12473$ $\color{red}-21.54%$	$\color{red}130/155$ $\color{red}-16.13%$	$\color{red}0/100$ $\color{red}-100.00%$	$\color{red}9786/12473$ $\color{red}-21.54%$

The above coverage report was generated for the changes in this PR.

[kgryte]

It would be good to actually show matrices which are upper and lower triangular here.

[kgryte]

Overall, this is looking good. It's a more readable version of https://netlib.org/lapack/lug/node124.html. That stated, one thing which needs to be changed is that currently everything is written using one-based indexing. We need to update everything to zero-based indexing.

[kgryte]

This comment still needs to be resolved.

[kgryte]

One question I have: how do row- and column-major orders factor in here. For example, in dgbbrd, the routine accepts a layout argument.

It may be good to also show how the band storage matrices are stored in linear memory for the two supported layouts. IIRC, this was a point of confusion during development.

[aayush0325]

One question I have: how do row- and column-major orders factor in here. For example, in dgbbrd, the routine accepts a layout argument.

this is how i understand the row and column major interpretations of these matrices:

for column major:

{
  "order": "column-major",
  "type": "banded",
  "M": 5,
  "N": 5,
  "KL": 2,
  "KU": 1,
  "A": [
    0.0, 0.0, 0.0, 1.1, 2.1, 3.1,
    0.0, 0.0, 1.2, 2.2, 3.2, 4.2,
    0.0, 0.0, 2.3, 3.3, 4.3, 5.3,
    0.0, 0.0, 3.4, 4.4, 5.4, 0.0,
    0.0, 0.0, 4.5, 5.5, 0.0, 0.0
  ],
  "A_mat": [
    [ 1.1, 1.2, 0.0, 0.0, 0.0 ],
    [ 2.1, 2.2, 2.3, 0.0, 0.0 ],
    [ 3.1, 3.2, 3.3, 3.4, 0.0 ],
    [ 0.0, 4.2, 4.3, 4.4, 4.5 ],
    [ 0.0, 0.0, 5.3, 5.4, 5.5 ]
  ],
}

for row major:

{
  "order": "row-major",
  "type": "banded",
  "M": 5,
  "N": 5,
  "KL": 2,
  "KU": 1,
  "A": [
    0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.0, 0.0,
    0.0, 1.2, 2.3, 3.4, 4.5,
    1.1, 2.2, 3.3, 4.4, 5.5,
    2.1, 3.2, 4.3, 5.4, 0.0,
    3.1, 4.2, 5.3, 0.0, 0.0
  ],
  "A_mat": [
    [ 1.1, 1.2, 0.0, 0.0, 0.0 ],
    [ 2.1, 2.2, 2.3, 0.0, 0.0 ],
    [ 3.1, 3.2, 3.3, 3.4, 0.0 ],
    [ 0.0, 4.2, 4.3, 4.4, 4.5 ],
    [ 0.0, 0.0, 5.3, 5.4, 5.5 ]
  ]
}

i am basically reading this column-by-column for col-major and row-by-row for row-major

[kgryte]

@aayush0325 What you are displaying is specific for LU correct with the extra cols/rows?

Regardless, I suggest making explicit with a couple of brief examples showing how banded matrices are stored according to layouts.

[kgryte]

A few comments still need to be resolved before we can finally get this one in. :)

[kgryte]

cc @DivitJain26 Would you mind taking a look at this? Cheers!

[DivitJain26]

hello

MKL describes two methods of packing a band matrix according to the layout in Matrix Storage Schemes for LAPACK Routines.

A more descriptive explanation can be found in the cblas_?gbmv

Layout = CblasColMajor: Before entry, the leading (kl + ku + 1) by n part of the array a must contain the matrix of coefficients. This matrix must be supplied column-by-column, with the leading diagonal of the matrix in row (ku) of the array, the first super-diagonal starting at position 1 in row (ku - 1), the first sub-diagonal starting at position 0 in row (ku + 1), and so on. Elements in the array a that do not correspond to elements in the band matrix (such as the top left ku by ku triangle) are not referenced.

Layout = CblasRowMajor: Before entry, the leading (kl + ku + 1) by m part of the array a must contain the matrix of coefficients. This matrix must be supplied row-by-row, with the leading diagonal of the matrix in column (kl) of the array, the first super-diagonal starting at position 0 in column (kl + 1), the first sub-diagonal starting at position 1 in row (kl - 1), and so on. Elements in the array a that do not correspond to elements in the band matrix (such as the top left kl by kl triangle) are not referenced.

This is true for BLAS routines as I have demonstrated here Layout Dependent Band Matrix Packing in cblas_?gbmv.

However, it does not hold for LAPACK routines: routines that operate on band matrices expect the input matrix to be packed as defined for Layout = CblasColMajor, regardless of the chosen layout, as demonstrated LAPACKE band routines always expect column-major style band packing; even when layout = LAPACK_ROW_MAJOR.

• For layout = LAPACK_COL_MAJOR, this compact matrix should be supplied column-wise.
• For layout = LAPACK_ROW_MAJOR, the compact matrix should be supplied row-wise.

This contradicts Matrix Storage Schemes for LAPACK Routines and warrants further investigation.

The reason for this behaviour becomes clear when tracing the call chain. When dgbbrd is called, its wrapper lapacke_dgbbrd.c validates the inputs, and calls lapacke_dgbbrd_work.c which in turn calls the native implementation dgbbrd.c. The native implementation has no layout parameter, as it operates exclusively in column-major style. Consequently, the wrapper simply transposes the input matrix when layout = LAPACK_ROW_MAJOR. This is why, even in the row-major case, the matrix must be packed according to the column-major convention and then supplied row-wise.

As mentioned in Matrix Layout (Row Major and Column Major).

Matrix Layout (Row Major and Column Major)
The standard BLAS and LAPACK APIs require a Fortran layout for matrices (column major), where matrices are stored column-by-column in memory and the entries in each column are stored in consecutive memory locations. By default, oneMKL for DPC++ likewise assumes this matrix layout. The oneapi::mkl::blas::row_major namespace must be used for row major layout for BLAS. Row major layout is not supported directly for LAPACK, but you can transpose row major input matrices, call the desired DPC++ LAPACK routines, and then transpose the output matrices back to row major layout.

References:
[1] Matrix Storage Schemes for LAPACK Routines.
[2] cblas_?gbmv
[3] Matrix Layout (Row Major and Column Major).
[4] writing a banded matrix in a row major layout for lapack solver dgbsv

[DivitJain26]

LAPACKE band routines always expect column-major style band packing; even when layout = LAPACK_ROW_MAJOR

Demonstrates that LAPACKE band-matrix routines always expect the input packed using the column-major convention, regardless of the layout parameter.

#include <stdio.h>
#include <lapacke.h>

int main() {

	lapack_int M = 4;
	lapack_int N = 4;

	lapack_int KL = 2;
	lapack_int KU = 1;

	/*
		A_mat: [
			[ 1.0,  2.0,  0.0,  0.0 ],
			[ 3.0,  4.0,  5.0,  0.0 ],
			[ 6.0,  7.0,  8.0,  9.0 ],
			[ 0.0, 10.0, 11.0, 12.0 ]
		]
	*/

	/*
		Compact band storage (column-major convention):
		A_compact: [
			[ 0.0,  2.0,  5.0,  9.0 ],
			[ 1.0,  4.0,  8.0, 12.0 ],
			[ 3.0,  7.0, 11.0,  0.0 ],
			[ 6.0, 10.0,  0.0,  0.0 ]
		]
	*/

	lapack_int LDAB = KL + KU + 1;

	double AB_col[4*4] = { 0.0, 1.0, 3.0, 6.0, 2.0, 4.0, 7.0, 10.0, 5.0, 8.0, 11.0, 0.0, 9.0, 12.0, 0.0, 0.0 };
	double AB_row[4*4] = { 0.0, 2.0, 5.0, 9.0, 1.0, 4.0, 8.0, 12.0, 3.0, 7.0, 11.0, 0.0, 6.0, 10.0, 0.0, 0.0 };

	double D_col[4];
	double E_col[4];
	double Q_col[4*4];
	double PT_col[4*4];
	double C_col[4*4];
	lapack_int info_col;
	double D_row[4];
	double E_row[4];
	double Q_row[4*4];
	double PT_row[4*4];
	double C_row[4*4];
	lapack_int info_row;

	info_col = LAPACKE_dgbbrd( LAPACK_COL_MAJOR, 'B', M, N, 0, KL, KU, AB_col, LDAB, D_col, E_col, Q_col, 4, PT_col, 4, C_col, 4 );
	info_row = LAPACKE_dgbbrd( LAPACK_ROW_MAJOR, 'B', M, N, 0, KL, KU, AB_row, LDAB, D_row, E_row, Q_row, 4, PT_row, 4, C_row, 4 );

	printf("\n=======Column Major========\n");
	printf("INFO = %d\n", info_col);

	printf("\nDiagonal D:\n");
	for(int i = 0; i < 4; i++) printf("%8.3f ", D_col[i]);

	printf("\n\nOff-diagonal E:\n");
	for(int i = 0; i < 3; i++) printf("%8.3f ", E_col[i]);

	printf("\n\nQ matrix:\n");
	for(int i = 0; i < 4; i++) {
		for(int j = 0; j < 4; j++) {
			printf("%8.3f ", Q_col[i + j*4]);
		}
		printf("\n");
	}

	printf("\nPT matrix:\n");
	for(int i = 0; i < 4; i++) {
		for(int j = 0; j < 4; j++) {
			printf("%8.3f ", PT_col[i + j*4]);
		}
		printf("\n");
	}

	printf("\n=======Row Major========\n");
	printf("INFO = %d\n", info_row);

	printf("\nDiagonal D:\n");
	for(int i = 0; i < 4; i++) printf("%8.3f ", D_row[i]);

	printf("\n\nOff-diagonal E:\n");
	for(int i = 0; i < 3; i++) printf("%8.3f ", E_row[i]);

	printf("\n\nQ matrix:\n");
	for(int i = 0; i < 4; i++) {
		for(int j = 0; j < 4; j++) {
			printf("%8.3f ", Q_row[i*4 + j]);
		}
		printf("\n");
	}

	printf("\nPT matrix:\n");
	for(int i = 0; i < 4; i++) {
		for(int j = 0; j < 4; j++) {
			printf("%8.3f ", PT_row[i*4 + j]);
		}
		printf("\n");
	}

	return 0;
}

Compile

gcc file_name.c -lmkl_rt -llapacke -lpthread -lm -ldl

Output

=======Column Major========
INFO = 0

Diagonal D:
   6.782   19.069    3.795    1.540 

Off-diagonal E:
  14.760    2.398   -0.123 

Q matrix:
   0.147   -0.055   -0.137    0.978 
   0.442   -0.060   -0.874   -0.193 
   0.885    0.039    0.460   -0.067 
   0.000    0.996   -0.078    0.045 

PT matrix:
   1.000    0.000    0.000    0.000 
   0.000    0.559    0.629    0.539 
   0.000   -0.327   -0.430    0.841 
   0.000    0.762   -0.647   -0.035 

=======Row Major========
INFO = 0

Diagonal D:
   6.782   19.069    3.795    1.540 

Off-diagonal E:
  14.760    2.398   -0.123 

Q matrix:
   0.147   -0.055   -0.137    0.978 
   0.442   -0.060   -0.874   -0.193 
   0.885    0.039    0.460   -0.067 
   0.000    0.996   -0.078    0.045 

PT matrix:
   1.000    0.000    0.000    0.000 
   0.000    0.559    0.629    0.539 
   0.000   -0.327   -0.430    0.841 
   0.000    0.762   -0.647   -0.035