feat: add blas/base/ctbsv

lib/node_modules/@stdlib/blas/base/ctbsv/lib/base.js

@@ -0,0 +1,259 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2026 The Stdlib Authors.
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*    http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+'use strict';
+
+// MODULES //
+
+var reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
+var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
+var f32 = require( '@stdlib/number/float64/base/to-float32' );
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+
+
+// MAIN //
+
+/**
+* Solves one of the systems of equations `A*x = b` or `A^T*x = b` or `x = A^H*x` where `b` and `x` are `N` element complex vectors and `A` is an `N` by `N` unit, or non-unit, upper or lower triangular complex band matrix, with ( `K` + 1 ) diagonals.
+*
+* @private
+* @param {string} uplo - specifies whether `A` is an upper or lower triangular matrix
+* @param {string} trans - specifies whether `A` should be transposed, conjugate-transposed, or not transposed
+* @param {string} diag - specifies whether `A` has a unit diagonal
+* @param {NonNegativeInteger} N - number of elements along each dimension of `A`
+* @param {NonNegativeInteger} K - number of super-diagonals or sub-diagonals of the matrix `A`
+* @param {Complex64Array} A - input complex matrix
+* @param {integer} strideA1 - stride of the first dimension of `A`
+* @param {integer} strideA2 - stride of the second dimension of `A`
+* @param {NonNegativeInteger} offsetA - starting index for `A`
+* @param {Complex64Array} x - input complex vector
+* @param {integer} strideX - stride length for `x`
+* @param {NonNegativeInteger} offsetX - starting index for `x`
+* @returns {Complex64Array} `x`
+*
+* @example
+* var Complex64Array = require( '@stdlib/array/complex64' );
+*
+* var A = new Complex64Array( [ 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0 ] );
+* var x = new Complex64Array( [ 0.0, 2.0, 0.0, 16.0, 0.0, 46.0 ] );
+*
+* ctbsv( 'lower', 'no-transpose', 'non-unit', 3, 1, A, 2, 1, 0, x, 1, 0 );
+* // x => <Complex64Array>[ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+*/
+function ctbsv( uplo, trans, diag, N, K, A, strideA1, strideA2, offsetA, x, strideX, offsetX ) { // eslint-disable-line max-params, max-len
+	var nonunit;
+	var viewA;
+	var viewX;
+	var retmp;
+	var imtmp;
+	var magsq;
+	var remul;
+	var immul;
+	var isrm;
+	var sign;
+	var rex;
+	var imx;
+	var rea;
+	var ima;
+	var ix0;
+	var ix1;
+	var sa0;
+	var sa1;
+	var oa2;
+	var ox;
+	var sx;
+	var oa;
+	var i0;
+	var i1;
+	var ia;
+
+	// Layout
+	isrm = isRowMajor( [ strideA1, strideA2 ] );
+	nonunit = ( diag === 'non-unit' );
+
+	// Reinterpret arrays to raw numeric views
+	viewA = reinterpret( A, 0 );
+	viewX = reinterpret( x, 0 );
+
+	// Set sign to handle conjugation: flip the imaginary part for conjugate-transpose
+	if ( trans === 'conjugate-transpose' ) {
+		sign = -1;
+	} else {
+		sign = 1;
+	}
+
+	if ( isrm ) {
+		// For row-major matrices, the last dimension has the fastest changing index...
+		sa0 = strideA2 * 2;                     // offset increment for innermost loop
+		sa1 = strideA1 * 2;                     // offset increment for outermost loop
+	} else { // isColMajor
+		// For column-major matrices, the first dimension has the fastest changing index...
+		sa0 = strideA1 * 2;                     // offset increment for innermost loop
+		sa1 = strideA2 * 2;                     // offset increment for outermost loop
+	}
+
+	// Vector indexing base
+	oa = offsetA * 2;
+	ox = offsetX * 2;
+
+	// Vector strides
+	sx = strideX * 2;
+
+	if (
+		( !isrm && uplo === 'upper' && trans === 'no-transpose' ) ||
+		( isrm && uplo === 'lower' && trans !== 'no-transpose' )
+	) {
+		for ( i1 = N - 1; i1 >= 0; i1-- ) {
+			ix1 = ox + ( i1 * sx );
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1 + 1 ];
+			retmp = rex;
+			imtmp = imx;
+			oa2 = oa + ( i1 * sa1 ) + ( K * sa0 );
+			for ( i0 = i1 + 1; i0 <= min( N - 1, i1 + K ); i0++ ) {
+				ix0 = ox + ( i0 * sa0);
+				ia = oa2 + ( ( i0 - i1 ) * ( sa1 - sa0 ) );
+				rea = viewA[ ia ];
+				ima = sign * viewA[ ia + 1 ];
+				rex = viewX[ ix0 ];
+				imx = viewX[ ix0 + 1 ];
+				retmp -= f32( ( rea * rex ) - ( ima * imx ) );
+				imtmp -= f32( ( rea * imx ) + ( ima * rex ) );
+			}
+			if ( nonunit ) {
+				rea = viewA[ oa2 ];
+				ima = sign * viewA[ oa2 + 1 ];
+				magsq = f32( ( rea * rea ) + ( ima * ima ) );
+				remul = f32( ( retmp * rea ) + ( imtmp * ima ) );
+				immul = f32( ( imtmp * rea ) - ( retmp * ima ) );
+				retmp = f32( remul / magsq );
+				imtmp = f32( immul / magsq );
+			}
+			viewX[ ix1 ] = retmp;
+			viewX[ ix1 + 1 ] = imtmp;
+		}
+		return x;
+	}
+	if (
+		( !isrm && uplo === 'lower' && trans === 'no-transpose' ) ||
+		( isrm && uplo === 'upper' && trans !== 'no-transpose' )
+	) {
+		for ( i1 = 0; i1 < N; i1++ ) {
+			ix1 = ox + ( i1 * sx );
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1 + 1 ];
+			retmp = rex;
+			imtmp = imx;
+			oa2 = oa + ( i1 * sa1 );
+			for ( i0 = max( 0, i1 - K ); i0 < i1; i0++ ) {
+				ix0 = ox + ( i0 * sx );
+				ia = oa2 + ( ( i0 - i1 ) * ( sa1 - sa0 ) );
+				rea = viewA[ ia ];
+				ima = sign * viewA[ ia + 1 ];
+				rex = viewX[ ix0 ];
+				imx = viewX[ ix0 + 1 ];
+				retmp -= f32( ( rea * rex ) - ( ima * imx ) );
+				imtmp -= f32( ( rea * imx ) + ( ima * rex ) );
+			}
+			if ( nonunit ) {
+				rea = viewA[ oa2 ];
+				ima = sign * viewA[ oa2 + 1 ];
+				magsq = f32( ( rea * rea ) + ( ima * ima ) );
+				remul = f32( ( retmp * rea ) + ( imtmp * ima ) );
+				immul = f32( ( imtmp * rea ) - ( retmp * ima ) );
+				retmp = f32( remul / magsq );
+				imtmp = f32( immul / magsq );
+			}
+			viewX[ ix1 ] = retmp;
+			viewX[ ix1 + 1 ] = imtmp;
+		}
+		return x;
+	}
+	if (
+		( !isrm && uplo === 'upper' && trans !== 'no-transpose' ) ||
+		( isrm && uplo === 'lower' && trans === 'no-transpose' )
+	) {
+		for ( i1 = 0; i1 < N; i1++ ) {
+			ix1 = ox + ( i1 * sx );
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1 + 1 ];
+			retmp = rex;
+			imtmp = imx;
+			oa2 = oa + ( i1 * sa1 ) + ( K * sa0 );
+			for ( i0 = max( 0, i1 - K ); i0 < i1; i0++ ) {
+				ix0 = ox + ( i0 * sx );
+				ia = oa2 + ( ( i0 - i1 ) * sa0 );
+				rea = viewA[ ia ];
+				ima = sign * viewA[ ia + 1 ];
+				rex = viewX[ ix0 ];
+				imx = viewX[ ix0 + 1 ];
+				retmp -= f32( ( rea * rex ) - ( ima * imx ) );
+				imtmp -= f32( ( rea * imx ) + ( ima * rex ) );
+			}
+			if ( nonunit ) {
+				rea = viewA[ oa2 ];
+				ima = sign * viewA[ oa2 + 1 ];
+				magsq = f32( ( rea * rea ) + ( ima * ima ) );
+				remul = f32( ( retmp * rea ) + ( imtmp * ima ) );
+				immul = f32( ( imtmp * rea ) - ( retmp * ima ) );
+				retmp = f32( remul / magsq );
+				imtmp = f32( immul / magsq );
+			}
+			viewX[ ix1 ] = retmp;
+			viewX[ ix1 + 1 ] = imtmp;
+		}
+		return x;
+	}
+	// ( !isrm && uplo === 'lower' && trans !== 'no-transpose' ) || ( isrm && uplo === 'upper' && trans === 'no-transpose' )
+	for ( i1 = N - 1; i1 >= 0; i1-- ) {
+		ix1 = ox + ( i1 * sx );
+		rex = viewX[ ix1 ];
+		imx = viewX[ ix1 + 1 ];
+		retmp = rex;
+		imtmp = imx;
+		oa2 = oa + ( i1 * sa1 );
+		for ( i0 = min( N - 1, i1 + K ); i0 > i1; i0-- ) {
+			ix0 = ox + ( i0 * sx );
+			ia = oa2 + ( ( i0 - i1 ) * sa0 );
+			rea = viewA[ ia ];
+			ima = sign * viewA[ ia + 1 ];
+			rex = viewX[ ix0 ];
+			imx = viewX[ ix0 + 1 ];
+			retmp -= f32( ( rea * rex ) - ( ima * imx ) );
+			imtmp -= f32( ( rea * imx ) + ( ima * rex ) );
+		}
+		if ( nonunit ) {
+			rea = viewA[ oa2 ];
+			ima = sign * viewA[ oa2 + 1 ];
+			magsq = f32( ( rea * rea ) + ( ima * ima ) );
+			remul = f32( ( retmp * rea ) + ( imtmp * ima ) );
+			immul = f32( ( imtmp * rea ) - ( retmp * ima ) );
+			retmp = f32( remul / magsq );
+			imtmp = f32( immul / magsq );
+		}
+		viewX[ ix1 ] = retmp;
+		viewX[ ix1 + 1 ] = imtmp;
+	}
+	return x;
+}
+
+
+// EXPORTS //
+
+module.exports = ctbsv;

