feat: add lapack/base/dorgbr

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

@@ -0,0 +1,171 @@
+/**
+* @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-len, max-params */
+
+'use strict';
+
+// MODULES //
+
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var dorglq = require( './dorglq.js' );
+var dorgqr = require( './dorgqr.js' );
+
+
+// MAIN //
+
+/**
+* Generates one of the real orthogonal matrices `Q` or `P**T` determined by `DGEBRD` when reducing a real matrix `A` to bi-diagonal form: `A = Q * B * P^T`.
+* `Q` and `P^T` are defined as products of elementary reflectors `H(i)` or `G(i)` respectively.
+*
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows of `Q`
+* @param {NonNegativeInteger} N - number of columns of `Q`
+* @param {NonNegativeInteger} K - number of elementary reflectors
+* @param {Float64Array} A - input/output matrix
+* @param {integer} strideA1 - stride length of the first dimension of `A`
+* @param {integer} strideA2 - stride length of the second dimension of `A`
+* @param {NonNegativeInteger} offsetA - starting index for `A`
+* @param {Float64Array} TAU - scalar factors of reflectors
+* @param {integer} strideTAU - stride length of `TAU`
+* @param {NonNegativeInteger} offsetTAU - starting index for `TAU`
+* @param {Float64Array} WORK - workspace array
+* @param {integer} strideWORK - stride length of `WORK`
+* @param {NonNegativeInteger} offsetWORK - starting index for `WORK`
+* @param {NonNegativeInteger} LWORK - dimension of the array `WORK`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1, 0, 0, -1, -2, 0 ] );
+* var TAU = new Float64Array( [ 1, 1 ] );
+* var WORK = new Float64Array( 10 );
+*
+* var info = dorgbr( 2, 3, 2, A, 1, 2, 0, TAU, 1, 0, WORK, 1, 0, 10 );
+* // A => <Float64Array>[ 0, 0, 0, 0, 2, 0 ]
+* // info => 0
+* // WORK[ 0 ] => 2
+*/
+function dorgbr( M, N, K, A, strideA1, strideA2, offsetA, TAU, strideTAU, offsetTAU, WORK, strideWORK, offsetWORK, LWORK ) {
+	var ldwork;
+	var lquery;
+	var lwkopt;
+	var nbmin;
+	var iws;
+	var nb;
+	var nx;
+	var kk;
+	var ki;
+	var ib;
+	var i;
+	var j;
+	var l;
+
+	nb = 32;  // Optimal block size derived from the Fortran LAPACK call `ilaenv( 1, 'DORGLQ', ' ', m, n, k, -1 )`
+	lwkopt = max( 1, M ) * nb;
+	WORK[ offsetWORK ] = lwkopt;
+	lquery = ( LWORK === -1 );
+
+	if ( lquery ) {
+		return 0;
+	}
+
+	// Quick return if possible
+	if ( M === 0 ) {
+		WORK[ offsetWORK ] = 1;
+		return 0;
+	}
+
+	nbmin = 2;
+	nx = 0;
+	iws = M;
+
+	if ( ( nb > 1 ) && ( nb < K ) ) {
+		// Determine when to cross over from blocked to unblocked code.
+		nx = 128; // The crossover point derived from the Fortran LAPACK call `ilaenv( 3, 'DORGLQ', ' ', m, n, k, -1 )`
+
+		if ( nx < K ) {
+			// Determine if workspace is large enough for blocked code.
+			ldwork = M;
+			iws = ldwork * nb;
+			if ( LWORK < iws ) {
+				// Not enough workspace to use optimal NB: reduce NB and determine the minimum value of NB.
+				nb = floor( LWORK / ldwork );
+				nbmin = 2; // The minimum block size derived from the Fortran LAPACK call `ilaenv( 2, 'DORGLQ', ' ', m, n, k, -1 )`
+			}
+		}
+	}
+
+	if ( ( nb >= nbmin ) && ( nb < K ) && ( nx < K ) ) {
+		// Use blocked code after the last block.
+		// The first kk rows are handled by the block method.
+		ki = floor( ( K - nx - 1 ) / nb ) * nb;
+		kk = min( K, ki + nb );
+
+		// Set A(kk+1:m,1:kk) to zero.
+		for ( j = 0; j < kk; j++ ) {
+			for ( i = kk; i < M; i++ ) {
+				A[ offsetA + ( i * strideA1 ) + ( j * strideA2 ) ] = 0;
+			}
+		}
+	} else {
+		kk = 0;
+	}
+
+	// Use unblocked code for the last or only block.
+	if ( kk < M ) {
+		dorgl2( M - kk, N - kk, K - kk, A, strideA1, strideA2, offsetA + ( kk * ( strideA1 + strideA2 ) ), TAU, strideTAU, offsetTAU + ( kk * strideTAU ), WORK, strideWORK, offsetWORK );
+	}
+
+	if ( kk > 0 ) {
+		// Use blocked code
+		for ( i = ki; i >= 0; i -= nb ) {
+			ib = min( nb, K - i );
+
+			if ( i + ib < M ) {
+				// Form the triangular factor of the block reflector, H = H(i) H(i+1) . . . H(i+ib-1)
+				dlarft( 'forward', 'rows', N - i, ib, A, strideA1, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), TAU, strideTAU, offsetTAU + ( i * strideTAU ), WORK, 1, ldwork, offsetWORK );
+
+				// Apply H**T to A(i+ib:m,i:n) from the right
+				dlarfb( 'right', 'transpose', 'forward', 'rows', M - i - ib, N - i, ib, A, strideA1, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), WORK, 1, ldwork, offsetWORK, A, strideA1, strideA2, offsetA + ( ( i + ib ) * strideA1 ) + ( i * strideA2 ), WORK, 1, ldwork, offsetWORK + ( ib * strideWORK ) );
+			}
+
+			// Apply H**T to columns i:n of current block
+			dorgl2( ib, N - i, ib, A, strideA1, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), TAU, strideTAU, offsetTAU + ( i * strideTAU ), WORK, strideWORK, offsetWORK );
+
+			// Set columns 1:i-1 of current block to zero
+			for ( j = 0; j < i; j++ ) {
+				for ( l = i; l < i + ib; l++ ) {
+					A[ offsetA + ( l * strideA1 ) + ( j * strideA2 ) ] = 0.0;
+				}
+			}
+		}
+	}
+
+	WORK[ offsetWORK ] = iws;
+	return 0;
+}
+
+
+// EXPORTS //
+
+module.exports = dorgbr;

