feat: add blas/base/cgbmv

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

@@ -0,0 +1,292 @@
+/**
+* @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.
+*/
+
+/* eslint-disable max-params */
+
+'use strict';
+
+// MODULES //
+
+var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
+var cfill = require( '@stdlib/blas/ext/base/cfill' ).ndarray;
+var cscal = require( '@stdlib/blas/base/cscal' ).ndarray;
+var realf = require( '@stdlib/complex/float32/real' );
+var imagf = require( '@stdlib/complex/float32/imag' );
+var f32 = require( '@stdlib/number/float64/base/to-float32' );
+var reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
+var muladd = require( '@stdlib/complex/float32/base/mul-add' ).assign;
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+
+
+// FUNCTIONS //
+
+/**
+* Tests whether a provided string indicates to transpose a matrix.
+*
+* @private
+* @param {string} str - input string
+* @returns {boolean} boolean indicating whether to transpose a matrix
+*
+* @example
+* var bool = isTransposed( 'transpose' );
+* // returns true
+*
+* @example
+* var bool = isTransposed( 'conjugate-transpose' );
+* // returns true
+*
+* @example
+* var bool = isTransposed( 'no-transpose' );
+* // returns false
+*/
+function isTransposed( str ) { // NOTE: consider moving to a separate helper utility package
+	return ( str !== 'no-transpose' );
+}
+
+
+// MAIN //
+
+/**
+* Performs one of the matrix-vector operations `y = α*A*x + β*y` or `y = α*A^T*x + β*y` or `y = α*A^H*x + β*y`, where `α` and `β` are complex scalars, `x` and `y` are complex vectors, and `A` is an `M` by `N` band matrix with `KL` sub-diagonals and `KU` super-diagonals.
+*
+* @private
+* @param {string} trans - specifies whether `A` should be transposed, conjugate-transposed, or not transposed
+* @param {NonNegativeInteger} M - number of rows in the matrix `A`
+* @param {NonNegativeInteger} N - number of columns in the matrix `A`
+* @param {NonNegativeInteger} KL - number of sub-diagonals of matrix `A`
+* @param {NonNegativeInteger} KU - number of super-diagonals of matrix `A`
+* @param {Complex64} alpha - complex scalar constant
+* @param {Complex64Array} A - input banded complex matrix
+* @param {integer} strideA1 - stride of first dimension of `A`
+* @param {integer} strideA2 - stride of second dimension of `A`
+* @param {NonNegativeInteger} offsetA - starting index into `A`
+* @param {Complex64Array} x - first input complex vector `x`
+* @param {integer} strideX - `x` stride length
+* @param {NonNegativeInteger} offsetX - starting index for `x`
+* @param {Complex64} beta - complex scalar constant
+* @param {Complex64Array} y - second input complex vector y
+* @param {integer} strideY - `y` stride length
+* @param {NonNegativeInteger} offsetY - starting index for `y`
+* @returns {Complex64Array} `y`
+*
+* @example
+* var Complex64Array = require( '@stdlib/array/complex64' );
+* var Complex64 = require( '@stdlib/complex/float32/ctor' );
+*
+* var A = new Complex64Array( [ 0.0, 0.0, 3.0, 3.0, 6.0, 6.0, 1.0, 1.0, 4.0, 4.0, 7.0, 7.0, 2.0, 2.0, 5.0, 5.0, 0.0, 0.0 ] );
+* var x = new Complex64Array( [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ] );
+* var y = new Complex64Array( [ 3.0, 3.0, 2.0, 2.0, 1.0, 1.0 ] );
+* var alpha = new Complex64( 0.5, 0.5 );
+* var beta = new Complex64( 0.5, -0.5 );
+*
+* cgbmv( 'no-transpose', 3, 3, 1, 1, alpha, A, 3, 1, 0, x, 1, 0, beta, y, 1, 0 );
+* // y => <Complex64Array>[ -4.0, 7.0, -26.0, 28.0, -30.0, 31.0 ]
+*/
+function cgbmv( trans, M, N, KL, KU, alpha, A, strideA1, strideA2, offsetA, x, strideX, offsetX, beta, y, strideY, offsetY ) {
+	var i0start;
+	var realpha;
+	var imalpha;
+	var rebeta;
+	var imbeta;
+	var viewA;
+	var viewX;
+	var viewY;
+	var i0end;
+	var retmp;
+	var imtmp;
+	var isrm;
+	var xlen;
+	var ylen;
+	var sign;
+	var kup1;
+	var sa0;
+	var sa1;
+	var rea;
+	var ima;
+	var rex;
+	var imx;
+	var oa2;
+	var oa;
+	var ox;
+	var oy;
+	var sx;
+	var sy;
+	var i0;
+	var i1;
+	var ia;
+	var ix;
+	var iy;
+
+	// layout
+	isrm = isRowMajor( [ strideA1, strideA2 ] );
+
+	// lengths
+	if ( isTransposed( trans ) ) {
+		xlen = M;
+		ylen = N;
+	} else {
+		xlen = N;
+		ylen = M;
+	}
+
+	// decompose complex scalars
+	rebeta = realf( beta );
+	imbeta = imagf( beta );
+	realpha = realf( alpha );
+	imalpha = imagf( alpha );
+
+	// y = beta*y
+	if ( rebeta === 0.0 && imbeta === 0.0 ) {
+		cfill( ylen, 0.0, y, strideY, offsetY );
+	} else if ( rebeta !== 1.0 || imbeta !== 0.0 ) {
+		cscal( ylen, beta, y, strideY, offsetY );
+	}
+
+	// if alpha is zero, early return y
+	if ( realpha === 0.0 && imalpha === 0.0 ) {
+		return y;
+	}
+
+	// reinterpret arrays to raw numeric views
+	viewA = reinterpret( A, 0 );
+	viewX = reinterpret( x, 0 );
+	viewY = reinterpret( y, 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
+	}
+	console.log('sa0', sa0);
+	console.log('sa1', sa1);
+
+	// Vector indexing base
+	oa = offsetA * 2;
+	ox = offsetX * 2;
+	oy = offsetY * 2;
+
+	// Vector strides
+	sx = strideX * 2;
+	sy = strideY * 2;
+
+	// Form: y = α*A*x + y
+	if (
+		( !isrm && !isTransposed( trans ) ) ||
+		( isrm && isTransposed( trans ) )
+	) {
+		console.log('1nd branch');
+		for ( i1 = 0; i1 < xlen; i1++ ) {
+			ix = ox + ( i1 * sx );
+			rex = viewX[ ix ];
+			imx = viewX[ ix + 1 ];
+			retmp = f32( ( realpha * rex ) - ( imalpha * imx ) );
+			imtmp = f32( ( realpha * imx ) + ( imalpha * rex ) );
+			oa2 = oa + ( i1 * sa1 );
+			i0start = max( 0, i1 - KU );
+			i0end = min( ylen, i1 + KL + 1 );
+			for ( i0 = i0start; i0 < i0end; i0++ ) {
+				kup1 = i0 - i1;
+				ia = oa2 + ( ( kup1 + KU ) * sa0 );
+				if ( ia >= 0 ) {
+					rea = viewA[ ia ];
+					ima = sign * viewA[ ia + 1 ];
+					iy = oy + ( i0 * sy );
+					muladd( rea, ima, retmp, imtmp, viewY[ iy ], viewY[ iy + 1 ], viewY, 1, iy ); // eslint-disable-line max-len
+				}
+			}
+		}
+		return y;
+	}
+
+	// Form: y = α*A^T*x + y
+	// ( !isrm && isTransposed( trans ) ) || ( isrm && !isTransposed( trans ) )
+	console.log('2nd branch');
+	for ( i1 = 0; i1 < ylen; i1++ ) {
+		retmp = 0.0;
+		imtmp = 0.0;
+		oa2 = oa + ( i1 * sa1 );
+		i0start = max( 0, i1 - KU );
+		i0end = min( xlen, i1 + KL + 1 );
+		for ( i0 = i0start; i0 < i0end; i0++ ) {
+			kup1 = i0 - i1;
+			ia = oa2 + ( ( kup1 + KU ) * sa0 );
+			if ( i0 < xlen && ia >= 0 ) {
+				ix = ox + ( i0 * sx );
+				rea = viewA[ ia ];
+				ima = sign * viewA[ ia + 1 ];
+				rex = viewX[ ix ];
+				imx = viewX[ ix + 1 ];
+				retmp = f32( retmp + f32( ( rea * rex ) - ( ima * imx ) ) );
+				imtmp = f32( imtmp + f32( ( rea * imx ) + ( ima * rex ) ) );
+			}
+		}
+		iy = oy + ( i1 * sy );
+		muladd( realpha, imalpha, retmp, imtmp, viewY[ iy ], viewY[ iy + 1 ], viewY, 1, iy ); // eslint-disable-line max-len
+	}
+	return y;
+}
+
+
+// EXPORTS //
+
+module.exports = cgbmv;
+
+var Complex64Array = require( '@stdlib/array/complex64' );
+var Complex64 = require( '@stdlib/complex/float32/ctor' );
+
+// A_Canonical = [
+// 	[ 1.0, 1.0, 3.0, 3.0, 0.0, 0.0 ],
+// 	[ 2.0, 2.0, 4.0, 4.0, 6.0, 6.0 ],
+// 	[ 0.0, 0.0, 5.0, 5.0, 7.0, 7.0 ]
+// ]
+
+// A_Compact = [
+// 	[ 0.0, 0.0, 3.0, 3.0, 6.0, 6.0 ],
+//  [ 1.0, 1.0, 4.0, 4.0, 7.0, 7.0 ],
+//  [ 2.0, 2.0, 5.0, 5.0, 0.0, 0.0 ],
+// ]
+
+var alpha = new Complex64( 0.5, 0.5 );
+var beta = new Complex64( 0.5, -0.5 );
+var x = new Complex64Array( [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ] );
+var y1 = new Complex64Array( [ 3.0, 3.0, 2.0, 2.0, 1.0, 1.0 ] );
+var y2 = new Complex64Array( [ 3.0, 3.0, 2.0, 2.0, 1.0, 1.0 ] );
+
+// row-major
+var A_rm = new Complex64Array( [ 0.0, 0.0, 3.0, 3.0, 6.0, 6.0, 1.0, 1.0, 4.0, 4.0, 7.0, 7.0, 2.0, 2.0, 5.0, 5.0, 0.0, 0.0 ] );
+cgbmv( 'no-transpose', 3, 3, 1, 1, alpha, A_rm, 3, 1, 0, x, 1, 0, beta, y1, 1, 0 );
+console.log( 'row-major', y1.toString() );
+// y => <Complex64Array>[ -12.0, 15.0, -28.0, 30.0,-18.0, 19.0 ]
+
+// col-major
+var A_cm = 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, 6.0, 6.0, 7.0, 7.0, 0.0, 0.0 ] );
+cgbmv( 'no-transpose', 3, 3, 1, 1, alpha, A_cm, 1, 3, 0, x, 1, 0, beta, y2, 1, 0 );
+// y => <Complex64Array>[ -4.0, 7.0, -26.0, 28.0, -30.0, 31.0 ]
+console.log( 'col-major', y2.toString() );