lib/node_modules/@stdlib/blas/base/ctbsv/test/fixtures/column_major_l_nt_nu.json

@@ -0,0 +1,26 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "no-transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "K": 1,
+  "A": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 0.0, 0.0 ],
+  "A_compact": [
+    [ 1.0, 1.0, 3.0, 3.0, 5.0, 5.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ]
+  ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 3.0, 3.0, 0.0, 0.0 ],
+    [ 0.0, 0.0, 4.0, 4.0, 5.0, 5.0 ]
+  ],
+  "LDA": 2,
+  "strideA1": 1,
+  "strideA2": 2,
+  "offsetA": 0,
+  "x": [ 0.0, 2.0, 0.0, 16.0, 0.0, 46.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}

lib/node_modules/@stdlib/blas/base/ctbsv/test/fixtures/column_major_l_nt_u.json

@@ -0,0 +1,26 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "no-transpose",
+  "diag": "unit",
+  "N": 3,
+  "K": 1,
+  "A": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 0.0, 0.0 ],
+  "A_compact": [
+    [ 1.0, 1.0, 3.0, 3.0, 5.0, 5.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ]
+  ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 3.0, 3.0, 0.0, 0.0 ],
+    [ 0.0, 0.0, 4.0, 4.0, 5.0, 5.0 ]
+  ],
+  "LDA": 2,
+  "strideA1": 1,
+  "strideA2": 2,
+  "offsetA": 0,
+  "x": [ 1.0, 1.0, 2.0, 6.0, 3.0, 19.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}