lib/node_modules/@stdlib/lapack/base/dorgbr/lib/dlarf1f.js

@@ -0,0 +1,184 @@
+/**
+* @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';
+
+/* eslint-disable max-len */
+
+// MODULES //
+
+var iladlc = require( '@stdlib/lapack/base/iladlc' ).ndarray;
+var iladlr = require( '@stdlib/lapack/base/iladlr' ).ndarray;
+var dgemv = require( '@stdlib/blas/base/dgemv' ).ndarray;
+var dger = require( '@stdlib/blas/base/dger' ).ndarray;
+var daxpy = require( '@stdlib/blas/base/daxpy' ).ndarray;
+var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
+
+
+// FUNCTIONS //
+
+/**
+* Tests whether an operation should be applied to the left side.
+*
+* @private
+* @param {string} side - operation side
+* @returns {boolean} boolean indicating if an operation should be applied to the left side
+*/
+function isLeftSide( side ) {
+	return side === 'left';
+}
+
+
+// MAIN //
+
+/**
+* Applies a real elementary reflector `H = I - tau * v * v^T` to a real M by N matrix `C`.
+*
+* ## Notes
+*
+* -   If `side = 'left'`,
+*
+*     -   `work` should have `N` indexed elements.
+*     -   `V` should have `1 + (M-1) * abs(strideV)` indexed elements.
+*     -   `C` is overwritten by `H * C`.
+*
+* -   If `side = 'right'`,
+*
+*     -   `work` should have `M` indexed elements.
+*     -   `V` should have `1 + (N-1) * abs(strideV)` indexed elements.
+*     -   `C` is overwritten by `C * H`.
+*
+* @private
+* @param {string} side - specifies the side of multiplication with `C`
+* @param {NonNegativeInteger} M - number of rows in `C`
+* @param {NonNegativeInteger} N - number of columns in `C`
+* @param {Float64Array} V - the vector `v`
+* @param {integer} strideV - stride length for `V`
+* @param {NonNegativeInteger} offsetV - starting index for `V`
+* @param {number} tau - scalar constant
+* @param {Float64Array} C - input matrix
+* @param {integer} strideC1 - stride of the first dimension of `C`
+* @param {integer} strideC2 - stride of the second dimension of `C`
+* @param {NonNegativeInteger} offsetC - starting index for `C`
+* @param {Float64Array} work - workspace array
+* @param {integer} strideWork - stride length for `work`
+* @param {NonNegativeInteger} offsetWork - starting index for `work`
+* @returns {Float64Array} `C * H` or `H * C`
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var C = new Float64Array( [ 1.0, 5.0, 9.0, 2.0, 6.0, 10.0, 3.0, 7.0, 11.0, 4.0, 8.0, 12.0 ] );
+* var V = new Float64Array( [ 0.5, 0.5, 0.5, 0.5 ] );
+* var work = new Float64Array( 3 );
+*
+* var out = dlarf1f( 'left', 4, 3, V, 1, 0, 1.0, C, 3, 1, 0, work, 1, 0 );
+* // returns <Float64Array>[ -4.5, -10.5, -16.5, -0.75, -1.75, -2.75, 0.25, -0.75, -1.75, 1.25,  0.25, -0.75 ]
+*/
+function dlarf1f( side, M, N, V, strideV, offsetV, tau, C, strideC1, strideC2, offsetC, work, strideWork, offsetWork ) { // eslint-disable-line max-params
+	var lastv;
+	var lastc;
+	var i;
+
+	if ( tau === 0.0 ) {
+		return C;
+	}
+	lastc = 0;
+	if ( isLeftSide( side ) ) {
+		lastv = M;
+	} else {
+		lastv = N;
+	}
+	// Initialize `i` to point to the last element in `V`:
+	i = offsetV + ( ( lastv-1 ) * strideV );
+
+	// Move `i` to the last non-zero element in `V`, where we assume that V[0] = 1, and it is not stored, so we shouldn't access it...
+	while ( lastv > 1 && V[ i ] === 0.0 ) {
+		lastv -= 1;
+		i -= strideV;
+	}
+	if ( isLeftSide( side ) ) {
+		// Scan for the last non-zero column in `C`:
+		lastc = iladlc( lastv + 1, N, C, strideC1, strideC2, offsetC ) + 1; // adjust by `+1` to account for the difference between zero-based and one-based indexing
+	} else {
+		// Scan for the last non-zero row in `C`:
+		lastc = iladlr( M, lastv + 1, C, strideC1, strideC2, offsetC ) + 1; // // adjust by `+1` to account for the difference between zero-based and one-based indexing
+	}
+	// Return `C` unchanged if all elements in `C` are zero...
+	if ( lastc === 0 ) {
+		return C;
+	}
+	if ( isLeftSide( side ) ) {
+		// Form: H*C
+
+		// If `lastv = 1`, this means `V = 1`, so we just need to compute `C = H*C = (1-tau)*C`...
+		if ( lastv === 1 ) {
+			// C[0,0:lastc] = (1-tau)*C[0,0:lastc]
+			dscal( lastc, 1.0-tau, C, strideC2, offsetC ); // scale the first row
+		} else {
+			// work[0:lastc,0] = C[0:lastv,0:lastc]^T * V[0:lastv,0]
+
+			// work[0:lastc,0] = C[1:lastv,0:lastc]^T * V[1:lastv,0]
+			dgemv( 'transpose', lastv-1, lastc, 1.0, C, strideC1, strideC2, offsetC+strideC1, V, strideV, offsetV+strideV, 0.0, work, strideWork, offsetWork );
+
+			// work[0:lastc,0] += C[0,0:lastc]^T * V[0,0] = C[0,0:lastc]^T
+			daxpy( lastc, 1.0, C, strideC2, offsetC, work, strideWork, offsetWork ); // operates on the first row of C
+
+			// C[0:lastv,0:lastc] = C[...] - ( tau * V[0:lastv,0] * work[0:lastc,0]^T)
+
+			// C[0,0:lastc] = C[...] - ( tau * V[0,0] * work[0:lastc,0]^^T ) = C[...] - ( tau * work[0:lastc,0]^T )
+			daxpy( lastc, -tau, work, strideWork, offsetWork, C, strideC2, offsetC ); // operates on the first row of C
+
+			// C[1:lastv,0:lastc] = C[...] - ( tau * V[1:lastv,0] * work[0:lastc,0]^T )
+			dger( lastv-1, lastc, -tau, V, strideV, offsetV+strideV, work, strideWork, offsetWork, C, strideC1, strideC2, offsetC+strideC1 );
+		}
+		return C;
+	}
+	// side === 'right'
+
+	// Form: C*H
+
+	// If `lastv = 1`, then `V = 1`, so we just need to compute `C = CH = C*(1-tau)`...
+	if ( lastv === 1 ) {
+		// C[0:lastc,0] = ( 1-tau ) * C[0:lastc,0]
+		dscal( lastc, 1.0-tau, C, strideC1, offsetC ); // scale the first column
+		return C;
+	}
+	// work[0:lastc,0] = ( 1-tau ) * C[0:lastc,0]
+
+	// work[0:lastc,0] = C[0:lastc,1:lastv] * V[1:lastv,0]
+	dgemv( 'no-transpose', lastc, lastv-1, 1.0, C, strideC1, strideC2, offsetC+strideC2, V, strideV, offsetV+strideV, 0.0, work, strideWork, offsetWork );
+
+	// work[0:lastc,0] += C[0:lastc,0] * V[0,0] = C[0:lastc,0]
+	daxpy( lastc, 1.0, C, strideC1, offsetC, work, strideWork, offsetWork ); // operates on the first column of C
+
+	// C[0:lastc,0:lastv] = C[...] - ( tau * work[0:lastc,0] * V[0:lastv,0]^T )
+
+	// C[0:lastc,0] = C[...] - ( tau * work[0:lastc,0] * V[0,0]^T ) = C[...] - ( tau * work[0:lastc,0] )
+	daxpy( lastc, -tau, work, strideWork, offsetWork, C, strideC1, offsetC ); // operates on the first column of C
+
+	// C[0:lastc,1:lastv] = C[...] - ( tau * work[0:lastc,0] * V[1:lastv]^T )
+	dger( lastc, lastv-1, -tau, work, strideWork, offsetWork, V, strideV, offsetV+strideV, C, strideC1, strideC2, offsetC+strideC2 );
+
+	return C;
+}
+
+
+// EXPORTS //
+
+module.exports = dlarf1f;