diff --git a/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/base.js b/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/base.js
new file mode 100644
index 000000000000..803057d1a4f9
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/base.js
@@ -0,0 +1,198 @@
+/**
+* @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 dgemv = require( '@stdlib/blas/base/dgemv' ).ndarray;
+var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
+var min = require( '@stdlib/math/base/special/min' );
+var dlarfg = require( './dlarfg.js' );
+
+
+// MAIN //
+
+/**
+* Reduces the first `NB` rows and columns of a real general matrix `A` to upper or lower bi-diagonal form by an orthogonal transformation `Q**T*A*P`.
+*
+* ## Notes
+*
+* -   If `M >= N`,
+*
+*     -   `A` is reduced to upper bi-diagonal form.
+*     -   Elements on and below the diagonal in the first `NB` columns, with the array `TAUQ`, represent the orthogonal matrix `Q` as a product of elementary reflectors.
+*     -   Elements above the diagonal in the first `NB` rows, with the array `TAUP`, represent the orthogonal matrix `P` as a product of elementary reflectors.
+*
+* -   If `M < N`,
+*
+*     -   `A` is reduced to lower bi-diagonal form.
+*     -   Elements below the diagonal in the first `NB` columns, with the array `TAUQ`, represent the orthogonal matrix `Q` as a product of elementary reflectors.
+*     -   Elements on and above the diagonal in the first `NB` rows, with the array `TAUP`, represent the orthogonal matrix `P` as a product of elementary reflectors.
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows in `A`
+* @param {NonNegativeInteger} N - number of columns in `A`
+* @param {integer} NB - number of leading rows and columns of `A` to reduce
+* @param {Float64Array} A - input 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 {Float64Array} D - real diagonal elements (length `NB`)
+* @param {integer} strideD - stride length for `D`
+* @param {NonNegativeInteger} offsetD - starting index for `D`
+* @param {Float64Array} E - real off-diagonal elements (length `NB`)
+* @param {integer} strideE - stride length for `E`
+* @param {NonNegativeInteger} offsetE - starting index for `E`
+* @param {Float64Array} TAUQ - scalars factors of the elementary reflectors that represent the orthogonal matrix `Q` (length `NB`)
+* @param {integer} strideTAUQ - stride length for `TAUQ`
+* @param {NonNegativeInteger} offsetTAUQ - starting index for `TAUQ`
+* @param {Float64Array} TAUP - scalars factors of the elementary reflectors that represent the orthogonal matrix `P` (length `NB`)
+* @param {integer} strideTAUP - stride length for `TAUP`
+* @param {NonNegativeInteger} offsetTAUP - starting index for `TAUP`
+* @param {Float64Array} X - output matrix
+* @param {integer} strideX1 - stride of the first dimension of `X`
+* @param {integer} strideX2 - stride of the second dimension of `X`
+* @param {NonNegativeInteger} offsetX - starting index for `X`
+* @param {Float64Array} Y - output matrix
+* @param {integer} strideY1 - stride of the first dimension of `Y`
+* @param {integer} strideY2 - stride of the second dimension of `Y`
+* @param {NonNegativeInteger} offsetY - starting index for `Y`
+*/
+function dlabrd( M, N, NB, A, strideA1, strideA2, offsetA, D, strideD, offsetD, E, strideE, offsetE, TAUQ, strideTAUQ, offsetTAUQ, TAUP, strideTAUP, offsetTAUP, X, strideX1, strideX2, offsetX, Y, strideY1, strideY2, offsetY ) {
+	var i;
+
+	if ( M <= 0 || N <= 0 ) {
+		return;
+	}
+
+	if ( M >= N ) {
+		// Reduce to upper bi-diagonal form
+		for ( i = 0; i < NB; i++ ) {
+			// Update A(i:M-1, i)
+			dgemv( 'no-transpose', M - i, i, -1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ), Y, strideY2, offsetY + ( i * strideY1 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) );
+
+			dgemv( 'no-transpose', M - i, i, -1.0, X, strideX1, strideX2, offsetX + ( i * strideX1 ), A, strideA1, offsetA + ( i * strideA2 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) );
+
+			// Generate elementary reflector H(i) to annihilate A(i+1:M-1, i)
+			dlarfg( M - i, A, offsetA + ( i * strideA1 ) + ( i * strideA2 ), A, strideA1, offsetA + ( ( i + 1, M - 1 ) * strideA1 ) + ( i * strideA2 ), TAUQ, offsetTAUQ + ( i * strideTAUQ ) );
+			A[ offsetD + ( i * strideD ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ];
+
+			if ( i < N - 1 ) {
+				A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ] = 1.0;
+
+				// Compute Y(i+1:N-1, i)
+				dgemv( 'transpose', M - i, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ) + ( ( i + 1 ) * strideA2 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dgemv( 'transpose', M - i, i, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );
+
+				dgemv( 'no-transpose', N - i - 1, i, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dgemv( 'transpose', M - i, i, 1.0, X, strideX1, strideX2, offsetX + ( i * strideX1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );
+
+				dgemv( 'transpose', i, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dscal( N - i - 1, TAUQ[ offsetTAUQ + ( i * strideTAUQ ) ], Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				// Update A(i, i+1:N-1)
+				dgemv( 'no-transpose', N - i - 1, i + 1, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), A, strideA2, offsetA + ( i * strideA1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 );
+
+				dgemv( 'transpose', i, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), X, strideX2, offsetX + ( i * strideX1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 );
+
+				// Generate elementary reflector P(i) to annihilate A(i, i+2:N-1)
+				dlarfg( N - i - 1, A, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, A, strideA2, offsetA + ( i * strideA1 ) + ( min( i + 2, N - 1 ) * strideA2 ), TAUP, offsetTAUP + ( i * strideTAUP ) );
+				E[ offsetE + ( i * strideE ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 ];
+				A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 ] = 1.0;
+
+				// Compute X(i+1:M-1, i)
+
+				dgemv( 'no-transpose', M - i - 1, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * ( strideA1 + strideA2 ) ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dgemv( 'transpose', N - i - 1, i + 1, 1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * strideX2 ) );
+
+				dgemv( 'no-transpose', M - i - 1, i + 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dgemv( 'no-transpose', i, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * strideX2 ) );
+
+				dgemv( 'no-transpose', M - i - 1, i, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dscal( M - i - 1, TAUP[ offsetTAUP + ( i * strideTAUP ) ], X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+			}
+		}
+	} else {
+		// Reduce to lower bi-diagonal form (M < N)
+		for ( i = 0; i < NB; i++ ) {
+			// Update A(i, i:N-1)
+			dgemv( 'no-transpose', N - i, i, -1.0, Y, strideY1, strideY2, offsetY + ( i * strideY1 ), A, strideA2, offsetA + ( i * strideA1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) );
+
+			dgemv( 'transpose', i, N - i, -1.0, A, strideA1, strideA2, offsetA + ( i * strideA2 ), X, strideX2, offsetX + ( i * strideX1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) );
+
+			// Generate elementary reflector P(i) to annihilate A(i, i+1:N-1)
+			dlarfg( N - i, A, offsetA + ( i * ( strideA1 + strideA2 ) ), A, strideA2, offsetA + ( i * strideA1 ) + ( min( i + 1, N - 1 ) * strideA2 ), TAUP, offsetTAUP + ( i * strideTAUP ) );
+			D[ offsetD + ( i * strideD ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ];
+
+			if ( i < M - 1 ) {
+				// Set A(i,i) = 1 for use as reflector vector
+				A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ] = 1.0;
+
+				// Compute X(i+1:M-1, i)
+				dgemv( 'no-transpose', M - i - 1, N - i, 1.0, A, strideA1, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dgemv( 'transpose', N - i, i, 1.0, Y, strideY1, strideY2, offsetY + ( i * strideY1 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * strideX2 ) );
+
+				dgemv( 'no-transpose', M - i - 1, i, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dgemv( 'no-transpose', i, N - i, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA2 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * strideX2 ) );
+
+				dgemv( 'no-transpose', M - i - 1, i, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				dscal( M - i - 1, TAUP[ offsetTAUP + ( i * strideTAUP ) ], X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
+
+				// Update A(i+1:M-1, i)
+				dgemv( 'no-transpose', M - i - 1, i, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), Y, strideY2, offsetY + ( i * strideY1 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 );
+
+				dgemv( 'no-transpose', M - i - 1, i + 1, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), A, strideA1, offsetA + ( i * strideA2 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 );
+
+				dlarfg( M - i - 1, A, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, A, strideA1, offsetA + ( min( i + 2, M - 1 ) * strideA1 ) + ( i * strideA2 ), TAUQ, offsetTAUQ + ( i * strideTAUQ ) );
+
+				E[ offsetE + ( i * strideE ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 ];
+				A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 ] = 1.0;
+
+				// Compute Y(i+1:N-1, i)
+				dgemv( 'transpose', M - i - 1, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * ( strideA1 + strideA2 ) ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dgemv( 'transpose', M - i - 1, i, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );
+
+				dgemv( 'no-transpose', N - i - 1, i, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dgemv( 'transpose', M - i - 1, i + 1, 1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );
+
+				dgemv( 'transpose', i + 1, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+
+				dscal( N - i - 1, TAUQ[ offsetTAUQ + ( i * strideTAUQ ) ], Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
+			}
+		}
+	}
+}
+
+
+// EXPORTS //
+
+module.exports = dlabrd;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/dlarfg.js b/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/dlarfg.js
new file mode 100644
index 000000000000..04ab125dc7e7
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlabrd/lib/dlarfg.js
@@ -0,0 +1,142 @@
+/**
+* @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 dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray;
+var sign = require( '@stdlib/math/base/special/copysign' );
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
+var dlapy2 = require( '@stdlib/lapack/base/dlapy2' );
+
+
+// MAIN //
+
+/**
+* Generates a real elementary reflector `H` of order `N` such that applying `H` to a vector `[alpha; X]` zeros out `X`.
+*
+* `H` is a Householder matrix with the form:
+*
+* ```tex
+* H \cdot \begin{bmatrix} \alpha \\ x \end{bmatrix} = \begin{bmatrix} \beta \\ 0 \end{bmatrix}, \quad \text{and} \quad H^T H = I
+* ```
+*
+* where:
+*
+* -   `tau` is a scalar
+* -   `X` is a vector of length `N-1`
+* -   `beta` is a scalar value
+* -   `H` is an orthogonal matrix known as a Householder reflector.
+*
+* The reflector `H` is constructed in the form:
+*
+* ```tex
+* H = I - \tau \begin{bmatrix}1 \\ v \end{bmatrix} \begin{bmatrix}1 & v^T \end{bmatrix}
+* ```
+*
+* where:
+*
+* -   `tau` is a real scalar
+* -   `V` is a real vector of length `N-1` that defines the Householder vector
+* -   The vector `[1; V]` is the Householder direction.
+*
+* The values of `tau` and `V` are chosen so that applying `H` to the vector `[alpha; X]` results in a new vector `[beta; 0]`, i.e., only the first component remains nonzero. The reflector matrix `H` is symmetric and orthogonal, satisfying `H^T = H` and `H^T H = I`
+*
+* ## Special cases
+*
+* -   If all elements of `X` are zero, then `tau = 0` and `H = I`, the identity matrix.
+* -   Otherwise, `tau` satisfies `1 ≤ tau ≤ 2`, ensuring numerical stability in transformations.
+*
+* ## Notes
+*
+* -   `X` should have `N-1` indexed elements
+* -   The output array contains the following two elements: `alpha` and `tau`
+*
+* @private
+* @param {NonNegativeInteger} N - number of rows/columns of the elementary reflector `H`
+* @param {Float64Array} X - input vector
+* @param {integer} strideX - stride length for `X`
+* @param {NonNegativeInteger} offsetX - starting index of `X`
+* @param {Float64Array} out - output array
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index of `out`
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var X = new Float64Array( [ 2.0, 3.0, 4.0 ] );
+* var out = new Float64Array( [ 4.0, 0.0 ] );
+*
+* dlarfg( 4, X, 1, 0, out, 1, 0 );
+* // X => <Float64Array>[ ~0.19, ~0.28, ~0.37 ]
+* // out => <Float64Array>[ ~-6.7, ~1.6 ]
+*/
+function dlarfg( N, X, strideX, offsetX, out, strideOut, offsetOut ) {
+	var safemin;
+	var rsafmin;
+	var xnorm;
+	var alpha;
+	var beta;
+	var tau;
+	var knt;
+	var i;
+
+	if ( N <= 1 ) {
+		out[ offsetOut + strideOut ] = 0.0;
+		return;
+	}
+
+	xnorm = dnrm2( N - 1, X, strideX, offsetX );
+	alpha = out[ offsetOut ];
+
+	if ( xnorm === 0.0 ) {
+		out[ strideOut + offsetOut ] = 0.0;
+	} else {
+		beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
+		safemin = dlamch( 'safemin' ) / dlamch( 'epsilon' );
+		knt = 0;
+		if ( abs( beta ) < safemin ) {
+			rsafmin = 1.0 / safemin;
+			while ( abs( beta ) < safemin && knt < 20 ) {
+				knt += 1;
+				dscal( N-1, rsafmin, X, strideX, offsetX );
+				beta *= rsafmin;
+				alpha *= rsafmin;
+			}
+			xnorm = dnrm2( N - 1, X, strideX, offsetX );
+			beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
+		}
+		tau = ( beta - alpha ) / beta;
+		dscal( N-1, 1.0 / ( alpha - beta ), X, strideX, offsetX );
+		for ( i = 0; i < knt; i++ ) {
+			beta *= safemin;
+		}
+
+		out[ offsetOut ] = beta;
+		out[ strideOut + offsetOut ] = tau;
+	}
+}
+
+
+// EXPORTS //
+
+module.exports = dlarfg;