feat: add blas/base/ctpsv

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

@@ -0,0 +1,262 @@
+/*
+* @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 isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major-string' );
+var f32 = require( '@stdlib/number/float64/base/to-float32' );
+var reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
+
+
+// MAIN //
+
+/**
+* Solves one of the systems of equations `A*x = b` or `A^T*x = b` or `A^H*x = b` where `b` and `x` are `N` element vectors and `A` is an `N` by `N` unit, or non-unit, upper or lower triangular matrix, supplied in packed form.
+*
+* @private
+* @param {string} order - storage layout
+* @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 {Complex64Array} AP - packed form of a symmetric matrix `A`
+* @param {integer} strideAP - `AP` stride length
+* @param {NonNegativeInteger} offsetAP - starting index for `AP`
+* @param {Complex64Array} x - input vector
+* @param {integer} strideX - `x` stride length
+* @param {NonNegativeInteger} offsetX - starting index for `x`
+* @returns {Complex64Array} `x`
+*
+* @example
+* var Complex64Array = require( '@stdlib/array/complex64' );
+*
+* var AP = new Complex64Array( [ 1.0, 1.0, 2.0, 2.0, 4.0, 4.0, 3.0, 3.0, 5.0, 5.0, 6.0, 6.0 ] );
+* var x = new Complex64Array( [ 0.0, 2.0, 0.0, 20.0, 0.0, 62.0 ] );
+*
+* ctpsv( 'row-major', 'lower', 'no-transpose', 'non-unit', 3, AP, 1, 0, x, 1, 0 );
+* // x => <Complex64Array>[ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+*/
+function ctpsv( order, uplo, trans, diag, N, AP, strideAP, offsetAP, x, strideX, offsetX ) { // eslint-disable-line max-params, max-len
+	var nonunit;
+	var viewAP;
+	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 iap;
+	var ix0;
+	var ix1;
+	var sap;
+	var kk;
+	var ox;
+	var sx;
+	var i1;
+	var i0;
+
+	// Note on variable naming convention: sa#, ix#, i# where # corresponds to the loop number, with `0` being the innermost loop...
+
+	isrm = isRowMajor( order );
+	nonunit = ( diag === 'non-unit' );
+
+	// Reinterpret arrays as real-valued views of interleaved real and imaginary components:
+	viewAP = reinterpret( AP, 0 );
+	viewX = reinterpret( x, 0 );
+	if ( trans === 'conjugate-transpose' ) {
+		sign = -1;
+	} else {
+		sign = 1;
+	}
+	ox = offsetX * 2;
+	kk = offsetAP * 2;
+	sx = strideX * 2;
+	sap = strideAP * 2;
+
+	if (
+		( !isrm && uplo === 'upper' && trans === 'no-transpose' ) ||
+		( isrm && uplo === 'lower' && trans !== 'no-transpose' )
+	) {
+		ix1 = ox + ((N-1)*sx);
+		kk += ((N*(N+1)/2)-1) * sap;
+		for ( i1 = N - 1; i1 >= 0; i1-- ) {
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1+1 ];
+			if ( rex !== 0.0 || imx !== 0.0 ) {
+				iap = kk;
+				if ( nonunit ) {
+					rea = viewAP[ iap ];
+					ima = sign * viewAP[ iap+1 ];
+					magsq = f32(f32(rea*rea) + f32(ima*ima));
+					retmp = f32(f32(f32(rex*rea) + f32(imx*ima)) / magsq );
+					imtmp = f32(f32(f32(imx*rea) - f32(rex*ima)) / magsq );
+					viewX[ ix1 ] = retmp;
+					viewX[ ix1+1 ] = imtmp;
+				} else {
+					retmp = rex;
+					imtmp = imx;
+				}
+				ix0 = ix1 - sx;
+				iap -= sap;
+				for ( i0 = i1 - 1; i0 >= 0; i0-- ) {
+					rea = viewAP[ iap ];
+					ima = sign * viewAP[ iap+1 ];
+					rex = viewX[ ix0 ];
+					imx = viewX[ ix0+1 ];
+					remul = f32(f32(retmp*rea) - f32(imtmp*ima));
+					immul = f32(f32(retmp*ima) + f32(imtmp*rea));
+					viewX[ ix0 ] = f32(rex - remul);
+					viewX[ ix0+1 ] = f32(imx - immul);
+					ix0 -= sx;
+					iap -= sap;
+				}
+			}
+			ix1 -= sx;
+			kk -= (i1+1) * sap;
+		}
+		return x;
+	}
+	if (
+		( !isrm && uplo === 'lower' && trans === 'no-transpose' ) ||
+		( isrm && uplo === 'upper' && trans !== 'no-transpose' )
+	) {
+		ix1 = ox;
+		for ( i1 = 0; i1 < N; i1++ ) {
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1+1 ];
+			if ( rex !== 0.0 || imx !== 0.0 ) {
+				iap = kk;
+				if ( nonunit ) {
+					rea = viewAP[ iap ];
+					ima = sign * viewAP[ iap+1 ];
+					magsq = f32(f32(rea*rea) + f32(ima*ima));
+					retmp = f32(f32(f32(rex*rea) + f32(imx*ima)) / magsq );
+					imtmp = f32(f32(f32(imx*rea) - f32(rex*ima)) / magsq );
+					viewX[ ix1 ] = retmp;
+					viewX[ ix1+1 ] = imtmp;
+				} else {
+					retmp = rex;
+					imtmp = imx;
+				}
+				ix0 = ix1 + sx;
+				iap += sap;
+				for ( i0 = i1 + 1; i0 < N; i0++ ) {
+					rea = viewAP[ iap ];
+					ima = sign * viewAP[ iap+1 ];
+					rex = viewX[ ix0 ];
+					imx = viewX[ ix0+1 ];
+					remul = f32(f32(retmp*rea) - f32(imtmp*ima));
+					immul = f32(f32(retmp*ima) + f32(imtmp*rea));
+					viewX[ ix0 ] = f32(rex - remul);
+					viewX[ ix0+1 ] = f32(imx - immul);
+					ix0 += sx;
+					iap += sap;
+				}
+			}
+			ix1 += sx;
+			kk += (N-i1) * sap;
+		}
+		return x;
+	}
+	if (
+		( !isrm && uplo === 'upper' && trans !== 'no-transpose' ) ||
+		( isrm && uplo === 'lower' && trans === 'no-transpose' )
+	) {
+		ix1 = ox;
+		for ( i1 = 0; i1 < N; i1++ ) {
+			rex = viewX[ ix1 ];
+			imx = viewX[ ix1+1 ];
+			retmp = rex;
+			imtmp = imx;
+			ix0 = ox;
+			iap = kk;
+			for ( i0 = 0; i0 < i1; i0++ ) {
+				rea = viewAP[ iap ];
+				ima = sign * viewAP[ iap+1 ];
+				rex = viewX[ ix0 ];
+				imx = viewX[ ix0+1 ];
+				retmp = f32(retmp - f32(f32(rex*rea) - f32(imx*ima)));
+				imtmp = f32(imtmp - f32(f32(rex*ima) + f32(imx*rea)));
+				ix0 += sx;
+				iap += sap;
+			}
+			if ( nonunit ) {
+				rea = viewAP[ iap ];
+				ima = sign * viewAP[ iap+1 ];
+				magsq = f32(f32(rea*rea) + f32(ima*ima));
+				remul = f32(f32(retmp*rea) + f32(imtmp*ima));
+				immul = f32(f32(imtmp*rea) - f32(retmp*ima));
+				retmp = f32(remul / magsq);
+				imtmp = f32(immul / magsq);
+			}
+			viewX[ ix1 ] = retmp;
+			viewX[ ix1+1 ] = imtmp;
+			ix1 += sx;
+			kk += (i1+1) * sap;
+		}
+		return x;
+	}
+	// ( !isrm && uplo === 'lower' && trans !== 'no-transpose' ) || ( isrm && uplo === 'upper' && trans === 'no-transpose' )
+	ix1 = ox + ((N-1)*sx);
+	kk += ((N*(N+1)/2)-1) * sap;
+	for ( i1 = N - 1; i1 >= 0; i1-- ) {
+		rex = viewX[ ix1 ];
+		imx = viewX[ ix1+1 ];
+		retmp = rex;
+		imtmp = imx;
+		ix0 = ox + ((N-1)*sx);
+		iap = kk;
+		for ( i0 = N - 1; i0 > i1; i0-- ) {
+			rea = viewAP[ iap ];
+			ima = sign * viewAP[ iap+1 ];
+			rex = viewX[ ix0 ];
+			imx = viewX[ ix0+1 ];
+			retmp = f32(retmp - f32(f32(rex*rea) - f32(imx*ima)));
+			imtmp = f32(imtmp - f32(f32(rex*ima) + f32(imx*rea)));
+			ix0 -= sx;
+			iap -= sap;
+		}
+		if ( nonunit ) {
+			rea = viewAP[ iap ];
+			ima = sign * viewAP[ iap+1 ];
+			magsq = f32(f32(rea*rea) + f32(ima*ima));
+			remul = f32(f32(retmp*rea) + f32(imtmp*ima));
+			immul = f32(f32(imtmp*rea) - f32(retmp*ima));
+			retmp = f32(remul / magsq);
+			imtmp = f32(immul / magsq);
+		}
+		viewX[ ix1 ] = retmp;
+		viewX[ ix1+1 ] = imtmp;
+		ix1 -= sx;
+		kk -= (N-i1) * sap;
+	}
+	return x;
+}
+
+
+// EXPORTS //
+
+module.exports = ctpsv;

