feat: add lapack/base/dgebal

lib/node_modules/@stdlib/lapack/base/dgebal/lib/base.js

@@ -0,0 +1,283 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 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-len, max-params, max-statements */
+
+'use strict';
+
+// MODULES //
+
+var dswap = require( '@stdlib/blas/base/dswap' ).ndarray;
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray;
+var idamax = require( '@stdlib/blas/base/idamax' ).ndarray;
+var abs = require( '@stdlib/math/base/special/abs' );
+var isnan = require( '@stdlib/assert/is-nan' );
+var max = require( '@stdlib/math/base/special/maxn' );
+var min = require( '@stdlib/math/base/special/minn' );
+var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
+
+
+// MAIN //
+
+/**
+* Balances a general real matrix `A`.
+*
+* ## Notes
+*
+* The job parameter can be one of the following:
+*
+* -   'N': none, return immediately
+* -   'P': permute only
+* -   'S': scale only
+* -   'B': both permute and scale
+* -   The matrix `A` is overwritten by the balanced matrix.
+*
+* @private
+* @param {string} job - indicates the operations to be performed
+* @param {NonNegativeInteger} N - number of rows/columns in matrix `A`
+* @param {Float64Array} A - input matrix to be balanced
+* @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 {Float64Array} out - stores the first and last row/column of the balanced submatrix
+* @param {integer} strideOut - stride of `out`
+* @param {NonNegativeInteger} offsetOut - starting index for `out`
+* @param {Float64Array} scale - array containing permutation and scaling information
+* @param {integer} strideScale - stride of `scale`
+* @param {NonNegativeInteger} offsetScale  - starting index for `scale`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 100.0, 0.0, 2.0, 200.0, 0.0, 0.0, 0.0, 3.0 ] );
+* var out = new Float64Array( 2 );
+* var scale = new Float64Array( 3 );
+*
+* dgebal( 'B', 3, A, 3, 1, 0, out, 1, 0, scale, 1, 0 );
+* // A => <Float64Array>[ 1, 12.5, 0, 16, 200, 0, 0, 0, 3 ]
+* // out => <Float64Array>[ 0, 1 ]
+* // scale => <Float64Array>[ 8, 1, 2 ]
+*/
+function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetOut, scale, strideScale, offsetScale ) {
+	var canSwap;
+	var noconv;
+	var sfmin1;
+	var sfmin2;
+	var sfmax1;
+	var sfmax2;
+	var sclfac;
+	var factor;
+	var ica;
+	var ira;
+	var ca;
+	var ra;
+	var is;
+	var c;
+	var r;
+	var k;
+	var l;
+	var i;
+	var j;
+	var g;
+	var f;
+	var s;
+
+	sclfac = 2.0;
+	factor = 0.95;
+
+	// Quick return if possible
+	if ( N === 0 ) {
+		out[ offsetOut ] = 0.0; // ilo
+		out[ offsetOut + strideOut ] = -1.0; // ihi (invalid)
+		return 0;
+	}
+
+	if ( job === 'N' ) {
+		is = offsetScale;
+		for ( i = 0; i < N; i++ ) {
+			scale[ is ] = 1.0;
+			is += strideScale;
+		}
+
+		out[ offsetOut ] = 0.0; // ilo
+		out[ offsetOut + strideOut ] = N - 1; // ihi
+		return 0;
+	}
+
+	// Permutation to isolate eigenvalues if possible
+	k = 0;
+	l = N - 1;
+
+	if ( job !== 'S' ) {
+		// Row and column exchange
+		noconv = true;
+		while ( noconv ) {
+			// Search for rows isolating an eigenvalue and push them down
+			noconv = false;
+			for ( i = l; i >= 0; i-- ) {
+				canSwap = true;
+				for ( j = 0; j <= l; j++ ) {
+					if ( i !== j && A[ offsetA + (i*strideA1) + (j*strideA2) ] !== 0.0 ) {
+						canSwap = false;
+						break;
+					}
+				}
+
+				if ( canSwap ) {
+					scale[ offsetScale + (l*strideScale) ] = i;
+					if ( i !== l ) {
+						dswap( l+1, A, strideA1, offsetA + (i*strideA2), A, strideA1, offsetA + (l*strideA2) );
+						dswap( N - k, A, strideA2, offsetA + (i*strideA1) + (k*strideA2), A, strideA2, offsetA + (l*strideA1) + (k*strideA2) );
+					}
+					noconv = true;
+
+					if ( l === 0.0 ) {
+						out[ offsetOut ] = 0.0; // ilo
+						out[ offsetOut + strideOut ] = 0.0; // ihi
+						return 0;
+					}
+					l -= 1;
+				}
+			}
+		}
+
+		noconv = true;
+		while ( noconv ) {
+			// Search for columns isolating an eigenvalue and push them left
+			noconv = false;
+			for ( j = k; j <= l; j++ ) {
+				canSwap = true;
+				for ( i = k; i <= l; i++ ) {
+					if ( i !== j && A[ offsetA + (i*strideA1) + (j*strideA2) ] !== 0.0 ) {
+						canSwap = false;
+						break;
+					}
+				}
+
+				if ( canSwap ) {
+					scale[ offsetScale + (k*strideScale) ] = j;
+					if ( j !== k ) {
+						dswap( l+1, A, strideA1, offsetA + (j*strideA2), A, strideA1, offsetA + (k*strideA2) );
+						dswap( N-k, A, strideA2, offsetA + (j*strideA1), A, strideA2, offsetA + (k*strideA1) );
+					}
+					noconv = true;
+					k += 1;
+				}
+			}
+		}
+	}
+
+	// Initialize `scale` for non-permuted submatrix
+	is = offsetScale;
+	for ( i = k; i <= l; i++ ) {
+		scale[ is ] = 1.0;
+		is += strideScale;
+	}
+
+	if ( job === 'P' ) {
+		out[ offsetOut ] = k; // ilo
+		out[ offsetOut + strideOut ] = l; // ihi
+		return 0;
+	}
+
+	// Balance the submatrix in rows K to L, iterative loop for norm reduction
+	sfmin1 = dlamch( 'S' ) / dlamch( 'P' );
+	sfmax1 = 1.0 / sfmin1;
+	sfmin2 = sfmin1 * sclfac;
+	sfmax2 = 1.0 / sfmin2;
+
+	noconv = true;
+	while ( noconv ) {
+		noconv = false;
+		for ( i = k; i <= l; i++ ) {
+			c = dnrm2( l-k+1, A, strideA1, offsetA + (k*strideA1) + (i*strideA2) );
+			r = dnrm2( l-k+1, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) );
+			ica = idamax( l+1, A, strideA1, offsetA + (i*strideA2) );
+			ca = abs( A[ offsetA + (ica*strideA1) + (i*strideA2) ] );
+			ira = idamax( N-k+1, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) );
+			ra = abs( A[ offsetA + (i*strideA1) + ((ira+k)*strideA2) ] );
+
+			if ( c === 0.0 || r === 0.0 ) {
+				continue;
+			}
+
+			if ( isnan( c ) || isnan( r ) || isnan( ca ) || isnan( ra ) ) {
+				return -3;
+			}
+
+			g = r / sclfac;
+			f = 1.0;
+			s = c + r;
+
+			while ( c < g && max( f, c, ca ) < sfmax2 && min( r, g, ra ) > sfmin2 ) {
+				f *= sclfac;
+				c *= sclfac;
+				ca *= sclfac;
+				r /= sclfac;
+				g /= sclfac;
+				ra /= sclfac;
+			}
+
+			g = c / sclfac;
+
+			while ( g >= r && max( r, ra ) < sfmax2 && min( f, c, g, ca ) > sfmin2 ) {
+				f /= sclfac;
+				c /= sclfac;
+				g /= sclfac;
+				ca /= sclfac;
+				r *= sclfac;
+				ra *= sclfac;
+			}
+
+			// Now balance
+			if ( ( c + r ) >= factor * s ) {
+				continue;
+			}
+
+			if ( f < 1.0 && scale[ offsetScale + (i*strideScale) ] < 1.0 ) {
+				if ( f * scale[ offsetScale + (i*strideScale) ] <= sfmin1 ) {
+					continue;
+				}
+			}
+
+			if ( f > 1.0 && scale[ offsetScale + (i*strideScale) ] > 1.0 ) {
+				if ( scale[ offsetScale + (i*strideScale) ] >= sfmax1 / f ) {
+					continue;
+				}
+			}
+
+			g = 1.0 / f;
+			scale[ offsetScale + (i*strideScale) ] *= f;
+			noconv = true;
+
+			dscal( N-k, g, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) );
+			dscal( l+1, f, A, strideA1, offsetA + (i*strideA2) );
+		}
+	}
+
+	out[ offsetOut ] = k; // ilo
+	out[ offsetOut + strideOut ] = l; // ihi
+	return 0;
+}
+
+
+// EXPORTS //
+
+module.exports = dgebal;