lib/node_modules/@stdlib/blas/base/ctbsv/test/fixtures/column_major_l_t_nu.json

@@ -0,0 +1,26 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "K": 1,
+  "A": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 0.0, 0.0 ],
+  "A_compact": [
+    [ 1.0, 1.0, 3.0, 3.0, 5.0, 5.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ]
+  ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 3.0, 3.0, 0.0, 0.0 ],
+    [ 0.0, 0.0, 4.0, 4.0, 5.0, 5.0 ]
+  ],
+  "LDA": 2,
+  "strideA1": 1,
+  "strideA2": 2,
+  "offsetA": 0,
+  "x": [ 0.0, 10.0, 0.0, 36.0, 0.0, 30.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}

lib/node_modules/@stdlib/blas/base/ctbsv/test/fixtures/column_major_l_t_u.json

@@ -0,0 +1,26 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "transpose",
+  "diag": "unit",
+  "N": 3,
+  "K": 1,
+  "A": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 0.0, 0.0 ],
+  "A_compact": [
+    [ 1.0, 1.0, 3.0, 3.0, 5.0, 5.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ]
+  ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 3.0, 3.0, 0.0, 0.0 ],
+    [ 0.0, 0.0, 4.0, 4.0, 5.0, 5.0 ]
+  ],
+  "LDA": 2,
+  "strideA1": 1,
+  "strideA2": 2,
+  "offsetA": 0,
+  "x": [ 1.0, 9.0, 2.0, 26.0, 3.0, 3.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}

lib/node_modules/@stdlib/blas/base/ctbsv/test/fixtures/column_major_u_nt_nu.json

@@ -0,0 +1,26 @@
+{
+  "order": "column-major",
+  "uplo": "upper",
+  "trans": "no-transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "K": 1,
+  "A": [ 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0 ],
+  "A_compact": [
+    [ 0.0, 0.0, 2.0, 2.0, 4.0, 4.0 ],
+    [ 1.0, 1.0, 3.0, 3.0, 5.0, 5.0 ]
+  ],
+  "A_mat": [
+    [ 1.0, 1.0, 2.0, 2.0, 0.0, 0.0 ],
+    [ 0.0, 0.0, 3.0, 3.0, 4.0, 4.0 ],
+    [ 0.0, 0.0, 0.0, 0.0, 5.0, 5.0 ]
+  ],
+  "LDA": 2,
+  "strideA1": 1,
+  "strideA2": 2,
+  "offsetA": 0,
+  "x": [ 0.0, 10.0, 0.0, 36.0, 0.0, 30.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}