lib/node_modules/@stdlib/blas/base/ctpsv/package.json

@@ -0,0 +1,71 @@
+{
+  "name": "@stdlib/blas/base/ctpsv",
+  "version": "0.0.0",
+  "description": "Solves one of the systems of equations `A*x = b` or `A^T*x = b` or `A^H*x = b`.",
+  "license": "Apache-2.0",
+  "author": {
+    "name": "The Stdlib Authors",
+    "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
+  },
+  "contributors": [
+    {
+      "name": "The Stdlib Authors",
+      "url": "https://github.com/stdlib-js/stdlib/graphs/contributors"
+    }
+  ],
+  "main": "./lib",
+  "browser": "./lib/main.js",
+  "directories": {
+    "benchmark": "./benchmark",
+    "doc": "./docs",
+    "example": "./examples",
+    "lib": "./lib",
+    "test": "./test"
+  },
+  "types": "./docs/types",
+  "scripts": {},
+  "homepage": "https://github.com/stdlib-js/stdlib",
+  "repository": {
+    "type": "git",
+    "url": "git://github.com/stdlib-js/stdlib.git"
+  },
+  "bugs": {
+    "url": "https://github.com/stdlib-js/stdlib/issues"
+  },
+  "dependencies": {},
+  "devDependencies": {},
+  "engines": {
+    "node": ">=0.10.0",
+    "npm": ">2.7.0"
+  },
+  "os": [
+    "aix",
+    "darwin",
+    "freebsd",
+    "linux",
+    "macos",
+    "openbsd",
+    "sunos",
+    "win32",
+    "windows"
+  ],
+  "keywords": [
+    "stdlib",
+    "stdmath",
+    "mathematics",
+    "math",
+    "blas",
+    "level 2",
+    "ctpsv",
+    "hermitian",
+    "linear",
+    "algebra",
+    "subroutines",
+    "array",
+    "ndarray",
+    "complex",
+    "complex64",
+    "complex64array",
+    "single"
+  ]
+}