lib/node_modules/@stdlib/lapack/base/dgebal/lib/dgebal.js

@@ -0,0 +1,96 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 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 isLayout = require( '@stdlib/blas/base/assert/is-layout' );
+var isColumnMajor = require( '@stdlib/ndarray/base/assert/is-column-major-string' );
+var format = require( '@stdlib/string/format' );
+var base = require( './base.js' );
+
+
+// MAIN //
+
+/**
+* Balances a general real matrix `A`.
+*
+* ## Notes
+*
+* The job parameter can be one of the following:
+*
+* -   'N': none, return immediately
+* -   'P': permute only
+* -   'S': scale only
+* -   'B': both permute and scale
+* -   The matrix `A` is overwritten by the balanced matrix.
+*
+* @private
+* @param {string} order - storage layout of `A`
+* @param {string} job - indicates the operations to be performed
+* @param {NonNegativeInteger} N - number of rows/columns in matrix `A`
+* @param {Float64Array} A - input matrix to be balanced
+* @param {NonNegativeInteger} LDA - leading dimension of `A`
+* @param {Float64Array} out - stores the first and last row/column of the balanced submatrix
+* @param {Float64Array} scale - array containing permutation and scaling information
+* @throws {TypeError} first argument must be a valid order
+* @throws {TypeError} second argument must be a valid job
+* @throws {RangeError} fifth argument must be greater than or equal to `N`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 100.0, 0.0, 2.0, 200.0, 0.0, 0.0, 0.0, 3.0 ] );
+* var out = new Float64Array( 2 );
+* var scale = new Float64Array( 3 );
+*
+* dgebal( 'row-major', 'B', 3, A, 3, out, scale );
+* // A => <Float64Array>[ 1, 12.5, 0, 16, 200, 0, 0, 0, 3 ]
+* // out => <Float64Array>[ 0, 1 ]
+* // scale => <Float64Array>[ 8, 1, 2 ]
+*/
+function dgebal( order, job, N, A, LDA, out, scale ) {
+	var sa1;
+	var sa2;
+
+	if ( !isLayout( order ) ) {
+		throw new TypeError( format( 'invalid argument. First argument must be a valid order. Value: `%s`.', order ) );
+	}
+	if ( job !== 'B' && job !== 'S' && job !== 'P' && job !== 'N' ) {
+		throw new TypeError( format( 'invalid argument. Second argument must be one of the following: `B`, `S`, `P`, or `N`. Value: `%s`.', job ) );
+	}
+	if ( isColumnMajor( order ) ) {
+		sa1 = 1;
+		sa2 = LDA;
+	} else { // order === 'row-major'
+		if ( LDA < N ) {
+			throw new RangeError( format( 'invalid argument. Eighth argument must be greater than or equal to %d. Value: `%d`.', N, LDA ) );
+		}
+		sa1 = LDA;
+		sa2 = 1;
+	}
+
+	return base( job, N, A, sa1, sa2, 0, out, 1, 0, scale, 1, 0 );
+}
+
+
+// EXPORTS //
+
+module.exports = dgebal;