lib/node_modules/@stdlib/blas/base/ctpsv/test/fixtures/column_major_complex_access_pattern.json

@@ -0,0 +1,19 @@
+{
+  "order": "column-major",
+  "uplo": "upper",
+  "trans": "no-transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "strideAP": -3,
+  "offsetAP": 15,
+  "AP": [ 6.0, 6.0, 999.0, 999.0, 999.0, 999.0, 5.0, 5.0, 999.0, 999.0, 999.0, 999.0, 3.0, 3.0, 999.0, 999.0, 999.0, 999.0, 4.0, 4.0, 999.0, 999.0, 999.0, 999.0, 2.0, 2.0, 999.0, 999.0, 999.0, 999.0, 1.0, 1.0 ],
+  "A_mat": [
+    [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ],
+    [ 0.0, 0.0, 4.0, 4.0, 5.0, 5.0 ],
+    [ 0.0, 0.0, 0.0, 0.0, 6.0, 6.0 ]
+  ],
+  "x": [ 0.0, 36.0, 0.0, 46.0, 0.0, 28.0 ],
+  "strideX": -1,
+  "offsetX": 2,
+  "x_out": [ 3.0, 3.0, 2.0, 2.0, 1.0, 1.0 ]
+}

lib/node_modules/@stdlib/blas/base/ctpsv/test/fixtures/column_major_l_ct_nu.json

@@ -0,0 +1,19 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "conjugate-transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "strideAP": 1,
+  "offsetAP": 0,
+  "AP": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 6.0, 6.0 ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ],
+    [ 3.0, 3.0, 5.0, 5.0, 6.0, 6.0 ]
+  ],
+  "x": [ 28.0, 0.0, 46.0, 0.0, 36.0, 0.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/ctpsv/test/fixtures/column_major_l_nt_nu.json

@@ -0,0 +1,19 @@
+{
+  "order": "column-major",
+  "uplo": "lower",
+  "trans": "no-transpose",
+  "diag": "non-unit",
+  "N": 3,
+  "strideAP": 1,
+  "offsetAP": 0,
+  "AP": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 6.0, 6.0 ],
+  "A_mat": [
+    [ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 ],
+    [ 2.0, 2.0, 4.0, 4.0, 0.0, 0.0 ],
+    [ 3.0, 3.0, 5.0, 5.0, 6.0, 6.0 ]
+  ],
+  "x": [ 0.0, 2.0, 0.0, 20.0, 0.0, 62.0 ],
+  "strideX": 1,
+  "offsetX": 0,
+  "x_out": [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0 ]
+}