diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/README.md b/lib/node_modules/@stdlib/lapack/base/dlasq1/README.md
new file mode 100644
index 000000000000..f32acf791c13
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/README.md
@@ -0,0 +1,300 @@
+
+
+# dlasq1
+
+> Compute the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E`.
+
+
+
+The `dlasq1` routine computes all singular values of a real N-by-N bidiagonal matrix:
+
+```math
+B = \left[\begin{array}{rrrrr}d_1 & e_1 & 0 & 0 & 0 \\0 & d_2 & e_2 & 0 & 0 \\0 & 0 & d_3 & e_3 & 0 \\0 & 0 & 0 & d_4 & e_4 \\0 & 0 & 0 & 0 & d_5\end{array}\right]
+```
+
+where:
+
+- `d` contains the diagonal elements.
+- `e` contains the off-diagonal elements (`e_N` is ignored on input).
+
+The algorithm is based on the DQDS (Differential Quotient-Difference with Shift) method, which computes singular values to high relative accuracy. On normal exit, the diagonal array `D` is overwritten with the singular values in decreasing order.
+
+
+
+
+
+
+
+## Usage
+
+```javascript
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+```
+
+#### dlasq1( N, D, E, WORK )
+
+Computes the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E`.
+
+```javascript
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+
+var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+var WORK = new Float64Array( 12 );
+
+var info = dlasq1( 3, D, E, WORK );
+// D => [ ~139.377, ~56.064, ~15.997 ]
+// E => [ 90.0, 40.0, 0.0 ]
+// info => 0
+```
+
+The function has the following parameters:
+
+- **N**: number of rows/columns in the matrix.
+- **D**: the diagonal elements of the bidiagonal matrix as a [`Float64Array`][@stdlib/array/float64]. Should have `N` indexed elements. On normal exit, `D` is overwritten with the singular values in decreasing order.
+- **E**: the off-diagonal elements of the bidiagonal matrix as a [`Float64Array`][@stdlib/array/float64]. Should have `N` indexed elements. On entry, elements `E(1:N-1)` contain the off-diagonal elements. On exit, `E` is overwritten.
+- **WORK**: workspace array as a [`Float64Array`][@stdlib/array/float64]. Should have length at least `4*N`.
+
+Note that indexing is relative to the first index. To introduce an offset, use [`typed array`][mdn-typed-array] views.
+
+
+
+```javascript
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+
+// Initial arrays...
+var D0 = new Float64Array( [ 0.0, 100.0, 50.0, 25.0 ] );
+var E0 = new Float64Array( [ 0.0, 90.0, 40.0, 0.0 ] );
+var WORK0 = new Float64Array( 13 );
+
+// Create offset views...
+var D = new Float64Array( D0.buffer, D0.BYTES_PER_ELEMENT*1 ); // start at 2nd element
+var E = new Float64Array( E0.buffer, E0.BYTES_PER_ELEMENT*1 ); // start at 2nd element
+var WORK = new Float64Array( WORK0.buffer, WORK0.BYTES_PER_ELEMENT*1 ); // start at 2nd element
+
+var info = dlasq1( 3, D, E, WORK );
+// D0 => [ 0.0, ~139.377, ~56.064, ~15.997 ]
+// E0 => [ 0.0, 90.0, 40.0, 0.0 ]
+// info => 0
+```
+
+#### dlasq1.ndarray( N, D, sD, oD, E, sE, oE, WORK, sWORK, oWORK )
+
+Computes the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E` using alternative indexing semantics.
+
+```javascript
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+
+var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+var WORK = new Float64Array( 12 );
+
+var info = dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+// D => [ ~139.377, ~56.064, ~15.997 ]
+// E => [ 90.0, 40.0, 0.0 ]
+// info => 0
+```
+
+The function has the following additional parameters:
+
+- **sD**: stride length for `D`.
+- **oD**: starting index for `D`.
+- **sE**: stride length for `E`.
+- **oE**: starting index for `E`.
+- **sWORK**: stride length for `WORK`.
+- **oWORK**: starting index for `WORK`.
+
+While [`typed array`][mdn-typed-array] views mandate a view offset based on the underlying buffer, the offset parameters support indexing semantics based on starting indices. For example,
+
+```javascript
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+
+var D = new Float64Array( [ 0.0, 100.0, 50.0, 25.0 ] );
+var E = new Float64Array( [ 0.0, 90.0, 40.0, 0.0 ] );
+var WORK = new Float64Array( 13 );
+
+var info = dlasq1.ndarray( 3, D, 1, 1, E, 1, 1, WORK, 1, 1 );
+// D => [ 0.0, ~139.377, ~56.064, ~15.997 ]
+// E => [ 0.0, 90.0, 40.0, 0.0 ]
+// info => 0
+```
+
+
+
+
+
+
+
+## Notes
+
+- Both functions mutate the input arrays `D`, `E`, and `WORK`.
+
+- Both functions return a status code indicating success or failure. The status code indicates the following conditions:
+
+ - `0`: successful exit.
+
+ - `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value.
+
+ - `> 0`: the algorithm failed:
+ - `1`: a split was marked by a positive value in `E`.
+ - `2`: current block not diagonalized after `100*N` iterations (in inner while loop). On exit, `D` and `E` represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into `DLASQ1`.
+ - `3`: termination criterion of outer while loop not met (program created more than `N` unreduced blocks).
+
+- `dlasq1()` corresponds to the [LAPACK][LAPACK] routine [`dlasq1`][lapack-dlasq1].
+
+
+
+
+
+
+
+## Examples
+
+
+
+```javascript
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+
+var N = 4;
+
+/*
+ B = [
+ [ 5.0, 4.0, 0.0, 0.0 ],
+ [ 0.0, 2.0, 3.0, 0.0 ],
+ [ 0.0, 0.0, 6.0, 2.0 ],
+ [ 0.0, 0.0, 0.0, 1.0 ]
+ ]
+*/
+
+var D = new Float64Array( [ 5.0, 2.0, 6.0, 1.0 ] );
+var E = new Float64Array( [ 4.0, 3.0, 2.0, 0.0 ] );
+var WORK = new Float64Array( 16 );
+
+// Compute the singular values:
+var info = dlasq1( N, D, E, WORK );
+
+console.log( D );
+
+console.log( info );
+```
+
+
+
+
+
+
+
+* * *
+
+
+
+## C APIs
+
+
+
+
+
+
+
+
+
+
+
+### Usage
+
+```c
+TODO
+```
+
+#### TODO
+
+TODO.
+
+```c
+TODO
+```
+
+TODO
+
+```c
+TODO
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+### Examples
+
+```c
+TODO
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+[lapack]: https://www.netlib.org/lapack/explore-html/
+
+[lapack-dlasq1]: https://www.netlib.org/lapack/explore-html/d5/dce/group__lasq1_ga5a8c1474ef61ff7c59c17412ae456ca6.html
+
+[@stdlib/array/float64]: https://stdlib.io/docs/api/latest/@stdlib/array/float64
+
+[mdn-typed-array]: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/TypedArray
+
+
+
+
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.js
new file mode 100644
index 000000000000..72aada1b00ea
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.js
@@ -0,0 +1,105 @@
+/**
+* @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 bench = require( '@stdlib/bench' );
+var uniform = require( '@stdlib/random/array/uniform' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var isnan = require( '@stdlib/assert/is-nan' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var format = require( '@stdlib/string/format' );
+var pkg = require( './../package.json' ).name;
+var dlasq1 = require( './../lib/dlasq1.js' );
+
+
+// VARIABLES //
+
+var options = {
+ 'dtype': 'float64'
+};
+
+
+// FUNCTIONS //
+
+/**
+* Creates a benchmark function.
+*
+* @private
+* @param {PositiveInteger} len - array length
+* @returns {Function} benchmark function
+*/
+function createBenchmark( len ) {
+ var WORK = uniform( len*4, 0.0, 100.0, options );
+ var D = uniform( len, 0.0, 100.0, options );
+ var E = uniform( len, 0.0, 100.0, options );
+ return benchmark;
+
+ /**
+ * Benchmark function.
+ *
+ * @private
+ * @param {Benchmark} b - benchmark instance
+ */
+ function benchmark( b ) {
+ var i;
+
+ b.tic();
+ for ( i = 0; i < b.iterations; i++ ) {
+ dlasq1( len, D, E, WORK );
+ if ( isnan( D[ i%D.length ] ) ) {
+ b.fail( 'should not return NaN' );
+ }
+ }
+ b.toc();
+ if ( isnan( D[ i%D.length ] ) ) {
+ b.fail( 'should not return NaN' );
+ }
+ b.pass( 'benchmark finished' );
+ b.end();
+ }
+}
+
+
+// MAIN //
+
+/**
+* Main execution sequence.
+*
+* @private
+*/
+function main() {
+ var len;
+ var min;
+ var max;
+ var f;
+ var i;
+
+ min = 1; // 10^min
+ max = 6; // 10^max
+
+ for ( i = min; i <= max; i++ ) {
+ len = floor( pow( pow( 10, i ), 1.0/2.0 ) );
+ f = createBenchmark( len );
+ bench( format( '%s:len=%d', pkg, len ), f );
+ }
+}
+
+main();
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.ndarray.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.ndarray.js
new file mode 100644
index 000000000000..80b50a15cb5c
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/benchmark/benchmark.ndarray.js
@@ -0,0 +1,105 @@
+/**
+* @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 bench = require( '@stdlib/bench' );
+var uniform = require( '@stdlib/random/array/uniform' );
+var isnan = require( '@stdlib/assert/is-nan' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var format = require( '@stdlib/string/format' );
+var pkg = require( './../package.json' ).name;
+var dlasq1 = require( './../lib/ndarray.js' );
+
+
+// VARIABLES //
+
+var options = {
+ 'dtype': 'float64'
+};
+
+
+// FUNCTIONS //
+
+/**
+* Creates a benchmark function.
+*
+* @private
+* @param {PositiveInteger} len - array length
+* @returns {Function} benchmark function
+*/
+function createBenchmark( len ) {
+ var WORK = uniform( len*4, 0.0, 100.0, options );
+ var D = uniform( len, 0.0, 100.0, options );
+ var E = uniform( len, 0.0, 100.0, options );
+ return benchmark;
+
+ /**
+ * Benchmark function.
+ *
+ * @private
+ * @param {Benchmark} b - benchmark instance
+ */
+ function benchmark( b ) {
+ var i;
+
+ b.tic();
+ for ( i = 0; i < b.iterations; i++ ) {
+ dlasq1( len, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+ if ( isnan( D[ i%D.length ] ) ) {
+ b.fail( 'should not return NaN' );
+ }
+ }
+ b.toc();
+ if ( isnan( D[ i%D.length ] ) ) {
+ b.fail( 'should not return NaN' );
+ }
+ b.pass( 'benchmark finished' );
+ b.end();
+ }
+}
+
+
+// MAIN //
+
+/**
+* Main execution sequence.
+*
+* @private
+*/
+function main() {
+ var len;
+ var min;
+ var max;
+ var f;
+ var i;
+
+ min = 1; // 10^min
+ max = 6; // 10^max
+
+ for ( i = min; i <= max; i++ ) {
+ len = floor( pow( pow( 10, i ), 1.0/2.0 ) );
+ f = createBenchmark( len );
+ bench( format( '%s:ndarray:len=%d', pkg, len ), f );
+ }
+}
+
+main();
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/repl.txt b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/repl.txt
new file mode 100644
index 000000000000..792da3c50d99
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/repl.txt
@@ -0,0 +1,143 @@
+
+{{alias}}( N, D, E, WORK )
+ Computes the singular values of a real `N-by-N` bidiagonal matrix with
+ diagonal `D` and off-diagonal `E`.
+
+ Indexing is relative to the first index. To introduce an offset, use typed
+ array views.
+
+ The function mutates `D`, `E`, and `WORK`.
+
+ Parameters
+ ----------
+ N: integer
+ Number of rows/columns in the matrix.
+
+ D: Float64Array
+ The diagonal elements of the bidiagonal matrix. Should have `N`
+ indexed elements. On normal exit, `D` is overwritten with the
+ singular values in decreasing order.
+
+ E: Float64Array
+ The off-diagonal elements of the bidiagonal matrix. Should have `N`
+ indexed elements. On entry, elements E(1:N-1) contain the off-diagonal
+ elements. On exit, `E` is overwritten.
+
+ WORK: Float64Array
+ Workspace array. Should have length at least `4*N`.
+
+ Returns
+ -------
+ info: integer
+ Status code. The status code indicates the following conditions:
+
+ - if equal to zero, then the computation was successful.
+ - if less than zero, then the `i`-th argument where `i = -info` had an
+ illegal value.
+ - if equal to 1, then a split was marked by a positive value in `E`.
+ - if equal to 2, then the current block was not diagonalized after
+ `100*N` iterations.
+ - if equal to 3, then the termination criterion of outer while loop
+ was not met (program created more than `N` unreduced blocks).
+
+ Examples
+ --------
+ > var D = new {{alias:@stdlib/array/float64}}( [ 100.0, 50.0, 25.0 ] );
+ > var E = new {{alias:@stdlib/array/float64}}( [ 90.0, 40.0, 0.0 ] );
+ > var WORK = new {{alias:@stdlib/array/float64}}( 12 );
+ > {{alias}}( 3, D, E, WORK )
+ 0
+ > D
+ [ ~139.377, ~56.064, ~15.997 ]
+ > E
+ [ 90.0, 40.0, 0.0 ]
+
+ // Using typed array views:
+ > var D0 = new {{alias:@stdlib/array/float64}}( [ 0.0, 100.0, 50.0, 25.0 ] );
+ > var E0 = new {{alias:@stdlib/array/float64}}( [ 0.0, 90.0, 40.0, 0.0 ] );
+ > var WORK0 = new {{alias:@stdlib/array/float64}}( 13 );
+ > D = new Float64Array( D0.buffer, D0.BYTES_PER_ELEMENT*1 );
+ > E = new Float64Array( E0.buffer, E0.BYTES_PER_ELEMENT*1 );
+ > WORK = new Float64Array( WORK0.buffer, WORK0.BYTES_PER_ELEMENT*1 );
+ > {{alias}}( 3, D, E, WORK )
+ 0
+ > D0
+ [ 0.0, ~139.377, ~56.064, ~15.997 ]
+ > E0
+ [ 0.0, 90.0, 40.0, 0.0 ]
+
+
+{{alias}}.ndarray( N, D, sd, od, E, se, oe, WORK, sw, ow )
+ Computes the singular values of a real `N-by-N` bidiagonal matrix with
+ diagonal `D` and off-diagonal `E` using alternative indexing semantics.
+
+ While typed array views mandate a view offset based on the underlying
+ buffer, the offset parameters support indexing semantics based on starting
+ indices.
+
+ The function mutates `D`, `E`, and `WORK`.
+
+ Parameters
+ ----------
+ N: integer
+ Number of rows/columns in the matrix.
+
+ D: Float64Array
+ The diagonal elements of the bidiagonal matrix. Should have `N`
+ indexed elements. On normal exit, `D` is overwritten with the
+ singular values in decreasing order.
+
+ sd: integer
+ Stride length for `D`.
+
+ od: integer
+ Starting index for `D`.
+
+ E: Float64Array
+ The off-diagonal elements of the bidiagonal matrix. Should have `N`
+ indexed elements. On entry, elements E(1:N-1) contain the off-diagonal
+ elements. On exit, `E` is overwritten.
+
+ se: integer
+ Stride length for `E`.
+
+ oe: integer
+ Starting index for `E`.
+
+ WORK: Float64Array
+ Workspace array. Should have length at least `4*N`.
+
+ sw: integer
+ Stride length for `WORK`.
+
+ ow: integer
+ Starting index for `WORK`.
+
+ Returns
+ -------
+ info: integer
+ Status code. The status code indicates the following conditions:
+
+ - if equal to zero, then the computation was successful.
+ - if less than zero, then the `i`-th argument where `i = -info` had an
+ illegal value.
+ - if equal to 1, then a split was marked by a positive value in `E`.
+ - if equal to 2, then the current block was not diagonalized after
+ `100*N` iterations.
+ - if equal to 3, then the termination criterion of outer while loop
+ was not met (program created more than `N` unreduced blocks).
+
+ Examples
+ --------
+ > var D = new {{alias:@stdlib/array/float64}}( [ 100.0, 50.0, 25.0 ] );
+ > var E = new {{alias:@stdlib/array/float64}}( [ 90.0, 40.0, 0.0 ] );
+ > var WORK = new {{alias:@stdlib/array/float64}}( 12 );
+ > {{alias}}.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 )
+ 0
+ > D
+ [ ~139.377, ~56.064, ~15.997 ]
+ > E
+ [ 90.0, 40.0, 0.0 ]
+
+ See Also
+ --------
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/index.d.ts b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/index.d.ts
new file mode 100644
index 000000000000..90a34ce7e905
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/index.d.ts
@@ -0,0 +1,145 @@
+/*
+* @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.
+*/
+
+// TypeScript Version: 4.1
+
+/**
+* Status code.
+*
+* ## Notes
+*
+* The status code indicates the following conditions:
+*
+* - if equal to zero, then the computation was successful.
+* - if less than zero, then the `i`-th argument where `i = -INFO` had an illegal value.
+* - if equal to 1, then a split was marked by a positive value in `E`.
+* - if equal to 2, then the current block was not diagonalized after `100*N` iterations.
+* - if equal to 3, then the termination criterion of outer while loop was not met (program created more than `N` unreduced blocks).
+*/
+type StatusCode = number;
+
+/**
+* Interface describing `dlasq1`.
+*/
+interface Routine {
+ /**
+ * Computes the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E`.
+ *
+ * ## Notes
+ *
+ * - `D` should have `N` indexed elements and on normal exit is overwritten with the singular values in decreasing order.
+ * - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements. On exit, `E` is overwritten.
+ * - `WORK` should have length at least `4*N`.
+ *
+ * @param N - number of rows/columns in the matrix
+ * @param D - the diagonal elements of the bidiagonal matrix
+ * @param E - the off-diagonal elements of the bidiagonal matrix
+ * @param WORK - workspace array
+ * @returns status code
+ *
+ * @example
+ * var Float64Array = require( '@stdlib/array/float64' );
+ *
+ * var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+ * var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+ * var WORK = new Float64Array( 12 );
+ *
+ * dlasq1( 3, D, E, WORK );
+ * // D => [ ~139.377, ~56.064, ~15.997 ]
+ * // E => [ 90.0, 40.0, 0.0 ]
+ */
+ ( N: number, D: Float64Array, E: Float64Array, WORK: Float64Array ): StatusCode;
+
+ /**
+ * Computes the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E` using alternative indexing semantics.
+ *
+ * ## Notes
+ *
+ * - `D` should have `N` indexed elements and on normal exit is overwritten with the singular values in decreasing order.
+ * - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements. On exit, `E` is overwritten.
+ * - `WORK` should have length at least `4*N`.
+ *
+ * @param N - number of rows/columns in the matrix
+ * @param D - the diagonal elements of the bidiagonal matrix
+ * @param strideD - stride length for `D`
+ * @param offsetD - starting index of `D`
+ * @param E - the off-diagonal elements of the bidiagonal matrix
+ * @param strideE - stride length for `E`
+ * @param offsetE - starting index of `E`
+ * @param WORK - workspace array
+ * @param strideWORK - stride length for `WORK`
+ * @param offsetWORK - starting index of `WORK`
+ * @returns status code
+ *
+ * @example
+ * var Float64Array = require( '@stdlib/array/float64' );
+ *
+ * var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+ * var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+ * var WORK = new Float64Array( 12 );
+ *
+ * dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+ * // D => [ ~139.377, ~56.064, ~15.997 ]
+ * // E => [ 90.0, 40.0, 0.0 ]
+ */
+ ndarray( N: number, D: Float64Array, strideD: number, offsetD: number, E: Float64Array, strideE: number, offsetE: number, WORK: Float64Array, strideWORK: number, offsetWORK: number ): StatusCode;
+}
+
+/**
+* LAPACK routine to compute the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E`.
+*
+* ## Notes
+*
+* - `D` should have `N` indexed elements and on normal exit is overwritten with the singular values in decreasing order.
+* - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements. On exit, `E` is overwritten.
+* - `WORK` should have length at least `4*N`.
+*
+* @param N - number of rows/columns in the matrix
+* @param D - the diagonal elements of the bidiagonal matrix
+* @param E - the off-diagonal elements of the bidiagonal matrix
+* @param WORK - workspace array
+* @returns status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+* var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* dlasq1( 3, D, E, WORK );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90.0, 40.0, 0.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+* var E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90.0, 40.0, 0.0 ]
+*/
+declare var dlasq1: Routine;
+
+
+// EXPORTS //
+
+export = dlasq1;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/test.ts b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/test.ts
new file mode 100644
index 000000000000..57d46054de4a
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/docs/types/test.ts
@@ -0,0 +1,298 @@
+/*
+* @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.
+*/
+
+// TypeScript Version: 4.1
+
+import dlasq1 = require( './index' );
+
+
+// TESTS //
+
+// The function returns a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1( 3, D, E, WORK ); // $ExpectType number
+}
+
+// The compiler throws an error if the function is provided a first argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1( '5', D, E, WORK ); // $ExpectError
+ dlasq1( true, D, E, WORK ); // $ExpectError
+ dlasq1( false, D, E, WORK ); // $ExpectError
+ dlasq1( null, D, E, WORK ); // $ExpectError
+ dlasq1( undefined, D, E, WORK ); // $ExpectError
+ dlasq1( [], D, E, WORK ); // $ExpectError
+ dlasq1( {}, D, E, WORK ); // $ExpectError
+ dlasq1( ( x: number ): number => x, D, E, WORK ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a second argument which is not a Float64Array...
+{
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1( 3, '5', E, WORK ); // $ExpectError
+ dlasq1( 3, 5, E, WORK ); // $ExpectError
+ dlasq1( 3, true, E, WORK ); // $ExpectError
+ dlasq1( 3, false, E, WORK ); // $ExpectError
+ dlasq1( 3, null, E, WORK ); // $ExpectError
+ dlasq1( 3, undefined, E, WORK ); // $ExpectError
+ dlasq1( 3, [], E, WORK ); // $ExpectError
+ dlasq1( 3, {}, E, WORK ); // $ExpectError
+ dlasq1( 3, ( x: number ): number => x, E, WORK ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a third argument which is not a Float64Array...
+{
+ const D = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1( 3, D, '5', WORK ); // $ExpectError
+ dlasq1( 3, D, 5, WORK ); // $ExpectError
+ dlasq1( 3, D, true, WORK ); // $ExpectError
+ dlasq1( 3, D, false, WORK ); // $ExpectError
+ dlasq1( 3, D, null, WORK ); // $ExpectError
+ dlasq1( 3, D, undefined, WORK ); // $ExpectError
+ dlasq1( 3, D, [], WORK ); // $ExpectError
+ dlasq1( 3, D, {}, WORK ); // $ExpectError
+ dlasq1( 3, D, ( x: number ): number => x, WORK ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a fourth argument which is not a Float64Array...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+
+ dlasq1( 3, D, E, '5' ); // $ExpectError
+ dlasq1( 3, D, E, 5 ); // $ExpectError
+ dlasq1( 3, D, E, true ); // $ExpectError
+ dlasq1( 3, D, E, false ); // $ExpectError
+ dlasq1( 3, D, E, null ); // $ExpectError
+ dlasq1( 3, D, E, undefined ); // $ExpectError
+ dlasq1( 3, D, E, [] ); // $ExpectError
+ dlasq1( 3, D, E, {} ); // $ExpectError
+ dlasq1( 3, D, E, ( x: number ): number => x ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided an unsupported number of arguments...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1(); // $ExpectError
+ dlasq1( 3 ); // $ExpectError
+ dlasq1( 3, D ); // $ExpectError
+ dlasq1( 3, D, E ); // $ExpectError
+ dlasq1( 3, D, E, WORK, 10 ); // $ExpectError
+}
+
+// Attached to main export is an `ndarray` method which returns a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectType number
+}
+
+// The compiler throws an error if the function is provided a first argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( '5', D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( true, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( false, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( null, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( undefined, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( [], D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( {}, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( ( x: number ): number => x, D, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a second argument which is not a Float64Array...
+{
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, '5', 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, 5, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, true, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, false, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, null, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, undefined, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, [], 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, {}, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, ( x: number ): number => x, 1, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a third argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, '5', 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, true, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, false, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, null, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, undefined, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, [], 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, {}, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, ( x: number ): number => x, 0, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a fourth argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, '5', E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, true, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, false, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, null, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, undefined, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, [], E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, {}, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, ( x: number ): number => x, E, 1, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a fifth argument which is not a Float64Array...
+{
+ const D = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, '5', 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, 5, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, true, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, false, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, null, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, undefined, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, [], 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, {}, 1, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, ( x: number ): number => x, 1, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a sixth argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, '5', 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, true, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, false, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, null, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, undefined, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, [], 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, {}, 0, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, ( x: number ): number => x, 0, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a seventh argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, '5', WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, true, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, false, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, null, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, undefined, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, [], WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, {}, WORK, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, ( x: number ): number => x, WORK, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided an eighth argument which is not a Float64Array...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, '5', 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, 5, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, true, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, false, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, null, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, undefined, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, [], 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, {}, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, ( x: number ): number => x, 1, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a ninth argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, '5', 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, true, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, false, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, null, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, undefined, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, [], 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, {}, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, ( x: number ): number => x, 0 ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided a tenth argument which is not a number...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, '5' ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, true ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, false ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, null ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, undefined ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, [] ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, {} ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, ( x: number ): number => x ); // $ExpectError
+}
+
+// The compiler throws an error if the function is provided an unsupported number of arguments...
+{
+ const D = new Float64Array( 3 );
+ const E = new Float64Array( 3 );
+ const WORK = new Float64Array( 12 );
+
+ dlasq1.ndarray(); // $ExpectError
+ dlasq1.ndarray( 3 ); // $ExpectError
+ dlasq1.ndarray( 3, D ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1 ); // $ExpectError
+ dlasq1.ndarray( 3, D, 1, 0, E, 1, 0, WORK, 1, 0, 10 ); // $ExpectError
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/examples/index.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/examples/index.js
new file mode 100644
index 000000000000..e87e246ad431
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/examples/index.js
@@ -0,0 +1,44 @@
+/**
+* @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';
+
+var Float64Array = require( '@stdlib/array/float64' );
+var dlasq1 = require( './../lib' );
+
+var N = 4;
+
+/*
+ B = [
+ [ 5.0, 4.0, 0.0, 0.0 ],
+ [ 0.0, 2.0, 3.0, 0.0 ],
+ [ 0.0, 0.0, 6.0, 2.0 ],
+ [ 0.0, 0.0, 0.0, 1.0 ]
+ ]
+*/
+
+var D = new Float64Array( [ 5.0, 2.0, 6.0, 1.0 ] );
+var E = new Float64Array( [ 4.0, 3.0, 2.0, 0.0 ] );
+var WORK = new Float64Array( 16 );
+
+// Compute the singular values:
+var info = dlasq1( N, D, E, WORK );
+
+console.log( D );
+
+console.log( info );
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/base.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/base.js
new file mode 100644
index 000000000000..30156236b983
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/base.js
@@ -0,0 +1,190 @@
+/**
+* @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 Float64Array = require( '@stdlib/array/float64' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var max = require( '@stdlib/math/base/special/max' );
+var sqrt = require( '@stdlib/math/base/special/sqrt' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var dcopy = require( '@stdlib/blas/base/dcopy' ).ndarray;
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var dlas2 = require( './dlas2.js' );
+var dlascl = require( './dlascl.js' );
+var dlasq2 = require( './dlasq2.js' );
+var dlasrt = require( './dlasrt.js' );
+
+
+// VARIABLES //
+
+var EPS = dlamch( 'Precision' );
+var SAFMIN = dlamch( 'Safe minimum' );
+var SCALE = sqrt( EPS / SAFMIN );
+
+
+// MAIN //
+
+/**
+* Computes the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E`.
+*
+* ## Notes
+*
+* - `D` should have `N` indexed elements. On entry, `D` contains the diagonal elements of the bi-diagonal matrix whose SVD is desired. On normal exit, `D` contains the singular values in decreasing order.
+*
+* - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements of the bi-diagonal matrix whose SVD is desired. On exit, E is overwritten.
+*
+* - The function returns a status code:
+*
+* - `= 0`: successful exit.
+*
+* - `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value.
+*
+* - `> 0`: the algorithm failed:
+* - `= 1`, a split was marked by a positive value in `E`.
+* - `= 2`, current block of `Z` not diagonalized after `100*N` iterations (in inner while loop). On exit, `D` and `E` represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into `DLASQ1`
+* - `= 3`, termination criterion of outer while loop not met (program created more than `N` unreduced blocks)
+*
+* @private
+* @param {integer} N - number of rows/columns in the matrix
+* @param {Float64Array} D - the array with diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {integer} strideD - stride length for `D`
+* @param {NonNegativeInteger} offsetD - starting index of `D`
+* @param {Float64Array} E - the array with off-diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {integer} strideE - stride length for `E`
+* @param {NonNegativeInteger} offsetE - starting index of `E`
+* @param {Float64Array} WORK - workspace array (length >= 4*N)
+* @param {integer} strideWORK - stride length for `WORK`
+* @param {NonNegativeInteger} offsetWORK - starting index of `WORK`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 100, 50, 25 ] );
+* var E = new Float64Array( [ 90, 40, 0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* var info = dlasq1( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90, 40, 0 ]
+* // WORK => [ ~1.939e292, ~3.137e291, ~2.554e290, ~1.575e286, ~3.137e291, ~2.881e291, ~2.278e292, ~2.278e292, 7.0, ~2.667, 0.0, ~2.881e291 ]
+* // info => 0
+*/
+function dlasq1( N, D, strideD, offsetD, E, strideE, offsetE, WORK, strideWORK, offsetWORK ) { // eslint-disable-line stdlib/jsdoc-doctest-decimal-point
+ var sigmx;
+ var sigmn;
+ var info;
+ var out;
+ var id;
+ var ie;
+ var iw;
+ var i;
+
+ if ( N === 0 ) {
+ return 0;
+ }
+ if ( N === 1 ) {
+ D[ offsetD ] = abs( D[ offsetD ] );
+ return 0;
+ }
+ if ( N === 2 ) {
+ out = new Float64Array( 2 );
+ dlas2( D[ offsetD ], E[ offsetE ], D[ offsetD + strideD ], out, 1, 0 );
+ sigmn = out[ 0 ];
+ sigmx = out[ 1 ];
+ D[ offsetD ] = sigmx;
+ D[ offsetD + strideD ] = sigmn;
+ return 0;
+ }
+
+ // Estimate the largest singular value
+ sigmx = 0;
+ id = offsetD;
+ ie = offsetE;
+ for ( i = 0; i < N - 1; i++ ) {
+ D[ id ] = abs( D[ id ] );
+ sigmx = max( sigmx, abs( E[ ie ] ) );
+ id += strideD;
+ ie += strideE;
+ }
+ D[ id ] = abs( D[ id ] );
+
+ // Early return if SIGMX is zero (matrix is already diagonal)
+ if ( sigmx === 0 ) {
+ dlasrt( 'D', N, D, strideD, offsetD );
+ return 0;
+ }
+
+ id = offsetD;
+ for ( i = 0; i < N; i++ ) {
+ sigmx = max( sigmx, D[ id ] );
+ id += strideD;
+ }
+
+ // Copy D and E into WORK (in the Z format) and SCALE (squaring the input data makes scaling by a power of the radix pointless)
+ dcopy( N, D, strideD, offsetD, WORK, 2 * strideWORK, offsetWORK );
+ dcopy( N - 1, E, strideE, offsetE, WORK, 2 * strideWORK, offsetWORK + strideWORK );
+ dlascl( 'G', 0, 0, sigmx, SCALE, ( 2 * N ) - 1, 1, WORK, strideWORK, ( ( 2 * N ) - 1 ) * strideWORK, offsetWORK );
+
+ // Compute the q's and e's.
+ iw = offsetWORK;
+ for ( i = 0; i < ( 2 * N ) - 1; i++ ) {
+ WORK[ iw ] = pow( WORK[ iw ], 2 );
+ iw += strideWORK;
+ }
+ WORK[ offsetWORK + ( ( ( 2 * N ) - 1 ) * strideWORK ) ] = 0;
+
+ // Call DLASQ2 to compute eigenvalues of the qd array
+ info = dlasq2( N, WORK, strideWORK, offsetWORK );
+
+ if ( info === 0 ) {
+ iw = offsetWORK;
+ id = offsetD;
+ for ( i = 0; i < N; i++ ) {
+ D[ id ] = sqrt( WORK[ iw ] );
+ id += strideD;
+ iw += strideWORK;
+ }
+ dlascl( 'G', 0, 0, SCALE, sigmx, N, 1, D, strideD, N * strideD, offsetD );
+ } else if ( info === 2 ) {
+ // Maximum number of iterations exceeded. Move data from WORK into D and E so the calling subroutine can try to finish
+ id = offsetD;
+ ie = offsetE;
+ iw = offsetWORK;
+ for ( i = 0; i < N; i++ ) {
+ D[ id ] = sqrt( WORK[ iw ] );
+ E[ ie ] = sqrt( WORK[ iw + strideWORK ] );
+ id += strideD;
+ ie += strideE;
+ iw += 2 * strideWORK;
+ }
+ dlascl( 'G', 0, 0, SCALE, sigmx, N, 1, D, strideD, N * strideD, offsetD );
+ dlascl( 'G', 0, 0, SCALE, sigmx, N, 1, E, strideE, N * strideE, offsetE );
+ }
+ return info;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq1;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlas2.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlas2.js
new file mode 100644
index 000000000000..87d9af7d3ddd
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlas2.js
@@ -0,0 +1,106 @@
+/**
+* @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 sqrt = require( '@stdlib/math/base/special/sqrt' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+var abs2 = require( '@stdlib/math/base/special/abs2' );
+
+
+// MAIN //
+
+/**
+* Computes the singular values of a `2x2` upper triangular matrix.
+*
+* @private
+* @param {number} F - the (0,0) element of a `2x2` matrix
+* @param {number} G - the (0,1) element of a `2x2` matrix
+* @param {number} H - the (1,1) element of a `2x2` matrix
+* @param {Float64Array} out - output array containing the smaller and larger singular values respectively
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index of `out`
+* @returns {Float64Array} output array
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var out = new Float64Array( 2 );
+* dlas2( 1.0, 2.0, 3.0, out, 1, 0 );
+* // out => [ ~0.822, ~3.65 ]
+*/
+function dlas2( F, G, H, out, strideOut, offsetOut ) {
+ var ssmax;
+ var ssmin;
+ var fhmn;
+ var fhmx;
+ var as;
+ var at;
+ var au;
+ var fa;
+ var ga;
+ var ha;
+ var c;
+
+ fa = abs( F );
+ ga = abs( G );
+ ha = abs( H );
+ fhmn = min( fa, ha );
+ fhmx = max( fa, ha );
+ if ( fhmn === 0.0 ) {
+ ssmin = 0.0;
+ if ( fhmx === 0.0 ) {
+ ssmax = ga;
+ } else {
+ ssmax = max( fhmx, ga ) * sqrt( 1.0 + abs2( min( fhmx, ga ) / max( fhmx, ga ) ) );
+ }
+ } else if ( ga < fhmx ) {
+ as = 1.0 + ( fhmn / fhmx );
+ at = ( fhmx - fhmn ) / fhmx;
+ au = abs2( ga / fhmx );
+ c = 2.0 / ( sqrt( ( as * as ) + au ) + sqrt( ( at * at ) + au ) );
+ ssmin = fhmn * c;
+ ssmax = fhmx / c;
+ } else {
+ au = fhmx / ga;
+ if ( au === 0.0 ) {
+ // Avoid possible harmful underflow if exponent range asymmetric.
+ ssmin = ( fhmn * fhmx ) / ga;
+ ssmax = ga;
+ } else {
+ as = 1.0 + ( fhmn / fhmx );
+ at = ( fhmx - fhmn ) / fhmx;
+ c = 1.0 / ( sqrt( 1.0 + abs2( as * au ) ) + sqrt( 1.0 + abs2( at * au ) ) );
+ ssmin = ( fhmn * c ) * au;
+ ssmin += ssmin;
+ ssmax = ga / ( c + c );
+ }
+ }
+ out[ offsetOut ] = ssmin;
+ out[ offsetOut + strideOut ] = ssmax;
+ return out;
+}
+
+
+// EXPORTS //
+
+module.exports = dlas2;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlascl.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlascl.js
new file mode 100644
index 000000000000..f4f79e971044
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlascl.js
@@ -0,0 +1,571 @@
+/**
+* @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 */
+
+'use strict';
+
+// MODULES //
+
+var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var loopOrder = require( '@stdlib/ndarray/base/nullary-loop-interchange-order' );
+var max = require( '@stdlib/math/base/special/fast/max' );
+var min = require( '@stdlib/math/base/special/fast/min' );
+
+
+// VARIABLES //
+
+var TINY = dlamch( 'safe minimum' );
+var HUGE = 1.0 / TINY;
+
+
+// FUNCTIONS //
+
+/**
+* Multiplies a double-precision floating-point M-by-N general matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 4.0, 2.0, 5.0, 3.0, 6.0 ] );
+*
+* scaleGeneral( 3, 2, A, 2, 1, 0, 2.0 );
+* // A => [ 2.0, 8.0, 4.0, 10.0, 6.0, 12.0 ]
+*/
+function scaleGeneral( M, N, A, strideA1, strideA2, offsetA, mul ) {
+ var da1;
+ var da0;
+ var S1;
+ var S0;
+ var ia;
+ var i0;
+ var i1;
+ var o;
+
+ // Resolve the loop interchange order:
+ o = loopOrder( [ M, N ], [ strideA1, strideA2 ] );
+ S0 = o.sh[ 0 ];
+ S1 = o.sh[ 1 ];
+ da0 = o.sx[ 0 ];
+ da1 = o.sx[ 1 ] - ( S0 * o.sx[ 0 ] );
+
+ // Set the pointer to the first indexed element:
+ ia = offsetA;
+
+ // Iterate over the matrix dimensions...
+ for ( i1 = 0; i1 < S1; i1++ ) {
+ for ( i0 = 0; i0 < S0; i0++ ) {
+ A[ ia ] *= mul;
+ ia += da0;
+ }
+ ia += da1;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N upper triangular matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 2.0, 4.0, 0.0, 3.0, 5.0, 0.0, 0.0, 6.0 ] );
+*
+* scaleUpper( 3, 3, A, 3, 1, 0, true, 2.0 );
+* // A => [ 2.0, 4.0, 8.0, 0.0, 6.0, 10.0, 0.0, 0.0, 12.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 0.0, 0.0, 2.0, 3.0, 0.0, 4.0, 5.0, 6.0 ] );
+*
+* scaleUpper( 3, 3, A, 1, 3, 0, false, 2.0 );
+* // A => [ 2.0, 0.0, 0.0, 4.0, 6.0, 0.0, 8.0, 10.0, 12.0 ]
+*/
+function scaleUpper( M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var idx;
+ var ia;
+ var i0;
+ var i1;
+
+ ia = offsetA;
+ if ( isrm ) {
+ for ( i1 = 0; i1 < M; i1++ ) {
+ idx = ia + ( i1*strideA2 );
+ for ( i0 = i1; i0 < N; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA2;
+ }
+ ia += strideA1;
+ }
+ return;
+ }
+ for ( i1 = 0; i1 < N; i1++ ) {
+ idx = ia;
+ for ( i0 = 0; i0 <= min( i1, M-1 ); i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N lower triangular matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 0.0, 0.0, 2.0, 4.0, 0.0, 3.0, 5.0, 6.0 ] );
+*
+* scaleLower( 3, 3, A, 3, 1, 0, true, 2.0 );
+* // A => [ 2.0, 0.0, 0.0, 4.0, 8.0, 0.0, 6.0, 10.0, 12.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 2.0, 3.0, 0.0, 4.0, 5.0, 0.0, 0.0, 6.0 ] );
+*
+* scaleLower( 3, 3, A, 1, 3, 0, false, 2.0 );
+* // A => [ 2.0, 4.0, 6.0, 0.0, 8.0, 10.0, 0.0, 0.0, 12.0 ]
+*/
+function scaleLower( M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var idx;
+ var ia;
+ var i0;
+ var i1;
+
+ ia = offsetA;
+ if ( isrm ) {
+ for ( i1 = 0; i1 < M; i1++ ) {
+ idx = ia;
+ for ( i0 = 0; i0 <= min( i1, N-1 ); i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA2;
+ }
+ ia += strideA1;
+ }
+ return;
+ }
+ for ( i1 = 0; i1 < N; i1++ ) {
+ idx = ia + ( i1*strideA1 );
+ for ( i0 = i1; i0 < M; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N upper Hessenberg matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 0.0, 9.0, 10.0, 11.0, 0.0, 0.0, 12.0, 13.0 ] );
+*
+* scaleUpperHessenberg( 4, 4, A, 4, 1, 0, true, 2.0 );
+* // A => [ 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 0.0, 18.0, 20.0, 22.0, 0.0, 0.0, 24.0, 26.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 5.0, 0.0, 0.0, 2.0, 6.0, 9.0, 0.0, 3.0, 7.0, 10.0, 12.0, 4.0, 8.0, 11.0, 13.0 ] );
+*
+* scaleUpperHessenberg( 4, 4, A, 1, 4, 0, false, 2.0 );
+* // A => [ 2.0, 10.0, 0.0, 0.0, 4.0, 12.0, 18.0, 0.0, 6.0, 14.0, 20.0, 24.0, 8.0, 16.0, 22.0, 26.0 ]
+*/
+function scaleUpperHessenberg( M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var idx;
+ var ia;
+ var i0;
+ var i1;
+
+ if ( isrm ) {
+ ia = offsetA;
+ for ( i1 = 0; i1 < M; i1++ ) {
+ idx = ia + ( max( i1-1, 0 ) * strideA2 );
+ for ( i0 = max( i1-1, 0 ); i0 < N; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA2;
+ }
+ ia += strideA1;
+ }
+ return;
+ }
+ ia = offsetA;
+ for ( i0 = 0; i0 < N; i0++ ) {
+ idx = ia;
+ for ( i1 = 0; i1 <= min( i0+1, M-1 ); i1++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N symmetric banded lower matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} KL - lower bandwidth of `A` (i.e., number of sub-diagonals)
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.1, 2.2, 3.3, 4.4, 5.5, 6.1, 7.2, 8.3, 9.4, 0.0, 10.1, 11.2, 12.3, 0.0, 0.0 ] );
+*
+* scaleSymmetricBandedLower( 2, 5, 5, A, 5, 1, 0, true, 10.0 );
+* // A => [ 11.0, 22.0, 33.0, 44.0, 55.0, 61.0, 72.0, 83.0, 94.0, 0.0, 101.0, 112.0, 123.0, 0.0, 0.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.1, 6.1, 10.1, 2.2, 7.2, 11.2, 3.3, 8.3, 12.3, 4.4, 9.4, 0.0, 5.5, 0.0, 0.0 ] );
+*
+* scaleSymmetricBandedLower( 2, 5, 5, A, 1, 3, 0, false, 10.0 );
+* // A => [ 11.0, 61.0, 101.0, 22.0, 72.0, 112.0, 33.0, 83.0, 123.0, 44.0, 94.0, 0.0, 55.0, 0.0, 0.0 ]
+*/
+function scaleSymmetricBandedLower( KL, M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var idx;
+ var ia;
+ var i0;
+ var i1;
+ var k3;
+ var k4;
+
+ ia = offsetA;
+ k3 = KL + 1;
+ k4 = N;
+ if ( isrm ) {
+ for ( i1 = 0; i1 < M; i1++ ) {
+ idx = ia;
+ for ( i0 = 0; i0 < N-i1; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA2;
+ }
+ ia += strideA1;
+ }
+ return;
+ }
+ for ( i1 = 0; i1 < N; i1++ ) {
+ idx = ia;
+ for ( i0 = 0; i0 < min( k3, k4-i1 ); i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N symmetric banded upper matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} KU - upper bandwidth of `A`
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 0.0, 0.0, 10.1, 11.2, 12.3, 0.0, 6.1, 7.2, 8.3, 9.4, 1.1, 2.2, 3.3, 4.4, 5.5 ] );
+*
+* scaleSymmetricBandedUpper( 2, 5, 5, A, 5, 1, 0, true, 10.0 );
+* // A => [ 0.0, 0.0, 101.0, 112.0, 123.0, 0.0, 61.0, 72.0, 83.0, 94.0, 11.0, 22.0, 33.0, 44.0, 55.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 0.0, 0.0, 1.1, 0.0, 6.1, 2.2, 10.1, 7.2, 3.3, 11.2, 8.3, 4.4, 12.3, 9.4, 5.5 ] );
+*
+* scaleSymmetricBandedUpper( 2, 5, 5, A, 1, 3, 0, false, 10.0 );
+* // A => [ 0.0, 0.0, 11.0, 0.0, 61.0, 22.0, 101.0, 72.0, 33.0, 112.0, 83.0, 44.0, 123.0, 94.0, 55.0 ]
+*/
+function scaleSymmetricBandedUpper( KU, M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var idx;
+ var ia;
+ var i0;
+ var i1;
+ var j;
+
+ ia = offsetA;
+ if ( isrm ) {
+ for ( i1 = 0; i1 <= KU; i1++ ) {
+ idx = ia + ( max( KU-i1, 0 ) * strideA2 );
+ for ( i0 = max( KU-i1, 0 ); i0 < N; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA2;
+ }
+ ia += strideA1;
+ }
+ return;
+ }
+ for ( i1 = 0; i1 < N; i1++ ) {
+ j = max( KU-i1, 0 );
+ idx = ia + ( j * strideA1 );
+ for ( i0 = j; i0 <= KU; i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+/**
+* Multiplies a double-precision floating-point M-by-N banded matrix `A` by a double-precision floating-point scalar `mul`.
+*
+* @private
+* @param {NonNegativeInteger} KL - lower bandwidth of `A`
+* @param {NonNegativeInteger} KU - upper bandwidth of `A`
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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 {boolean} isrm - boolean indicating if the matrix is row-major
+* @param {number} mul - scalar multiplier
+* @returns {void}
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.2, 2.3, 3.4, 4.5, 1.1, 2.2, 3.3, 4.4, 5.5, 2.1, 3.2, 4.3, 5.4, 0.0, 3.1, 4.2, 5.3, 0.0, 0.0 ] );
+*
+* scaleBanded( 2, 1, 5, 5, A, 5, 1, 0, true, 10.0 );
+* // A => [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.0, 23.0, 34.0, 45.0, 11.0, 22.0, 33.0, 44.0, 55.0, 21.0, 32.0, 43.0, 54.0, 0.0, 31.0, 42.0, 53.0, 0.0, 0.0 ]
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 0.0, 0.0, 0.0, 1.1, 2.1, 3.1, 0.0, 0.0, 1.2, 2.2, 3.2, 4.2, 0.0, 0.0, 2.3, 3.3, 4.3, 5.3, 0.0, 0.0, 3.4, 4.4, 5.4, 0.0, 0.0, 0.0, 4.5, 5.5, 0.0, 0.0 ] );
+*
+* scaleBanded( 2, 1, 5, 5, A, 1, 6, 0, false, 10.0 );
+* // A => [ 0.0, 0.0, 0.0, 11.0, 21.0, 31.0, 0.0, 0.0, 12.0, 22.0, 32.0, 42.0, 0.0, 0.0, 23.0, 33.0, 43.0, 53.0, 0.0, 0.0, 34.0, 44.0, 54.0, 0.0, 0.0, 0.0, 45.0, 55.0, 0.0, 0.0 ]
+*/
+function scaleBanded( KL, KU, M, N, A, strideA1, strideA2, offsetA, isrm, mul ) {
+ var start;
+ var end;
+ var idx;
+ var dsa;
+ var ia;
+ var i0;
+ var i1;
+ var k1;
+ var k2;
+ var k3;
+ var k4;
+ var j;
+
+ k1 = KL + KU;
+ k2 = KL;
+ k3 = ( 2*KL ) + KU;
+ k4 = KL + KU + M - 1;
+ if ( isrm ) {
+ idx = offsetA + ( ( KL+KU ) * strideA1 );
+ dsa = strideA2 - strideA1;
+ for ( i1 = 0; i1 < M; i1++ ) {
+ start = max( i1-KL, 0 );
+ end = min( i1+KU, N-1 );
+ for ( i0 = start; i0 <= end; i0++ ) {
+ A[ idx ] *= mul;
+ idx += dsa;
+ }
+ idx -= ( ( end-start+1 ) * dsa ) - strideA1;
+ if ( i1+1 > KL ) {
+ idx += dsa;
+ }
+ }
+ return;
+ }
+ ia = offsetA;
+ for ( i1 = 0; i1 < N; i1++ ) {
+ j = max( k1-i1, k2 );
+ idx = ia + ( j * strideA1 );
+ for ( i0 = j; i0 <= min( k3, k4-i1 ); i0++ ) {
+ A[ idx ] *= mul;
+ idx += strideA1;
+ }
+ ia += strideA2;
+ }
+}
+
+
+// MAIN //
+
+/**
+* Multiplies a double-precision floating-point M-by-N matrix `A` by a double-precision floating-point scalar.
+*
+* @private
+* @param {string} type - single-letter character indicating the type of matrix `A`
+* @param {NonNegativeInteger} KL - lower bandwidth of `A`
+* @param {NonNegativeInteger} KU - upper bandwidth of `A`
+* @param {number} gamma - the matrix `A` is multiplied by `β/γ`
+* @param {number} beta - the matrix `A` is multiplied by `β/γ`
+* @param {NonNegativeInteger} M - number of rows in matrix `A`
+* @param {NonNegativeInteger} N - number of columns in matrix `A`
+* @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`
+* @returns {Float64Array} scaled matrix `A`
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 ] ); // => [ [ 1.0, 2.0 ], [ 3.0, 4.0 ], [ 5.0, 6.0 ] ]
+*
+* dlascl( 'G', 0, 0, 1.0, 2.0, 3, 2, A, 2, 1, 0 );
+* // A => [ 2.0, 4.0, 6.0, 8.0, 10.0, 12.0 ]
+*/
+function dlascl( type, KL, KU, gamma, beta, M, N, A, strideA1, strideA2, offsetA ) { // eslint-disable-line max-params
+ var isrm;
+ var done;
+ var mul;
+ var a;
+ var b;
+ var c;
+ var d;
+
+ if ( N === 0 || M === 0 ) {
+ return A;
+ }
+ isrm = isRowMajor( [ strideA1, strideA2 ] );
+
+ b = gamma; // denominator
+ a = beta; // numerator
+ done = false;
+ while ( !done ) {
+ c = b * TINY;
+
+ // Check whether `b` is infinite...
+ if ( c === b ) {
+ // Multiply by a correctly signed zero for finite `a` or NaN if `a` is infinite...
+ mul = a / b;
+ d = a;
+ done = true;
+ } else {
+ d = a / HUGE;
+
+ // Check whether `a` is either zero or infinite...
+ if ( d === a ) {
+ // In this case, `a` itself is the correct multiplication factor...
+ mul = a;
+ b = 1.0;
+ done = true;
+ } else if ( abs( c ) > abs( a ) && a !== 0.0 ) {
+ mul = TINY;
+ b = c;
+ done = false;
+ } else if ( abs( d ) > abs( b ) ) {
+ mul = HUGE;
+ a = d;
+ done = false;
+ } else {
+ mul = a / b;
+ done = true;
+ if ( mul === 1.0 ) {
+ return A;
+ }
+ }
+ }
+ if ( type === 'G' ) {
+ scaleGeneral( M, N, A, strideA1, strideA2, offsetA, mul );
+ } else if ( type === 'U' ) {
+ scaleUpper( M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ } else if ( type === 'L' ) {
+ scaleLower( M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ } else if ( type === 'H' ) {
+ scaleUpperHessenberg( M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ } else if ( type === 'B' ) {
+ scaleSymmetricBandedLower( KL, M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ } else if ( type === 'Q' ) {
+ scaleSymmetricBandedUpper( KU, M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ } else { // type === 'Z'
+ scaleBanded( KL, KU, M, N, A, strideA1, strideA2, offsetA, isrm, mul );
+ }
+ }
+ return A;
+}
+
+
+// EXPORTS //
+
+module.exports = dlascl;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq1.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq1.js
new file mode 100644
index 000000000000..584045859128
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq1.js
@@ -0,0 +1,74 @@
+/**
+* @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 base = require( './base.js' );
+
+
+// MAIN //
+
+/**
+* Computes the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E`.
+*
+* ## Notes
+*
+* - `D` should have `N` indexed elements. On entry, `D` contains the diagonal elements of the bi-diagonal matrix whose SVD is desired. On normal exit, `D` contains the singular values in decreasing order.
+*
+* - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements of the bi-diagonal matrix whose SVD is desired. On exit, E is overwritten.
+*
+* - The function returns a status code:
+*
+* - `= 0`: successful exit.
+*
+* - `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value.
+*
+* - `> 0`: the algorithm failed:
+* - `= 1`, a split was marked by a positive value in `E`.
+* - `= 2`, current block of `Z` not diagonalized after `100*N` iterations (in inner while loop). On exit, `D` and `E` represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into `DLASQ1`
+* - `= 3`, termination criterion of outer while loop not met (program created more than `N` unreduced blocks)
+*
+* @param {integer} N - number of rows/columns in the matrix
+* @param {Float64Array} D - the array with diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {Float64Array} E - the array with off-diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {Float64Array} WORK - workspace array (length >= 4*N)
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 100, 50, 25 ] );
+* var E = new Float64Array( [ 90, 40, 0 ] );
+* var WORK = new Float64Array( 16 );
+*
+* var info = dlasq1( 3, D, E, WORK );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90, 40, 0 ]
+* // WORK => [ ~1.939e292, ~3.137e291, ~2.554e290, ~1.575e286, ~3.137e291, ~2.881e291, ~2.278e292, ~2.278e292, 7.0, ~2.667, 0.0, ~2.881e291, 0, 0, 0, 0 ]
+* // info => 0
+*/
+function dlasq1( N, D, E, WORK ) { // eslint-disable-line stdlib/jsdoc-doctest-decimal-point
+ return base( N, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq1;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq2.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq2.js
new file mode 100644
index 000000000000..565a98c7330c
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq2.js
@@ -0,0 +1,527 @@
+/**
+* @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, max-statements, max-lines-per-function */
+
+// MODULES //
+
+var Float64Array = require( '@stdlib/array/float64' );
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+var sqrt = require( '@stdlib/math/base/special/sqrt' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var dlasq3 = require( './dlasq3.js' );
+var dlasrt = require( './dlasrt.js' );
+
+
+// VARIABLES //
+
+var CBIAS = 1.50;
+var EPS = dlamch( 'P' );
+var TOL = EPS * 100;
+var TOL2 = pow( TOL, 2 );
+var SAFMIN = dlamch( 'S' );
+var IEEE = true;
+
+
+// MAIN //
+
+/**
+* Computes all the eigenvalues of the symmetric positive definite tri-diagonal matrix associated with the QD Array `Z` to high relative accuracy.
+*
+* @private
+* @param {integer} N - number of rows/columns in `Z`
+* @param {Float64Array} Z - qd array
+* @param {integer} strideZ - stride length for `Z`
+* @param {NonNegativeInteger} offsetZ - starting index of `Z`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var Z = new Float64Array( [ 100, 4, 81, 3, 64, 2.5, 49, 2, 36, 1.5, 25, 1, 16, 0.5, 9, 0 ] );
+*
+* var out = dlasq2( 4, Z, 1, 0 );
+* // Z => [ ~115.713, ~83.153, ~62.17, ~42.464, ~83.153, ~20.987, ~0.0, ~0.0, 303.5, ~303.5, 15, 4.0625, 0.0, ~0.0, ~0.172, ~0.096 ]
+* // out => 0
+*/
+function dlasq2( N, Z, strideZ, offsetZ ) { // eslint-disable-line stdlib/jsdoc-doctest-decimal-point
+ var deemin;
+ var oldemn;
+ var iwhila;
+ var iwhilb;
+ var desig;
+ var dmin1;
+ var dmin2;
+ var sigma;
+ var trace;
+ var tempq;
+ var tempe;
+ var ttype;
+ var nfail;
+ var emax;
+ var emin;
+ var iter;
+ var nbig;
+ var ndiv;
+ var ipn4;
+ var kmin;
+ var dmin;
+ var qmax;
+ var qmin;
+ var temp;
+ var zmax;
+ var splt;
+ var dee;
+ var tau;
+ var dn1;
+ var dn2;
+ var out;
+ var pp;
+ var i0;
+ var i1;
+ var i4;
+ var n0;
+ var n1;
+ var dn;
+ var d;
+ var e;
+ var g;
+ var k;
+ var s;
+ var t;
+
+ out = new Float64Array( 18 );
+
+ if ( N === 0 ) {
+ return 0;
+ }
+
+ if ( N === 1 ) {
+ // 1-by-1 case
+ if ( Z[ offsetZ ] < 0 ) {
+ return -201;
+ }
+ return 0;
+ }
+
+ if ( N === 2 ) {
+ // 2-by-2 case
+ if ( Z[ offsetZ ] < 0 ) {
+ return -201;
+ }
+
+ if ( Z[ offsetZ + strideZ ] < 0 ) {
+ return -202;
+ }
+
+ if ( Z[ offsetZ + ( 2 * strideZ ) ] < 0 ) {
+ return -203;
+ }
+
+ if ( Z[ offsetZ + ( 2 * strideZ ) ] > Z[ offsetZ ] ) {
+ d = Z[ offsetZ + ( 2 * strideZ ) ];
+ Z[ offsetZ + ( 2 * strideZ ) ] = Z[ offsetZ ];
+ Z[ offsetZ ] = d;
+ }
+
+ Z[ offsetZ + ( 4 * strideZ ) ] = Z[ offsetZ ] + Z[ offsetZ + strideZ ] + Z[ offsetZ + ( 2 * strideZ ) ];
+
+ if ( Z[ offsetZ + strideZ ] > Z[ offsetZ + ( 2 * strideZ ) ] * TOL2 ) {
+ t = 0.5 * ( ( Z[ offsetZ ] - Z[ offsetZ + ( 2 * strideZ ) ] ) + Z[ offsetZ + strideZ ] );
+ s = Z[ offsetZ + ( 2 * strideZ ) ] * ( Z[ offsetZ + strideZ ] / t );
+ if ( s <= t ) {
+ s = Z[ offsetZ + ( 2 * strideZ ) ] * ( Z[ offsetZ + strideZ ] / ( t * ( 1 + sqrt( 1 + ( s / t ) ) ) ) );
+ } else {
+ s = Z[ offsetZ + ( 2 * strideZ ) ] * ( Z[ offsetZ + strideZ ] / ( t + ( sqrt( t ) * sqrt( t + s ) ) ) );
+ }
+ t = Z[ offsetZ ] + ( s + Z[ offsetZ + strideZ ] );
+ Z[ offsetZ + ( 2 * strideZ ) ] = Z[ offsetZ + ( 2 * strideZ ) ] * ( Z[ offsetZ ] / t );
+ Z[ offsetZ ] = t;
+ }
+ Z[ offsetZ + strideZ ] = Z[ offsetZ + ( 2 * strideZ ) ];
+ Z[ offsetZ + ( 5 * strideZ ) ] = Z[ offsetZ + strideZ ] + Z[ offsetZ ];
+ return 0;
+ }
+
+ // Check for negative data and compute sums of q's and e's
+ Z[ offsetZ + ( ( ( 2 * N ) - 1 ) * strideZ ) ] = 0;
+ emin = Z[ offsetZ + strideZ ];
+ qmax = 0;
+ zmax = 0;
+ d = 0;
+ e = 0;
+
+ for ( k = 0; k < 2 * ( N - 1 ); k += 2 ) {
+ if ( Z[ offsetZ + ( k * strideZ ) ] < 0 ) {
+ return -( 200 + k );
+ }
+ if ( Z[ offsetZ + ( k * strideZ ) + strideZ ] < 0 ) {
+ return -( 200 + k + 1 );
+ }
+ d += Z[ offsetZ + ( k * strideZ ) ];
+ e += Z[ offsetZ + ( k * strideZ ) + strideZ ];
+ qmax = max( qmax, Z[ offsetZ + ( k * strideZ ) ] );
+ emin = min( emin, Z[ offsetZ + ( k * strideZ ) + strideZ ] );
+ zmax = max( qmax, max( zmax, Z[ offsetZ + ( k * strideZ ) + strideZ ] ) );
+ }
+ if ( Z[ offsetZ + ( ( ( 2 * N ) - 2 ) * strideZ ) ] < 0 ) {
+ return -( 200 + ( 2 * N ) - 1 );
+ }
+ d += Z[ offsetZ + ( ( ( 2 * N ) - 2 ) * strideZ ) ];
+ qmax = max( qmax, Z[ offsetZ + ( ( ( 2 * N ) - 2 ) * strideZ ) ] );
+ zmax = max( qmax, zmax );
+
+ // Check for diagonality
+ if ( e === 0 ) {
+ for ( k = 1; k < N; k++ ) {
+ Z[ offsetZ + ( k * strideZ ) ] = Z[ offsetZ + ( ( ( 2 * k ) + 1 ) * strideZ ) ];
+ }
+ dlasrt( 'decreasing', N, Z, strideZ, offsetZ );
+ Z[ offsetZ + ( ( ( 2 * N ) - 2 ) * strideZ ) ] = d;
+ return 0;
+ }
+
+ trace = d + e;
+
+ // Check for zero data
+ if ( trace === 0 ) {
+ Z[ offsetZ + ( ( ( 2 * N ) - 2 ) * strideZ ) ] = 0;
+ return 0;
+ }
+
+ // Check whether the machine is IEEE conformable (In JS, always true)
+ IEEE = true;
+
+ // Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...)
+ for ( k = ( 2 * N ) - 1; k >= 1; k -= 2 ) {
+ Z[ offsetZ + ( ( ( 2 * k ) + 1 ) * strideZ ) ] = 0;
+ Z[ offsetZ + ( ( 2 * k ) * strideZ ) ] = Z[ offsetZ + ( k * strideZ ) ];
+ Z[ offsetZ + ( ( ( 2 * k ) - 1 ) * strideZ ) ] = 0;
+ Z[ offsetZ + ( ( ( 2 * k ) - 2 ) * strideZ ) ] = Z[ offsetZ + ( k * strideZ ) - strideZ ];
+ }
+
+ i0 = 0;
+ n0 = N-1;
+
+ // Reverse the qd-array, if warranted
+ if ( CBIAS * Z[ offsetZ + ( ( 4 * i0 ) * strideZ ) ] < Z[ offsetZ + ( ( 4 * n0 ) * strideZ ) ] ) {
+ ipn4 = ( 4 * ( i0 + n0 ) ) + 6;
+ for ( i4 = ( 4 * i0 ) + 3; i4 <= ( 2 * ( i0 + n0 ) ) + 1; i4 += 4 ) {
+ temp = Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 3 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 3 ) * strideZ ) ] = temp;
+ temp = Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 5 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 5 ) * strideZ ) ] = temp;
+ }
+ }
+
+ // Initial split checking via DQD and Li's test
+ pp = 0;
+
+ for ( k = 0; k < 2; k++ ) {
+ d = Z[ offsetZ + ( ( ( 4 * n0 ) + pp ) * strideZ ) ];
+ for ( i4 = ( 4 * n0 ) - 1 + pp; i4 >= ( 4 * i0 ) + pp + 3; i4 -= 4 ) {
+ if ( Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] <= TOL2 * d ) {
+ Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] = -0;
+ d = Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ];
+ } else {
+ d = Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ] * ( d / ( d + Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] ) );
+ }
+ }
+
+ // DQD maps Z to ZZ plus Li's test
+ emin = Z[ offsetZ + ( ( ( 4 * i0 ) + pp + 4 ) * strideZ ) ];
+ d = Z[ offsetZ + ( ( ( 4 * i0 ) + pp ) * strideZ ) ];
+ for ( i4 = ( 4 * i0 ) + pp + 3; i4 <= ( 4 * n0 ) + pp - 1; i4 += 4 ) {
+ Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] = d + Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ];
+ if ( Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] <= TOL2 * d ) {
+ Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] = -0;
+ Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] = d;
+ Z[ offsetZ + ( ( i4 - ( 2 * pp ) ) * strideZ ) ] = 0;
+ d = Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ];
+ } else if ( ( SAFMIN * Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] < Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] ) && ( SAFMIN * Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] < Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] ) ) {
+ temp = Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] / Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - ( 2 * pp ) ) * strideZ ) ] = Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] * temp;
+ d *= temp;
+ } else {
+ Z[ offsetZ + ( ( i4 - ( 2 * pp ) ) * strideZ ) ] = Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] * ( Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] / Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] );
+ d = Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] * ( d / Z[ offsetZ + ( ( i4 - ( 2 * pp ) - 2 ) * strideZ ) ] );
+ }
+ emin = min( emin, Z[ offsetZ + ( ( i4 - ( 2 * pp ) ) * strideZ ) ] );
+ }
+ Z[ offsetZ + ( ( ( 4 * n0 ) - pp + 1 ) * strideZ ) ] = d;
+
+ // Now find qmax
+ qmax = Z[ offsetZ + ( ( ( 4 * i0 ) - pp + 1 ) * strideZ ) ];
+ for ( i4 = ( 4 * i0 ) - pp + 5; i4 <= ( 4 * n0 ) - pp + 1; i4 += 4 ) {
+ qmax = max( qmax, Z[ offsetZ + ( i4 * strideZ ) ] );
+ }
+
+ // Prepare for the next iteration on K
+ pp = 1 - pp;
+ }
+
+ // Initialize variables to pass to DLASQ3
+ ttype = 0;
+ dmin1 = 0;
+ dmin2 = 0;
+ dn = 0;
+ dn1 = 0;
+ dn2 = 0;
+ g = 0;
+ tau = 0;
+
+ iter = 2;
+ nfail = 0;
+ ndiv = 2 * ( n0 - i0 );
+
+ for ( iwhila = 0; iwhila <= N; iwhila++ ) {
+ if ( n0 < 0 ) {
+ // Move q's to the front
+ for ( k = 1; k < N; k++ ) {
+ Z[ offsetZ + ( k * strideZ ) ] = Z[ offsetZ + ( ( 4 * k ) * strideZ ) ];
+ }
+
+ // Sort and compute sum of eigenvalues
+ dlasrt( 'decreasing', N, Z, strideZ, offsetZ );
+
+ e = 0;
+ for ( k = N - 1; k >= 0; k-- ) {
+ e += Z[ offsetZ + ( k * strideZ ) ];
+ }
+
+ // Store trace, sum(eigenvalues) and information on performance
+ Z[ offsetZ + ( ( 2 * N ) * strideZ ) ] = trace;
+ Z[ offsetZ + ( ( ( 2 * N ) + 1 ) * strideZ ) ] = e;
+ Z[ offsetZ + ( ( ( 2 * N ) + 2 ) * strideZ ) ] = iter;
+ Z[ offsetZ + ( ( ( 2 * N ) + 3 ) * strideZ ) ] = ndiv / ( N * N );
+ Z[ offsetZ + ( ( ( 2 * N ) + 4 ) * strideZ ) ] = ( 100 * nfail ) / iter;
+ return 0;
+ }
+
+ // E(N0) holds the value of SIGMA when submatrix in I0:N0 splits from the rest of the array, but is negated.
+ desig = 0;
+ if ( n0 === N - 1 ) {
+ sigma = 0;
+ } else {
+ sigma = -Z[ offsetZ + ( ( ( 4 * n0 ) + 2 ) * strideZ ) ];
+ }
+ if ( sigma < 0 ) {
+ return 1;
+ }
+
+ // Find last unreduced submatrix's top index I0, find QMAX and EMIN. Find Gershgorin-type bound if Q's much greater than E's.
+ emax = 0;
+ if ( n0 > i0 ) {
+ emin = abs( Z[ offsetZ + ( ( ( 4 * n0 ) - 2 ) * strideZ ) ] );
+ } else {
+ emin = 0;
+ }
+ qmin = Z[ offsetZ + ( ( 4 * n0 ) * strideZ ) ];
+ qmax = qmin;
+
+ for ( i4 = ( 4 * n0 ) + 3; i4 >= 7; i4 -= 4 ) {
+ if ( Z[ offsetZ + ( ( i4 - 5 ) * strideZ ) ] <= 0 ) {
+ break;
+ }
+ if ( qmin >= 4 * emax ) {
+ qmin = min( qmin, Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ] );
+ emax = max( emax, Z[ offsetZ + ( ( i4 - 5 ) * strideZ ) ] );
+ }
+ qmax = max( qmax, Z[ offsetZ + ( ( i4 - 7 ) * strideZ ) ] + Z[ offsetZ + ( ( i4 - 5 ) * strideZ ) ] );
+ emin = min( emin, Z[ offsetZ + ( ( i4 - 5 ) * strideZ ) ] );
+ }
+
+ // If the loop completed without break, set i4 = 4
+ if ( i4 < 7 ) {
+ i4 = 3;
+ }
+
+ i0 = floor( ( ( i4 + 1 ) / 4 ) - 1 );
+ pp = 0;
+
+ if ( n0 - i0 > 1 ) {
+ dee = Z[ offsetZ + ( ( 4 * i0 ) * strideZ ) ];
+ deemin = dee;
+ kmin = i0;
+ for ( i4 = ( 4 * i0 ) + 4; i4 <= 4 * n0; i4 += 4 ) {
+ dee = Z[ offsetZ + ( i4 * strideZ ) ] * ( dee / ( dee + Z[ offsetZ + ( ( i4 - 2 ) * strideZ ) ] ) );
+ if ( dee <= deemin ) {
+ deemin = dee;
+ kmin = floor( i4 / 4 );
+ }
+ }
+ if ( ( ( kmin - i0 ) * 2 < n0 - kmin ) && ( deemin <= 0.5 * Z[ offsetZ + ( ( 4 * n0 ) * strideZ ) ] ) ) {
+ ipn4 = ( 4 * ( i0 + n0 ) ) + 7;
+ pp = 2;
+ for ( i4 = ( 4 * i0 ) + 3; i4 <= ( 2 * ( i0 + n0 ) ) + 1; i4 += 4 ) {
+ temp = Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 3 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 3 ) * strideZ ) ] = temp;
+ temp = Z[ offsetZ + ( ( i4 - 2 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - 2 ) * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 2 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 2 ) * strideZ ) ] = temp;
+ temp = Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ];
+ Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 5 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 5 ) * strideZ ) ] = temp;
+ temp = Z[ offsetZ + ( i4 * strideZ ) ];
+ Z[ offsetZ + ( i4 * strideZ ) ] = Z[ offsetZ + ( ( ipn4 - i4 - 4 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ipn4 - i4 - 4 ) * strideZ ) ] = temp;
+ }
+ }
+ }
+
+ // Put -(initial shift) into DMIN
+ dmin = -max( 0, qmin - ( 2 * sqrt( qmin ) * sqrt( emax ) ) );
+
+ /*
+ Now I0:N0 is unreduced.
+ PP = 0 for ping,
+ PP = 1 for pong,
+ PP = 2 indicates that flipping was applied to the Z array and that the tests for deflation upon entry in DLASQ3 should not be performed.
+ */
+ nbig = 100 * ( n0 - i0 + 1 );
+ for ( iwhilb = 0; iwhilb < nbig; iwhilb++ ) {
+ if ( i0 > n0 ) {
+ break;
+ }
+
+ out[ 0 ] = n0;
+ out[ 1 ] = pp;
+ out[ 2 ] = dmin;
+ out[ 3 ] = sigma;
+ out[ 4 ] = desig;
+ out[ 5 ] = qmax;
+ out[ 6 ] = nfail;
+ out[ 7 ] = iter;
+ out[ 8 ] = ndiv;
+ out[ 9 ] = ttype;
+ out[ 10 ] = dmin1;
+ out[ 11 ] = dmin2;
+ out[ 12 ] = dn;
+ out[ 13 ] = dn1;
+ out[ 14 ] = dn2;
+ out[ 15 ] = g;
+ out[ 16 ] = tau;
+
+ // Call dlasq3
+ dlasq3( i0, Z, strideZ, offsetZ, IEEE, out, 1, 0 );
+
+ n0 = out[ 0 ];
+ pp = out[ 1 ];
+ dmin = out[ 2 ];
+ sigma = out[ 3 ];
+ desig = out[ 4 ];
+ qmax = out[ 5 ];
+ nfail = out[ 6 ];
+ iter = out[ 7 ];
+ ndiv = out[ 8 ];
+ ttype = out[ 9 ];
+ dmin1 = out[ 10 ];
+ dmin2 = out[ 11 ];
+ dn = out[ 12 ];
+ dn1 = out[ 13 ];
+ dn2 = out[ 14 ];
+ g = out[ 15 ];
+ tau = out[ 16 ];
+
+ pp = 1 - pp;
+
+ // When EMIN is very small check for splits
+ if ( pp === 0 && n0 - i0 >= 3 ) {
+ if ( ( Z[ offsetZ + ( ( ( 4 * n0 ) + 3 ) * strideZ ) ] <= TOL2 * qmax ) || ( Z[ offsetZ + ( ( ( 4 * n0 ) + 2 ) * strideZ ) ] <= TOL2 * sigma ) ) {
+ splt = i0 - 1;
+ qmax = Z[ offsetZ + ( ( 4 * i0 ) * strideZ ) ];
+ emin = Z[ offsetZ + ( ( ( 4 * i0 ) + 2 ) * strideZ ) ];
+ oldemn = Z[ offsetZ + ( ( ( 4 * i0 ) + 3 ) * strideZ ) ];
+ for ( i4 = ( 4 * i0 ) + 3; i4 <= ( 4 * n0 ) - 9; i4 += 4 ) {
+ if ( ( Z[ offsetZ + ( i4 * strideZ ) ] <= TOL2 * Z[ offsetZ + ( ( i4 - 3 ) * strideZ ) ] ) || ( Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] <= TOL2 * sigma) ) {
+ Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] = -sigma;
+ splt = floor( ( i4 + 1 ) / 4 ) - 1;
+ qmax = 0;
+ emin = Z[ offsetZ + ( ( i4 + 3 ) * strideZ ) ];
+ oldemn = Z[ offsetZ + ( ( i4 + 4 ) * strideZ ) ];
+ } else {
+ qmax = max( qmax, Z[ offsetZ + ( ( i4 + 1 ) * strideZ ) ] );
+ emin = min( emin, Z[ offsetZ + ( ( i4 - 1 ) * strideZ ) ] );
+ oldemn = min( oldemn, Z[ offsetZ + ( i4 * strideZ ) ] );
+ }
+ }
+ Z[ offsetZ + ( ( ( 4 * n0 ) + 2 ) * strideZ ) ] = emin;
+ Z[ offsetZ + ( ( ( 4 * n0 ) + 3 ) * strideZ ) ] = oldemn;
+ i0 = splt + 1;
+ }
+ }
+ }
+
+ // If inner loop exhausted without breaking, we have INFO = 2
+ if ( iwhilb >= nbig ) {
+ // Maximum number of iterations exceeded, restore the shift SIGMA and place the new d's and e's in a qd array. This might need to be done for several blocks
+
+ i1 = i0;
+ n1 = n0;
+
+ // Label 145 loop
+ while ( true ) {
+ tempq = Z[ offsetZ + ( ( 4 * i0 ) * strideZ ) ];
+ Z[ offsetZ + ( ( 4 * i0 ) * strideZ ) ] += sigma;
+ for ( k = i0 + 1; k <= n0; k++ ) {
+ tempe = Z[ offsetZ + ( ( ( 4 * k ) - 2 ) * strideZ ) ];
+ Z[ offsetZ + ( ( ( 4 * k ) - 2 ) * strideZ ) ] *= tempq / Z[ offsetZ + ( ( ( 4 * k ) - 4 ) * strideZ ) ];
+ tempq = Z[ offsetZ + ( ( 4 * k ) * strideZ ) ];
+ Z[ offsetZ + ( ( 4 * k ) * strideZ ) ] += sigma + tempe - Z[ offsetZ + ( ( ( 4 * k ) - 2 ) * strideZ ) ];
+ }
+
+ // Prepare to do this on the previous block if there is one
+ if ( i1 > 0 ) {
+ n1 = i1 - 1;
+ while ( i1 >= 1 && Z[ offsetZ + ( ( ( 4 * i1 ) - 2 ) * strideZ ) ] >= 0 ) {
+ i1 -= 1;
+ }
+ sigma = -Z[ offsetZ + ( ( ( 4 * n1 ) + 2 ) * strideZ ) ];
+ } else {
+ break;
+ }
+ }
+
+ for ( k = 0; k < N; k++ ) {
+ Z[ offsetZ + ( ( 2 * k ) * strideZ ) ] = Z[ offsetZ + ( ( 4 * k ) * strideZ ) ];
+ if ( k < n0 ) {
+ Z[ offsetZ + ( ( ( 2 * k ) + 1 ) * strideZ ) ] = Z[ offsetZ + ( ( ( 4 * k ) + 2 ) * strideZ ) ];
+ } else {
+ Z[ offsetZ + ( ( ( 2 * k ) + 1 ) * strideZ ) ] = 0;
+ }
+ }
+ return 2;
+ }
+ }
+
+ return 3;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq2;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq3.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq3.js
new file mode 100644
index 000000000000..5b4afe58b8e1
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq3.js
@@ -0,0 +1,379 @@
+/**
+* @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, max-statements, max-lines-per-function */
+
+// MODULES //
+
+var Float64Array = require( '@stdlib/array/float64' );
+var dlamch = require( '@stdlib/lapack/base/dlamch' );
+var isnan = require( '@stdlib/assert/is-nan' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+var sqrt = require( '@stdlib/math/base/special/sqrt' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var dlasq4 = require( './dlasq4.js' );
+var dlasq5 = require( './dlasq5.js' );
+var dlasq6 = require( './dlasq6.js' );
+
+
+// VARIABLES //
+
+var CBIAS = 1.50;
+var EPS = dlamch( 'P' );
+var TOL = EPS * 100;
+var TOL2 = pow( TOL, 2 );
+
+
+// MAIN //
+
+/**
+* Checks for deflation, computes a shift (`TAU`) and calls DQDS. In case of failure it changes shifts, and tries again until output is positive.
+*
+* @private
+* @param {integer} I0 - first index
+* @param {Float64Array} Z - qd array
+* @param {integer} strideZ - stride length for `Z`
+* @param {NonNegativeInteger} offsetZ - starting index of `Z`
+* @param {boolean} IEEE - IEEE arithmetic flag
+* @param {Float64Array} out - output array containing `N0`, `PP`, `DMIN`, `SIGMA`, `DESIG`, `QMAX`, `NFAIL`, `ITER`, `NDIV`, `TTYPE`, `DMIN1`, `DMIN2`, `DN`, `DN1`, `DN2`, `G` and `TAU`
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index of `out`
+* @returns {Float64Array} output array
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var out = new Float64Array( 17 );
+* out[ 0 ] = 3;
+* var Z = new Float64Array( [ 100, 4, 81, 3, 64, 2.5, 49, 2, 36, 1.5, 25, 1, 16, 0.5, 9, 0 ] );
+*
+* dlasq3( 0, Z, 1, 0, true, out, 1, 0 );
+* // out => [ 3.0, 0.0, ~6.022, 0.0, 0.0, 0.0, 0.0, 1.0, 5.0, -1.0, ~15.089, ~35.359, ~6.022, ~15.089, ~35.359, 0.0, 0.0 ]
+* // Z => [ 100.0, 181.0, 81.0, ~28.641, 64.0, ~84.359, 49.0, ~20.911, ~36.0, ~40.089, 25.0, ~9.978, 16.0, ~6.022, 9.0, ~28.641 ]
+*/
+function dlasq3( I0, Z, strideZ, offsetZ, IEEE, out, strideOut, offsetOut ) {
+ var goto80;
+ var desig;
+ var dmin1;
+ var dmin2;
+ var nfail;
+ var sigma;
+ var ttype;
+ var dmin;
+ var idx1;
+ var idx2;
+ var idx3;
+ var idx4;
+ var ipn4;
+ var iter;
+ var n0in;
+ var ndiv;
+ var out1;
+ var temp;
+ var qmax;
+ var dn1;
+ var dn2;
+ var idx;
+ var tau;
+ var dn;
+ var j4;
+ var N0;
+ var nn;
+ var PP;
+ var g;
+ var s;
+ var t;
+
+ // Read input values from output array
+ idx = offsetOut;
+ N0 = out[ idx ];
+ idx += strideOut;
+ PP = out[ idx ];
+ idx += strideOut;
+ dmin = out[ idx ];
+ idx += strideOut;
+ sigma = out[ idx ];
+ idx += strideOut;
+ desig = out[ idx ];
+ idx += strideOut;
+ qmax = out[ idx ];
+ idx += strideOut;
+ nfail = out[ idx ];
+ idx += strideOut;
+ iter = out[ idx ];
+ idx += strideOut;
+ ndiv = out[ idx ];
+ idx += strideOut;
+ ttype = out[ idx ];
+ idx += strideOut;
+ dmin1 = out[ idx ];
+ idx += strideOut;
+ dmin2 = out[ idx ];
+ idx += strideOut;
+ dn = out[ idx ];
+ idx += strideOut;
+ dn1 = out[ idx ];
+ idx += strideOut;
+ dn2 = out[ idx ];
+ idx += strideOut;
+ g = out[ idx ];
+ idx += strideOut;
+ tau = out[ idx ];
+
+ n0in = N0;
+ idx3 = offsetZ + ( strideZ * ( 4 * N0 ) );
+ idx4 = offsetZ + ( strideZ * ( 4 * I0 ) );
+
+ // Check for deflation
+ while ( true ) {
+ idx3 = offsetZ + ( strideZ * ( 4 * N0 ) );
+ if ( N0 < I0 ) {
+ out[ offsetOut ] = N0;
+ out[ offsetOut + strideOut ] = PP;
+ return out;
+ }
+ if ( N0 === I0 ) {
+ // Deflate 1 eigenvalue
+ Z[ idx3 ] = Z[ idx3 + ( strideZ * PP ) ] + sigma;
+ N0 -= 1;
+ continue;
+ }
+ nn = offsetZ + ( strideZ * ( ( 4 * N0 ) + PP + 3 ) );
+ if ( N0 !== ( I0 + 1 ) ) {
+ if ( Z[ nn - ( 5 * strideZ ) ] > TOL2 * ( sigma + Z[ nn - ( 3 * strideZ ) ] ) && Z[ nn - ( strideZ * ( ( 2 * PP ) + 4 ) ) ] > TOL2 * Z[ nn - ( 7 * strideZ ) ] ) {
+ // Check whether E(N0-1) is negligible, 1 eigenvalue
+ if ( Z[ nn - ( 9 * strideZ ) ] > TOL2 * sigma &&
+ Z[ nn - ( strideZ * ( ( 2 * PP ) + 8 ) ) ] > TOL2 * Z[ nn - ( 11 * strideZ ) ] ) {
+ break; // No deflation possible
+ }
+ } else {
+ Z[ idx3 ] = Z[ idx3 + ( strideZ * PP ) ] + sigma;
+ N0 -= 1;
+ continue;
+ }
+ }
+
+ // Check whether E(N0-2) is negligible, 2 eigenvalues
+ if ( Z[ nn - ( 3 * strideZ ) ] > Z[ nn - ( 7 * strideZ ) ] ) {
+ s = Z[ nn - ( 3 * strideZ ) ];
+ Z[ nn - ( 3 * strideZ ) ] = Z[ nn - ( 7 * strideZ ) ];
+ Z[ nn - ( 7 * strideZ ) ] = s;
+ }
+ t = 0.5 * ( ( Z[ nn - ( 7 * strideZ ) ] - Z[ nn - ( 3 * strideZ ) ] ) + Z[ nn - ( 5 * strideZ ) ] );
+ if ( Z[ nn - ( 5 * strideZ ) ] > Z[ nn - ( 3 * strideZ ) ] * TOL2 && t !== 0.0 ) {
+ s = Z[ nn - ( 3 * strideZ ) ] * ( Z[ nn - ( 5 * strideZ ) ] / t );
+ if ( s <= t ) {
+ s = Z[ nn - ( 3 * strideZ ) ] * ( Z[ nn - ( 5 * strideZ ) ] / ( t * ( 1.0 + sqrt( 1.0 + ( s / t ) ) ) ) );
+ } else {
+ s = Z[ nn - ( 3 * strideZ ) ] * ( Z[ nn - ( 5 * strideZ ) ] / ( t + ( sqrt( t ) * sqrt( t + s ) ) ) );
+ }
+ t = Z[ nn - ( 7 * strideZ ) ] + ( s + Z[ nn - ( 5 * strideZ ) ] );
+ Z[ nn - ( 3 * strideZ ) ] = Z[ nn - ( 3 * strideZ ) ] * ( Z[ nn - ( 7 * strideZ ) ] / t );
+ Z[ nn - ( 7 * strideZ ) ] = t;
+ }
+ Z[ idx3 - ( strideZ * 4 ) ] = Z[ nn - ( 7 * strideZ ) ] + sigma;
+ Z[ idx3 ] = Z[ nn - ( 3 * strideZ ) ] + sigma;
+ N0 -= 2;
+ }
+
+ // 50: No deflation, need shift
+ if ( PP === 2 ) {
+ PP = 0;
+ }
+
+ // Reverse the qd-array, if warranted
+ if ( dmin <= 0.0 || N0 < n0in ) {
+ if ( CBIAS * Z[ idx4 + ( strideZ * PP ) ] < Z[ idx3 + ( strideZ * PP ) ] ) {
+ ipn4 = ( 4 * ( I0 + N0 ) ) + 7;
+ for ( j4 = ( 4 * I0 ) + 3; j4 <= ( 2 * ( I0 + N0 ) ) + 1; j4 += 4 ) {
+ idx1 = offsetZ + ( strideZ * ( j4 - 3 ) );
+ idx2 = offsetZ + ( strideZ * ( ipn4 - j4 - 4 ) );
+
+ temp = Z[ idx1 ];
+ Z[ idx1 ] = Z[ idx2 ];
+ Z[ idx2 ] = temp;
+
+ idx1 += strideZ;
+ idx2 += strideZ;
+
+ temp = Z[ idx1 ];
+ Z[ idx1 ] = Z[ idx2 ];
+ Z[ idx2 ] = temp;
+
+ idx1 += strideZ;
+ idx2 -= 3*strideZ;
+
+ temp = Z[ idx1 ];
+ Z[ idx1 ] = Z[ idx2 ];
+ Z[ idx2 ] = temp;
+
+ idx1 += strideZ;
+ idx2 += strideZ;
+
+ temp = Z[ idx1 ];
+ Z[ idx1 ] = Z[ idx2 ];
+ Z[ idx2 ] = temp;
+ }
+ if ( ( N0 - I0 ) <= 4 ) {
+ Z[ idx3 + ( strideZ * ( PP + 2 ) ) ] = Z[ idx4 + ( strideZ * ( PP + 2 ) ) ];
+ Z[ idx3 + ( strideZ * ( 3 - PP ) ) ] = Z[ idx4 + ( strideZ * ( 3 - PP ) ) ];
+ }
+ dmin2 = min( dmin2, Z[ idx3 + ( strideZ * ( PP + 2 ) ) ] );
+ Z[ idx3 + ( strideZ * ( PP + 2 ) ) ] = min( Z[ idx3 + ( strideZ * ( PP + 2 ) ) ], min( Z[ idx4 + ( strideZ * ( PP + 2 ) ) ], Z[ idx4 + ( strideZ * ( PP + 6 ) ) ] ) );
+ Z[ idx3 + ( strideZ * ( 3 - PP ) ) ] = min( Z[ idx3 + ( strideZ * ( 3 - PP ) ) ], min( Z[ idx4 + ( strideZ * ( 3 - PP ) ) ], Z[ idx4 + ( strideZ * ( 7 - PP ) ) ] ) );
+ qmax = max( qmax, Z[ idx4 + ( strideZ * ( PP ) ) ], Z[ idx4 + ( strideZ * ( PP + 4 ) ) ] );
+ dmin = -0.0;
+ }
+ }
+
+ // Choose a shift
+ out1 = new Float64Array( 3 );
+ out1[ 1 ] = ttype;
+ out1[ 2 ] = g;
+ dlasq4( I0, N0, Z, strideZ, offsetZ, PP, n0in, dmin, dmin1, dmin2, dn, dn1, dn2, out1, 1, 0 );
+ tau = out1[ 0 ];
+ ttype = out1[ 1 ];
+ g = out1[ 2 ];
+
+ // Call DQDS until DMIN > 0
+ out1 = new Float64Array( 6 );
+ while ( true ) {
+ dlasq5( I0, N0, Z, strideZ, offsetZ, PP, tau, sigma, IEEE, EPS, out1, 1, 0 );
+ dmin = out1[ 0 ];
+ dmin1 = out1[ 1 ];
+ dmin2 = out1[ 2 ];
+ dn = out1[ 3 ];
+ dn1 = out1[ 4 ];
+ dn2 = out1[ 5 ];
+
+ ndiv += ( N0 - I0 + 2 );
+ iter += 1;
+
+ // Check status
+ if ( dmin >= 0.0 && dmin1 >= 0.0 ) {
+ // Success
+ goto80 = false;
+ break;
+ } else if ( dmin < 0.0 && dmin1 > 0.0 && Z[ idx3 - ( strideZ * ( PP + 1 ) ) ] < TOL * ( sigma + dn1 ) && abs( dn ) < TOL * sigma ) {
+ // Convergence hidden by negative DN
+ Z[ idx3 + ( strideZ * ( 1 - PP ) ) ] = 0.0;
+ dmin = 0.0;
+ goto80 = false;
+ break;
+ } else if ( dmin < 0.0 ) {
+ // TAU too big. Select new TAU and try again.
+ nfail += 1;
+ if ( ttype < -22 ) {
+ // Failed twice. Play it safe.
+ tau = 0.0;
+ } else if ( dmin1 > 0.0 ) {
+ // Late failure. Gives excellent shift.
+ tau = ( tau + dmin ) * ( 1.0 - ( 2.0 * EPS ) );
+ ttype -= 11;
+ } else {
+ // Early failure. Divide by 4.
+ tau *= 0.25;
+ ttype -= 12;
+ }
+ continue;
+ } else if ( isnan( dmin ) ) {
+ // NaN
+ if ( tau === 0.0 ) {
+ goto80 = true;
+ break;
+ }
+ tau = 0.0;
+ continue;
+ }
+ // Possible underflow. Play it safe.
+ goto80 = true;
+ break;
+ }
+
+ // Risk of underflow
+ if ( goto80 ) {
+ out1 = new Float64Array( 6 );
+ dlasq6( I0, N0, Z, strideZ, offsetZ, PP, out1, 1, 0 );
+ dmin = out1[ 0 ];
+ dmin1 = out1[ 1 ];
+ dmin2 = out1[ 2 ];
+ dn = out1[ 3 ];
+ dn1 = out1[ 4 ];
+ dn2 = out1[ 5 ];
+ ndiv += ( N0 - I0 + 2 );
+ iter += 1;
+ tau = 0.0;
+ }
+
+ if ( tau < sigma ) {
+ desig += tau;
+ t = sigma + desig;
+ desig -= ( t - sigma );
+ } else {
+ t = sigma + tau;
+ desig = sigma - ( t - tau ) + desig;
+ }
+ sigma = t;
+
+ // Store output values
+ idx = offsetOut;
+ out[ idx ] = N0;
+ idx += strideOut;
+ out[ idx ] = PP;
+ idx += strideOut;
+ out[ idx ] = dmin;
+ idx += strideOut;
+ out[ idx ] = sigma;
+ idx += strideOut;
+ out[ idx ] = desig;
+ idx += strideOut;
+ out[ idx ] = qmax;
+ idx += strideOut;
+ out[ idx ] = nfail;
+ idx += strideOut;
+ out[ idx ] = iter;
+ idx += strideOut;
+ out[ idx ] = ndiv;
+ idx += strideOut;
+ out[ idx ] = ttype;
+ idx += strideOut;
+ out[ idx ] = dmin1;
+ idx += strideOut;
+ out[ idx ] = dmin2;
+ idx += strideOut;
+ out[ idx ] = dn;
+ idx += strideOut;
+ out[ idx ] = dn1;
+ idx += strideOut;
+ out[ idx ] = dn2;
+ idx += strideOut;
+ out[ idx ] = g;
+ idx += strideOut;
+ out[ idx ] = tau;
+
+ return out;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq3;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq4.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq4.js
new file mode 100644
index 000000000000..b1fd9010c09d
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq4.js
@@ -0,0 +1,413 @@
+/**
+* @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 max = require( '@stdlib/math/base/special/max' );
+var min = require( '@stdlib/math/base/special/min' );
+var sqrt = require( '@stdlib/math/base/special/sqrt' );
+
+
+// VARIABLES //
+
+var CNST1 = 0.563;
+var CNST2 = 1.010;
+var CNST3 = 1.050;
+
+
+// MAIN //
+
+/**
+* Computes an approximation to the smallest eigenvalue using values of d from the previous transform.
+*
+* ## Notes
+*
+* - `Z` is a 1-D array of length >= `4*N0` storing interleaved q/e values.
+* - `PP` is `0` for ping, `1` for pong.
+* - `TAU` is approximation to the smallest eigenvalue and is used as a shift to accelerate convergence of the DQDS iteration.
+* - `TTYPE` is an integer flag describing how TAU was computed.
+* - `G` is a state variable that is preserved across successive calls and is used to regulate the magnitude of fallback shifts.
+*
+* @private
+* @param {integer} I0 - first index
+* @param {integer} N0 - last index
+* @param {Float64Array} Z - the QD array
+* @param {integer} strideZ - stride length for `Z`
+* @param {NonNegativeInteger} offsetZ - starting index for `Z`
+* @param {integer} PP - ping-pong flag (0 or 1)
+* @param {integer} N0IN - value of `N0` at the start of `EIGTEST`
+* @param {number} DMIN - minimum value of `d`
+* @param {number} DMIN1 - minimum value of `d`, excluding `D(N0)`
+* @param {number} DMIN2 - minimum value of `d`, excluding `D(N0)` and `D(N0-1)`
+* @param {number} DN - `d(N)`
+* @param {number} DN1 - `d(N-1)`
+* @param {number} DN2 - `d(N-2)`
+* @param {Float64Array} out - output array containing `tau`, `ttype`, and `G`
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index for `out`
+* @returns {Float64Array} output array
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var out = new Float64Array( [ 0, 0, 0.25 ] );
+* var Z = new Float64Array( [ 5, 0, 7, 0, 9, 0, 3, 0, 10, 0, 4.5, 0, 18, 0, 0, 0 ] );
+*
+* dlasq4( 0, 3, Z, 1, 0, 0, 3, 0.2, 0.15, 0.1, 0.8, 0.7, 0.6, out, 1, 0 );
+* // out => [ ~0.05, -6, 0.25 ]
+*/
+function dlasq4( I0, N0, Z, strideZ, offsetZ, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, DN2, out, strideOut, offsetOut ) { // eslint-disable-line max-len, max-params
+ var ttype;
+ var idx2;
+ var gap1;
+ var gap2;
+ var gam;
+ var idx;
+ var inc;
+ var ix1;
+ var tau;
+ var a2;
+ var b1;
+ var b2;
+ var i4;
+ var nn;
+ var np;
+ var G;
+ var s;
+
+ idx2 = offsetOut;
+ tau = out[ idx2 ];
+ idx2 += strideOut;
+ ttype = out[ idx2 ];
+ idx2 += strideOut;
+ G = out[ idx2 ];
+
+ // A negative DMIN forces the shift to take that absolute value TTYPE records the type of shift.
+ if ( DMIN <= 0 ) {
+ tau = -DMIN;
+ ttype = -1;
+
+ idx2 = offsetOut;
+ out[ idx2 ] = tau;
+ idx2 += strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+
+ nn = offsetZ + ( strideZ * ( ( 4 * N0 ) + PP + 3 ) );
+
+ inc = 4 * strideZ;
+
+ if ( N0IN === N0 ) {
+ // No eigenvalues deflated.
+ if ( DMIN === DN || DMIN === DN1 ) {
+ b1 = sqrt( Z[ nn - ( 3 * strideZ ) ] ) * sqrt( Z[ nn - ( 5 * strideZ ) ] );
+ b2 = sqrt( Z[ nn - ( 7 * strideZ ) ] ) * sqrt( Z[ nn - ( 9 * strideZ ) ] );
+ a2 = Z[ nn - ( 7 * strideZ ) ] + Z[ nn - ( 5 * strideZ ) ];
+
+ // Cases 2 and 3.
+ if ( DMIN === DN && DMIN1 === DN1 ) {
+ gap2 = DMIN2 - a2 - ( DMIN2 * 0.25 );
+
+ if ( gap2 > 0 && gap2 > b2 ) {
+ gap1 = a2 - DN - ( ( b2 / gap2 ) * b2 );
+ } else {
+ gap1 = a2 - DN - ( b1 + b2 );
+ }
+ if ( gap1 > 0 && gap1 > b1 ) {
+ s = max( DN - ( ( b1 / gap1 ) * b1 ), 0.5 * DMIN );
+ ttype = -2;
+ } else {
+ s = 0;
+ if ( DN > b1 ) {
+ s = DN - b1;
+ }
+ if ( a2 > ( b1 + b2 ) ) {
+ s = min( s, a2 - ( b1 + b2 ) );
+ }
+ s = max( s, 0.333 * DMIN );
+ ttype = -3;
+ }
+ } else {
+ // Case 4.
+ ttype = -4;
+ s = 0.25 * DMIN;
+
+ if ( DMIN === DN ) {
+ gam = DN;
+ a2 = 0;
+ if ( Z[ nn - ( 5 * strideZ ) ] > Z[ nn - ( 7 * strideZ ) ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b2 = Z[ nn - ( 5 * strideZ ) ] / Z[ nn - ( 7 * strideZ ) ];
+ np = nn - ( 9 * strideZ );
+ } else {
+ np = nn - ( 2 * PP * strideZ );
+ gam = DN1;
+ if ( Z[ np - ( 4 * strideZ ) ] > Z[ np - ( 2 * strideZ ) ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ a2 = Z[ np - ( 4 * strideZ ) ] / Z[ np - ( 2 * strideZ ) ];
+ if ( Z[ nn - ( 9 * strideZ ) ] > Z[ nn - ( 11 * strideZ ) ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b2 = Z[ nn - ( 9 * strideZ ) ] / Z[ nn - ( 11 * strideZ ) ];
+ np = nn - ( 13 * strideZ );
+ }
+
+ // Approximate contribution to norm squared from I < NN-1.
+ a2 += b2;
+ i4 = np;
+ ix1 = i4 - ( 2 * strideZ );
+ for ( idx = ( np - offsetZ ) / strideZ; idx >= ( 4 * I0 ) + PP + 2; idx -= 4 ) {
+ if ( b2 === 0 ) {
+ break;
+ }
+ b1 = b2;
+ if ( Z[ i4 ] > Z[ ix1 ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b2 *= Z[ i4 ] / Z[ ix1 ];
+ a2 += b2;
+ if ( ( 100 * max( b2, b1 ) ) < a2 || CNST1 < a2 ) {
+ break;
+ }
+ i4 -= inc;
+ ix1 -= inc;
+ }
+ a2 *= CNST3;
+
+ // Rayleigh quotient residual bound.
+ if ( a2 < CNST1 ) {
+ s = gam * ( 1 - sqrt( a2 ) ) / ( 1 + a2 );
+ }
+ }
+ } else if ( DMIN === DN2 ) {
+ // Case 5.
+ ttype = -5;
+ s = 0.25 * DMIN;
+
+ // Compute contribution to norm squared from I > NN-2.
+ np = nn - ( 2 * PP * strideZ );
+ b1 = Z[ np - ( 2 * strideZ ) ];
+ b2 = Z[ np - ( 6 * strideZ ) ];
+ gam = DN2;
+ if ( Z[ np - ( 8 * strideZ ) ] > b2 || Z[ np - ( 4 * strideZ ) ] > b1 ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ a2 = ( Z[ np - ( 8 * strideZ ) ] / b2 ) * ( 1 + ( Z[ np - ( 4 * strideZ ) ] / b1 ) );
+
+ // Approximate contribution to norm squared from I < NN-2.
+ if ( ( N0 - I0 ) > 2 ) {
+ b2 = Z[ nn - ( 13 * strideZ ) ] / Z[ nn - ( 15 * strideZ ) ];
+ a2 += b2;
+
+ i4 = nn - ( 17 * strideZ );
+ ix1 = i4 - ( 2 * strideZ );
+ for ( idx = ( ( nn - offsetZ ) / strideZ ) - 17; idx >= ( 4 * I0 ) + PP + 2; idx -= 4 ) {
+ if ( b2 === 0 ) {
+ break;
+ }
+ b1 = b2;
+ if ( Z[ i4 ] > Z[ ix1 ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b2 *= ( Z[ i4 ] / Z[ ix1 ] );
+ a2 += b2;
+ if ( ( 100 * max( b2, b1 ) ) < a2 || CNST1 < a2 ) {
+ break;
+ }
+ i4 -= inc;
+ ix1 -= inc;
+ }
+ a2 *= CNST3;
+ }
+
+ if ( a2 < CNST1 ) {
+ s = gam * ( 1 - sqrt( a2 ) ) / ( 1 + a2 );
+ }
+ } else {
+ // Case 6, no information to guide us.
+ if ( ttype === -6 ) { // Case when `ttype` is previously set by another routine.
+ G += 0.333 * ( 1 - G );
+ } else if ( ttype === -18 ) { // Case when `ttype` is previously set by another routine.
+ G = 0.25 * 0.333;
+ } else {
+ G = 0.25;
+ }
+ s = G * DMIN;
+ ttype = -6;
+ }
+ } else if ( N0IN === ( N0 + 1 ) ) {
+ // 1 eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN.
+ if ( DMIN1 === DN1 && DMIN2 === DN2 ) {
+ // Cases 7 and 8.
+ ttype = -7;
+ s = 0.3330 * DMIN1;
+
+ if ( Z[ nn - ( 5 * strideZ ) ] > Z[ nn - ( 7 * strideZ ) ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+
+ b1 = Z[ nn - ( 5 * strideZ ) ] / Z[ nn - ( 7 * strideZ ) ];
+ b2 = b1;
+
+ if ( b2 !== 0 ) {
+ i4 = offsetZ + ( strideZ * ( ( 4 * N0 ) - 6 + PP ) );
+ ix1 = i4 - ( 2 * strideZ );
+ for ( idx = 2; idx >= 0; idx-- ) {
+ a2 = b1;
+ if ( Z[ i4 ] > Z[ ix1 ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b1 *= Z[ i4 ] / Z[ ix1 ];
+ b2 += b1;
+
+ if ( ( 100 * max( b1, a2 ) ) < b2 ) {
+ break;
+ }
+ i4 -= inc;
+ ix1 -= inc;
+ }
+ }
+
+ b2 = sqrt( CNST3 * b2 );
+ a2 = DMIN1 / ( 1 + ( b2 * b2 ) );
+ gap2 = ( 0.5 * DMIN2 ) - a2;
+
+ if ( gap2 > 0 && gap2 > ( b2 * a2 ) ) {
+ s = max( s, a2 * ( 1 - ( CNST2 * a2 * ( b2 / gap2 ) * b2 ) ) );
+ } else {
+ s = max( s, a2 * ( 1 - ( CNST2 * b2 ) ) );
+ ttype = -8;
+ }
+ } else {
+ // Case 9.
+ s = 0.25 * DMIN1;
+ if ( DMIN1 === DN1 ) {
+ s = 0.5 * DMIN1;
+ }
+ ttype = -9;
+ }
+ } else if ( N0IN === ( N0 + 2 ) ) {
+ // 2 eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN.
+
+ // Cases 10 and 11.
+ if ( DMIN2 === DN2 && ( 2 * Z[ nn - ( 5 * strideZ ) ] ) < Z[ nn - ( 7 * strideZ ) ] ) {
+ ttype = -10;
+ s = 0.3330 * DMIN2;
+ if ( Z[ nn - ( 5 * strideZ ) ] > Z[ nn - ( 7 * strideZ ) ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b1 = Z[ nn - ( 5 * strideZ ) ] / Z[ nn - ( 7 * strideZ ) ];
+ b2 = b1;
+ if ( b2 !== 0 ) {
+ i4 = offsetZ + ( strideZ * ( ( 4 * N0 ) - 6 + PP ) );
+ ix1 = i4 - ( 2 * strideZ );
+ for ( idx = 1; idx >= 0; idx-- ) {
+ if ( Z[ i4 ] > Z[ ix1 ] ) {
+ idx2 = offsetOut + strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+ }
+ b1 *= ( Z[ i4 ] / Z[ ix1 ] );
+ b2 += b1;
+ if ( ( 100 * b1 ) < b2 ) {
+ break;
+ }
+ i4 -= inc;
+ ix1 -= inc;
+ }
+ }
+
+ b2 = sqrt( CNST3 * b2 );
+ a2 = DMIN2 / ( 1 + ( b2 * b2 ) );
+ gap2 = Z[ nn - ( 7 * strideZ ) ] + Z[ nn - ( 9 * strideZ ) ] - ( sqrt( Z[ nn - ( 11 * strideZ ) ] ) * sqrt( Z[ nn - ( 9 * strideZ ) ] ) ) - a2;
+ if ( gap2 > 0 && gap2 > ( b2 * a2 ) ) {
+ s = max( s, a2 * ( 1 - ( CNST2 * a2 * ( b2 / gap2 ) * b2 ) ) );
+ } else {
+ s = max( s, a2 * ( 1 - ( CNST2 * b2 ) ) );
+ }
+ } else {
+ // Case 11.
+ s = 0.25 * DMIN2;
+ ttype = -11;
+ }
+ } else if ( N0IN > ( N0 + 2 ) ) {
+ // Case 12, more than 2 eigenvalues deflated. No information.
+ s = 0;
+ ttype = -12;
+ }
+
+ tau = s;
+
+ idx2 = offsetOut;
+ out[ idx2 ] = tau;
+ idx2 += strideOut;
+ out[ idx2 ] = ttype;
+ idx2 += strideOut;
+ out[ idx2 ] = G;
+ return out;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq4;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq5.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq5.js
new file mode 100644
index 000000000000..8c938fe35377
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq5.js
@@ -0,0 +1,331 @@
+/**
+* @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 min = require( '@stdlib/math/base/special/min' );
+
+
+// MAIN //
+
+/**
+* Computes one dqds transform in ping-pong form.
+*
+* ## Notes
+*
+* - `Z` is a 1-D array of length >= `4*N0` storing interleaved q/e values.
+* - `PP` is `0` for ping, `1` for pong.
+*
+* @private
+* @param {integer} I0 - the first index
+* @param {integer} N0 - the last index
+* @param {Float64Array} Z - the QD array
+* @param {integer} strideZ - stride length for `Z`
+* @param {NonNegativeInteger} offsetZ - starting index for `Z`
+* @param {integer} PP - ping-pong flag (0 or 1)
+* @param {number} TAU - the shift
+* @param {number} SIGMA - the accumulated shift
+* @param {boolean} IEEE - IEEE arithmetic flag
+* @param {number} EPS - epsilon used by the routine
+* @param {Float64Array} out - output array containing `DMIN`, `DMIN1`, `DMIN2`, `DN`, `DNM1`, and `DNM2` respectively
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index of `out`
+* @returns {Float64Array} output array
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var out = new Float64Array( 6 );
+* var Z = new Float64Array( [ 5, 0, 7, 0, 9, 0, 3, 0, 11, 0, 4, 0, 20, 0, 0, 0 ] );
+*
+* dlasq5( 0, 3, Z, 1, 0, 0, 0.1, 0.0, true, 2.220446049250313e-16, out, 1, 0 );
+* // out => [ ~3.606, ~3.606, ~3.606, ~11.823, ~5.904, ~3.606 ]
+*/
+function dlasq5( I0, N0, Z, strideZ, offsetZ, PP, TAU, SIGMA, IEEE, EPS, out, strideOut, offsetOut ) { // eslint-disable-line max-len, max-params
+ var dthresh;
+ var dmin1;
+ var dmin2;
+ var dnm1;
+ var dnm2;
+ var emin;
+ var dmin;
+ var temp;
+ var j4p2;
+ var idx;
+ var j4;
+ var dn;
+ var d;
+
+ // Quick return...
+ if ( ( N0 - I0 - 1 ) <= 0 ) {
+ return out;
+ }
+
+ dthresh = EPS * ( SIGMA + TAU );
+ if ( TAU < ( dthresh * 0.5 ) ) {
+ TAU = 0;
+ }
+
+ if ( TAU > 0 || TAU < 0 ) {
+ j4 = offsetZ + ( strideZ*( ( 4*I0 ) + PP ) );
+ emin = Z[ j4 + ( 4*strideZ ) ];
+ d = Z[ j4 ] - TAU;
+ dmin = d;
+ dmin1 = -Z[ j4 ];
+
+ if ( IEEE ) {
+ if ( PP === 0 ) {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 2*strideZ ) ] = d + Z[ j4 - strideZ ];
+ temp = Z[ j4 + strideZ ] / Z[ j4 - ( 2*strideZ ) ];
+ d = ( d*temp ) - TAU;
+ dmin = min( dmin, d );
+ Z[ j4 ] = Z[ j4 - strideZ ]*temp;
+ emin = min( Z[ j4 ], emin );
+ j4 += 4*strideZ;
+ }
+ } else {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 3*strideZ ) ] = d + Z[ j4 ];
+ temp = Z[ j4 + ( 2*strideZ ) ] / Z[ j4 - ( 3*strideZ ) ];
+ d = ( d*temp ) - TAU;
+ dmin = min( dmin, d );
+ Z[ j4 - strideZ ] = Z[ j4 ] * temp;
+ emin = min( Z[ j4 - strideZ ], emin );
+ j4 += 4*strideZ;
+ }
+ }
+
+ // Unroll last two steps.
+ dnm2 = d;
+ dmin2 = dmin;
+
+ j4 = offsetZ + ( strideZ * ( ( 4*N0 ) - 5 - PP ) );
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm2 + Z[ j4p2 ];
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dnm1 = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm2 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dnm1 );
+
+ dmin1 = dmin;
+ j4 += 4*strideZ;
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm1 + Z[ j4p2 ];
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dn = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm1 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dn );
+ } else {
+ if ( PP === 0 ) {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 2*strideZ ) ] = d + Z[ j4 - strideZ ];
+ if ( d < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4 + strideZ ] * ( Z[ j4 - strideZ ] / Z[ j4 - ( 2*strideZ ) ] );
+ d = ( Z[ j4 + strideZ ]*( d / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ j4 ] );
+ j4 += 4*strideZ;
+ }
+ } else {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 3*strideZ ) ] = d + Z[ j4 ];
+ if ( d < 0 ) {
+ return out;
+ }
+ Z[ j4 - strideZ ] = Z[ j4 + ( 2*strideZ ) ] * ( Z[ j4 ] / Z[ j4 - ( 3*strideZ ) ] );
+ d = ( Z[ j4 + ( 2*strideZ ) ]*( d / Z[ j4 - ( 3*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ j4 - strideZ ] );
+ j4 += 4*strideZ;
+ }
+ }
+
+ // Unroll last two steps.
+ dnm2 = d;
+ dmin2 = dmin;
+
+ j4 = offsetZ + ( strideZ * ( ( 4*N0 ) - 5 - PP ) );
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm2 + Z[ j4p2 ];
+ if ( dnm2 < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dnm1 = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm2 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dnm1 );
+
+ dmin1 = dmin;
+ j4 += 4*strideZ;
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm1 + Z[ j4p2 ];
+ if ( dnm1 < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dn = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm1 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dn );
+ }
+ } else {
+ // This is the version that sets d's to zero if they are small enough
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + PP ) );
+ emin = Z[ j4 + ( 4*strideZ ) ];
+ d = Z[ j4 ];
+ dmin = d;
+ dmin1 = -Z[ j4 ];
+
+ if ( IEEE ) {
+ if ( PP === 0 ) {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 2*strideZ ) ] = d + Z[ j4 - strideZ ];
+ temp = Z[ j4 + strideZ ] / Z[ j4 - ( 2*strideZ ) ];
+ d = ( d*temp ) - TAU;
+ if ( d < dthresh ) {
+ d = 0;
+ }
+ dmin = min( dmin, d );
+ Z[ j4 ] = Z[ j4 - strideZ ] * temp;
+ emin = min( Z[ j4 ], emin );
+ j4 += 4*strideZ;
+ }
+ } else {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 3*strideZ ) ] = d + Z[ j4 ];
+ temp = Z[ j4 + ( 2*strideZ ) ] / Z[ j4 - ( 3*strideZ ) ];
+ d = ( d*temp ) - TAU;
+ if ( d < dthresh ) {
+ d = 0;
+ }
+ dmin = min( dmin, d );
+ Z[ j4 - strideZ ] = Z[ j4 ] * temp;
+ emin = min( Z[ j4 - strideZ ], emin );
+ j4 += 4*strideZ;
+ }
+ }
+
+ // Unroll last two steps.
+ dnm2 = d;
+ dmin2 = dmin;
+
+ j4 = offsetZ + ( strideZ * ( ( 4*N0 ) - 5 - PP ) );
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm2 + Z[ j4p2 ];
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dnm1 = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm2 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dnm1 );
+
+ dmin1 = dmin;
+ j4 += 4*strideZ;
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm1 + Z[ j4p2 ];
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dn = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm1 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dn );
+ } else {
+ if ( PP === 0 ) {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 2*strideZ ) ] = d + Z[ j4 - strideZ ];
+ if ( d < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4 + strideZ ] * ( Z[ j4 - strideZ ] / Z[ j4 - ( 2*strideZ ) ] );
+ d = ( Z[ j4 + strideZ ]*( d / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ if ( d < dthresh ) {
+ d = 0;
+ }
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ j4 ] );
+ j4 += 4*strideZ;
+ }
+ } else {
+ j4 = offsetZ + ( strideZ * ( ( 4*I0 ) + 3 ) );
+ for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) {
+ Z[ j4 - ( 3*strideZ ) ] = d + Z[ j4 ];
+ if ( d < 0 ) {
+ return out;
+ }
+ Z[ j4 - strideZ ] = Z[ j4 + ( 2*strideZ ) ] * ( Z[ j4 ] / Z[ j4 - ( 3*strideZ ) ] );
+ d = ( Z[ j4 + ( 2*strideZ ) ]*( d / Z[ j4 - ( 3*strideZ ) ] ) ) - TAU;
+ if ( d < dthresh ) {
+ d = 0;
+ }
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ j4 - strideZ ] );
+ j4 += 4*strideZ;
+ }
+ }
+
+ // Unroll last two steps.
+ dnm2 = d;
+ dmin2 = dmin;
+
+ j4 = offsetZ + ( strideZ * ( ( 4*N0 ) - 5 - PP ) );
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm2 + Z[ j4p2 ];
+ if ( dnm2 < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dnm1 = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm2 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+
+ dmin = min( dmin, dnm1 );
+ dmin1 = dmin;
+ j4 += 4*strideZ;
+ j4p2 = j4 + ( strideZ * ( ( 2*PP ) - 1 ) );
+ Z[ j4 - ( 2*strideZ ) ] = dnm1 + Z[ j4p2 ];
+ if ( dnm1 < 0 ) {
+ return out;
+ }
+ Z[ j4 ] = Z[ j4p2 + ( 2*strideZ ) ] * ( Z[ j4p2 ] / Z[ j4 - ( 2*strideZ ) ] );
+ dn = ( Z[ j4p2 + ( 2*strideZ ) ]*( dnm1 / Z[ j4 - ( 2*strideZ ) ] ) ) - TAU;
+ dmin = min( dmin, dn );
+ }
+ }
+
+ Z[ j4 + ( 2*strideZ ) ] = dn;
+ Z[ offsetZ + ( strideZ * ( ( 4*N0 ) + 3 - PP ) ) ] = emin;
+
+ idx = offsetOut;
+ out[ idx ] = dmin;
+ idx += strideOut;
+ out[ idx ] = dmin1;
+ idx += strideOut;
+ out[ idx ] = dmin2;
+ idx += strideOut;
+ out[ idx ] = dn;
+ idx += strideOut;
+ out[ idx ] = dnm1;
+ idx += strideOut;
+ out[ idx ] = dnm2;
+
+ return out;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq5;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq6.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq6.js
new file mode 100644
index 000000000000..29f7e421ac78
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasq6.js
@@ -0,0 +1,219 @@
+/**
+* @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 dlamch = require( '@stdlib/lapack/base/dlamch' );
+var min = require( '@stdlib/math/base/special/min' );
+
+
+// VARIABLES //
+
+var SFMIN = dlamch( 'S' );
+
+
+// MAIN //
+
+/**
+* Computes one DQD transform in ping-pong form.
+*
+* ## Notes
+*
+* - `Z` is a 1-D array of length >= `4 * N0` storing interleaved Q/E values.
+* - `PP` is 0 for ping, 1 for pong.
+*
+* @private
+* @param {integer} I0 - the first index
+* @param {integer} N0 - the last index
+* @param {Float64Array} Z - the QD array
+* @param {integer} strideZ - stride length for `Z`
+* @param {NonNegativeInteger} offsetZ - starting index for `Z`
+* @param {boolean} PP - ping-pong flag (0 or 1)
+* @param {Float64Array} out - output array containing `DMIN`, `DMIN1`, `DMIN2`, `DN`, `DNM1`, and `DNM2` respectively
+* @param {integer} strideOut - stride length for `out`
+* @param {NonNegativeInteger} offsetOut - starting index of `out`
+* @returns {Float64Array} output array
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var out = new Float64Array( 6 );
+* var Z = new Float64Array( [ 5, 0, 7, 0, 9, 0, 3, 0, 11, 0, 4, 0, 20, 0, 0, 0 ] );
+*
+* dlasq6( 0, 3, Z, 1, 0, 0, out, 1, 0 );
+* // out => [ 3.75, 3.75, 3.75, ~12.088, ~6.111, 3.75 ]
+*/
+function dlasq6( I0, N0, Z, strideZ, offsetZ, PP, out, strideOut, offsetOut ) {
+ var dmin1;
+ var dmin2;
+ var dnm1;
+ var dnm2;
+ var emin;
+ var dmin;
+ var temp;
+ var j4p2;
+ var idx;
+ var ix1;
+ var ix2;
+ var ix3;
+ var ix4;
+ var j4;
+ var dn;
+ var d;
+
+ // Quick return...
+ if ( ( N0 - I0 - 1 ) <= 0 ) {
+ return out;
+ }
+
+ idx = offsetZ;
+
+ j4 = offsetZ + ( strideZ * ( ( 4 * I0 ) + PP ) );
+ emin = Z[ j4 + ( 4 * strideZ ) ];
+ d = Z[ j4 ];
+ dmin = d;
+
+ j4 = offsetZ + ( strideZ * ( ( 4 * I0 ) + 3 ) );
+
+ if ( PP === 0 ) {
+ for ( idx = ( 4 * I0 ) + 3; idx <= ( 4 * N0 ) - 9; idx += 4 ) {
+ ix1 = j4 + strideZ;
+ ix2 = j4 - strideZ;
+ ix3 = j4 - ( 2 * strideZ );
+
+ Z[ ix3 ] = d + Z[ ix2 ];
+ if ( Z[ ix3 ] === 0 ) {
+ Z[ j4 ] = 0;
+ d = Z[ ix1 ];
+ dmin = d;
+ emin = 0;
+ } else if ( ( SFMIN * Z[ ix1 ] < Z[ ix3 ] ) &&
+ ( SFMIN * Z[ ix3 ] < Z[ ix1 ] ) ) {
+ temp = Z[ ix1 ] / Z[ ix3 ];
+ Z[ j4 ] = Z[ ix2 ] * temp;
+ d *= temp;
+ } else {
+ Z[ j4 ] = Z[ ix1 ] * ( Z[ ix2 ] / Z[ ix3 ] );
+ d = Z[ ix1 ] * ( d / Z[ ix3 ] );
+ }
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ j4 ] );
+ j4 += 4 * strideZ;
+ }
+ } else {
+ for ( idx = ( 4 * I0 ) + 3; idx <= ( 4 * N0 ) - 9; idx += 4 ) {
+ ix1 = j4 + ( 2 * strideZ );
+ ix2 = j4 - strideZ;
+ ix3 = j4 - ( 2 * strideZ );
+ ix4 = j4 - ( 3 * strideZ );
+
+ Z[ ix4 ] = d + Z[ j4 ];
+ if ( Z[ ix4 ] === 0 ) {
+ Z[ ix2 ] = 0;
+ d = Z[ ix1 ];
+ dmin = d;
+ emin = 0;
+ } else if ( ( SFMIN * Z[ ix1 ] < Z[ ix4 ] ) &&
+ ( SFMIN * Z[ ix4 ] < Z[ ix1 ] ) ) {
+ temp = Z[ ix1 ] / Z[ ix4 ];
+ Z[ ix2 ] = Z[ j4 ] * temp;
+ d *= temp;
+ } else {
+ Z[ ix2 ] = Z[ ix1 ] * ( Z[ j4 ] / Z[ ix4 ] );
+ d = Z[ ix1 ] * ( d / Z[ ix4 ] );
+ }
+ dmin = min( dmin, d );
+ emin = min( emin, Z[ ix2 ] );
+ j4 += 4 * strideZ;
+ }
+ }
+
+ // Unroll last two steps.
+ dnm2 = d;
+ dmin2 = dmin;
+ j4 = offsetZ + ( strideZ * ( ( 4 * N0 ) - 5 - PP ) );
+ j4p2 = j4 + ( strideZ * ( ( 2 * PP ) - 1 ) );
+
+ ix1 = j4 - ( 2 * strideZ );
+ ix2 = j4p2 + ( 2 * strideZ );
+
+ Z[ ix1 ] = dnm2 + Z[ j4p2 ];
+ if ( Z[ ix1 ] === 0 ) {
+ Z[ j4 ] = 0;
+ dnm1 = Z[ ix2 ];
+ dmin = dnm1;
+ emin = 0;
+ } else if ( ( SFMIN * Z[ ix2 ] < Z[ ix1 ] ) && ( SFMIN * Z[ ix1 ] < Z[ ix2 ] ) ) {
+ temp = Z[ ix2 ] / Z[ ix1 ];
+ Z[ j4 ] = Z[ j4p2 ] * temp;
+ dnm1 = dnm2 * temp;
+ } else {
+ Z[ j4 ] = Z[ ix2 ] * ( Z[ j4p2 ] / Z[ ix1 ] );
+ dnm1 = Z[ ix2 ] * ( dnm2 / Z[ ix1 ] );
+ }
+ dmin = min( dmin, dnm1 );
+
+ dmin1 = dmin;
+ j4 += 4 * strideZ;
+ j4p2 += 4 * strideZ;
+
+ ix1 = j4 - ( 2 * strideZ );
+ ix2 = j4p2 + ( 2 * strideZ );
+
+ Z[ ix1 ] = dnm1 + Z[ j4p2 ];
+ if ( Z[ ix1 ] === 0 ) {
+ Z[ j4 ] = 0;
+ dn = Z[ ix2 ];
+ dmin = dn;
+ emin = 0;
+ } else if ( ( SFMIN * Z[ ix2 ] < Z[ ix1 ] ) &&
+ ( SFMIN * Z[ ix1 ] < Z[ ix2 ] ) ) {
+ temp = Z[ ix2 ] / Z[ ix1 ];
+ Z[ j4 ] = Z[ j4p2 ] * temp;
+ dn = dnm1 * temp;
+ } else {
+ Z[ j4 ] = Z[ ix2 ] * ( Z[ j4p2 ] / Z[ ix1 ] );
+ dn = Z[ ix2 ] * ( dnm1 / Z[ ix1 ] );
+ }
+ dmin = min( dmin, dn );
+
+ Z[ j4 + ( 2 * strideZ ) ] = dn;
+ Z[ offsetZ + ( strideZ * ( ( 4 * ( N0 ) ) + 3 - PP ) ) ] = emin;
+
+ idx = offsetOut;
+ out[ idx ] = dmin;
+ idx += strideOut;
+ out[ idx ] = dmin1;
+ idx += strideOut;
+ out[ idx ] = dmin2;
+ idx += strideOut;
+ out[ idx ] = dn;
+ idx += strideOut;
+ out[ idx ] = dnm1;
+ idx += strideOut;
+ out[ idx ] = dnm2;
+
+ return out;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq6;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasrt.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasrt.js
new file mode 100644
index 000000000000..01ad9afb80a9
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/dlasrt.js
@@ -0,0 +1,250 @@
+/**
+* @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 lowercase = require( '@stdlib/string/base/lowercase' );
+var floor = require( '@stdlib/math/base/special/floor' );
+var Float64Array = require( '@stdlib/array/float64' );
+
+
+// VARIABLES //
+
+var SELECT = 20;
+
+
+// MAIN //
+
+/**
+* Sort an array of doubles in increasing or decreasing order using quicksort, with insertion sort for small partitions (size <= 20).
+*
+* @private
+* @param {string} ID - sort direction: 'I' or 'D'
+* @param {NonNegativeInteger} N - number of elements to sort
+* @param {Float64Array} D - array to sort in-place
+* @param {integer} strideD - stride length for `d`
+* @param {NonNegativeInteger} offsetD - starting index for `d`
+* @returns {integer} status code (0 = success, -i = the i-th argument had an illegal value)
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 5, 7, 9, 3, 11 ] );
+*
+* dlasrt( 'I', 5, D, 1, 0 );
+* // D => [ 3, 5, 7, 9, 11 ]
+*/
+function dlasrt( ID, N, D, strideD, offsetD ) {
+ var stkpnt;
+ var stack;
+ var start;
+ var dmnmx;
+ var endd;
+ var idx1;
+ var idx2;
+ var idx3;
+ var dir;
+ var tmp;
+ var d1;
+ var d2;
+ var d3;
+ var i;
+ var j;
+
+ // Determine sort direction
+ dir = -1;
+ if ( lowercase( ID ) === 'd' ) {
+ dir = 0;
+ } else if ( lowercase( ID ) === 'i' ) {
+ dir = 1;
+ }
+
+ if ( dir === -1 ) {
+ return -1;
+ }
+ if ( N < 0 ) {
+ return -2;
+ }
+
+ // Quick return
+ if ( N <= 1 ) {
+ return 0;
+ }
+
+ // Initialize the stack with the full range (using 0-based indices)...
+ stkpnt = 0;
+ stack = new Float64Array( 64 );
+ stack[ 0 ] = 0;
+ stack[ 1 ] = N - 1;
+
+ while ( stkpnt >= 0 ) {
+ start = stack[ 2*stkpnt ];
+ endd = stack[ ( 2*stkpnt ) + 1 ];
+ stkpnt -= 1;
+
+ if ( endd - start <= SELECT && endd - start > 0 ) {
+ // Insertion sort for small partitions
+ if ( dir === 0 ) {
+ // Sort in decreasing order
+ idx1 = offsetD + ( ( start + 1 )*strideD );
+ for ( i = start + 1; i <= endd; i++ ) {
+ idx2 = idx1; // index of i
+ for ( j = i; j >= start + 1; j-- ) {
+ idx3 = idx2 - strideD; // index of j - 1
+ if ( D[ idx2 ] > D[ idx3 ] ) {
+ dmnmx = D[ idx2 ];
+ D[ idx2 ] = D[ idx3 ];
+ D[ idx3 ] = dmnmx;
+ } else {
+ break;
+ }
+ idx2 -= strideD;
+ }
+ idx1 += strideD;
+ }
+ } else {
+ // Sort in increasing order
+ idx1 = offsetD + ( ( start + 1 )*strideD );
+ for ( i = start + 1; i <= endd; i++ ) {
+ idx2 = idx1; // index of i
+ for ( j = i; j >= start + 1; j-- ) {
+ idx3 = idx2 - strideD; // index of j - 1
+ if ( D[ idx2 ] < D[ idx3 ] ) {
+ dmnmx = D[ idx2 ];
+ D[ idx2 ] = D[ idx3 ];
+ D[ idx3 ] = dmnmx;
+ } else {
+ break;
+ }
+ idx2 -= strideD;
+ }
+ idx1 += strideD;
+ }
+ }
+ } else if ( endd - start > SELECT ) {
+ // Quicksort partition using median-of-three pivot
+ d1 = D[ offsetD + ( start*strideD ) ];
+ d2 = D[ offsetD + ( endd*strideD ) ];
+ i = floor( ( start + endd ) / 2 );
+ d3 = D[ offsetD + ( i*strideD ) ];
+
+ // Find median of d1, d2, d3
+ if ( d1 < d2 ) {
+ if ( d3 < d1 ) {
+ dmnmx = d1;
+ } else if ( d3 < d2 ) {
+ dmnmx = d3;
+ } else {
+ dmnmx = d2;
+ }
+ } else if ( d3 < d2 ) {
+ dmnmx = d2;
+ } else if ( d3 < d1 ) {
+ dmnmx = d3;
+ } else {
+ dmnmx = d1;
+ }
+
+ if ( dir === 0 ) {
+ // Partition for decreasing order
+ i = start;
+ j = endd;
+ while ( true ) {
+ idx1 = offsetD + ( j*strideD );
+ while ( D[ idx1 ] < dmnmx ) {
+ idx1 -= strideD;
+ j -= 1;
+ }
+
+ idx2 = offsetD + ( i*strideD );
+ while ( D[ idx2 ] > dmnmx ) {
+ idx2 += strideD;
+ i += 1;
+ }
+ if ( i < j ) {
+ tmp = D[ idx2 ];
+ D[ idx2 ] = D[ idx1 ];
+ D[ idx1 ] = tmp;
+ i += 1;
+ j -= 1;
+ } else {
+ break;
+ }
+ }
+ } else {
+ // Partition for increasing order
+ i = start;
+ j = endd;
+ while ( true ) {
+ idx1 = offsetD + ( j*strideD );
+ while ( D[ idx1 ] > dmnmx ) {
+ idx1 -= strideD;
+ j -= 1;
+ }
+
+ idx2 = offsetD + ( i*strideD );
+ while ( D[ idx2 ] < dmnmx ) {
+ idx2 += strideD;
+ i += 1;
+ }
+ if ( i < j ) {
+ tmp = D[ idx2 ];
+ D[ idx2 ] = D[ idx1 ];
+ D[ idx1 ] = tmp;
+ i += 1;
+ j -= 1;
+ } else {
+ break;
+ }
+ }
+ }
+
+ // Push sub-partitions onto stack (larger first for bounded stack depth)
+ idx1 = 2*stkpnt;
+ if ( j - start > endd - j - 1 ) {
+ stkpnt += 1;
+ idx1 += 2;
+ stack[ idx1 ] = start;
+ stack[ idx1 + 1 ] = j;
+
+ stkpnt += 1;
+ idx1 += 2;
+ stack[ idx1 ] = j + 1;
+ stack[ idx1 + 1 ] = endd;
+ } else {
+ stkpnt += 1;
+ idx1 += 2;
+ stack[ idx1 ] = j + 1;
+ stack[ idx1 + 1 ] = endd;
+
+ stkpnt += 1;
+ idx1 += 2;
+ stack[ idx1 ] = start;
+ stack[ idx1 + 1 ] = j;
+ }
+ }
+ }
+ return 0;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasrt;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/index.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/index.js
new file mode 100644
index 000000000000..5fafd6748cb6
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/index.js
@@ -0,0 +1,80 @@
+/**
+* @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';
+
+/**
+* Compute the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E`.
+*
+* @module @stdlib/lapack/base/dlasq1
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+* var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+*
+* var D = new Float64Array( [ 100, 50, 25 ] );
+* var E = new Float64Array( [ 90, 40, 0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* var info = dlasq1( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90, 40, 0 ]
+* // WORK => [ ~1.939e292, ~3.137e291, ~2.554e290, ~1.575e286, ~3.137e291, ~2.881e291, ~2.278e292, ~2.278e292, 7.0, ~2.667, 0.0, ~2.881e291 ]
+* // info => 0
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+* var dlasq1 = require( '@stdlib/lapack/base/dlasq1' );
+*
+* var D = new Float64Array( [ 100, 50, 25 ] );
+* var E = new Float64Array( [ 90, 40, 0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* var info = dlasq1( 3, D, E, WORK );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90, 40, 0 ]
+* // WORK => [ ~1.939e292, ~3.137e291, ~2.554e290, ~1.575e286, ~3.137e291, ~2.881e291, ~2.278e292, ~2.278e292, 7.0, ~2.667, 0.0, ~2.881e291 ]
+* // info => 0
+*
+*/
+
+
+// MODULES //
+
+var join = require( 'path' ).join;
+var tryRequire = require( '@stdlib/utils/try-require' );
+var isError = require( '@stdlib/assert/is-error' );
+var main = require( './main.js' );
+
+
+// MAIN //
+
+var dlasq1;
+var tmp = tryRequire( join( __dirname, './native.js' ) );
+if ( isError( tmp ) ) {
+ dlasq1 = main;
+} else {
+ dlasq1 = tmp;
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq1;
+
+// exports: { "ndarray": "dlasq1.ndarray" }
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/main.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/main.js
new file mode 100644
index 000000000000..f6d2b228e5b6
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/main.js
@@ -0,0 +1,35 @@
+/**
+* @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 setReadOnly = require( '@stdlib/utils/define-nonenumerable-read-only-property' );
+var dlasq1 = require( './dlasq1.js' );
+var ndarray = require( './ndarray.js' );
+
+
+// MAIN //
+
+setReadOnly( dlasq1, 'ndarray', ndarray );
+
+
+// EXPORTS //
+
+module.exports = dlasq1;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/ndarray.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/ndarray.js
new file mode 100644
index 000000000000..e377024fee92
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/lib/ndarray.js
@@ -0,0 +1,80 @@
+/**
+* @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 base = require( './base.js' );
+
+
+// MAIN //
+
+/**
+* Computes the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E` using alternative indexing semantics.
+*
+* ## Notes
+*
+* - `D` should have `N` indexed elements. On entry, `D` contains the diagonal elements of the bi-diagonal matrix whose SVD is desired. On normal exit, `D` contains the singular values in decreasing order.
+*
+* - `E` should have `N` indexed elements. On entry, elements E(1:N-1) contain the off-diagonal elements of the bi-diagonal matrix whose SVD is desired. On exit, E is overwritten.
+*
+* - The function returns a status code:
+*
+* - `= 0`: successful exit.
+*
+* - `< 0`: if `INFO = -i`, the `i`-th argument had an illegal value.
+*
+* - `> 0`: the algorithm failed:
+* - `= 1`, a split was marked by a positive value in `E`.
+* - `= 2`, current block of `Z` not diagonalized after `100*N` iterations (in inner while loop). On exit, `D` and `E` represent a matrix with the same singular values which the calling subroutine could use to finish the computation, or even feed back into `DLASQ1`
+* - `= 3`, termination criterion of outer while loop not met (program created more than `N` unreduced blocks)
+*
+* @param {integer} N - number of rows/columns in the matrix
+* @param {Float64Array} D - the array with diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {integer} strideD - stride length for `D`
+* @param {NonNegativeInteger} offsetD - starting index of `D`
+* @param {Float64Array} E - the array with off-diagonal elements of the bi-diagonal matrix whose SVD is desired
+* @param {integer} strideE - stride length for `E`
+* @param {NonNegativeInteger} offsetE - starting index of `E`
+* @param {Float64Array} WORK - workspace array (length >= 4*N)
+* @param {integer} strideWORK - stride length for `WORK`
+* @param {NonNegativeInteger} offsetWORK - starting index of `WORK`
+* @returns {integer} status code
+*
+* @example
+* var Float64Array = require( '@stdlib/array/float64' );
+*
+* var D = new Float64Array( [ 100, 50, 25 ] );
+* var E = new Float64Array( [ 90, 40, 0 ] );
+* var WORK = new Float64Array( 12 );
+*
+* var info = dlasq1( 3, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+* // D => [ ~139.377, ~56.064, ~15.997 ]
+* // E => [ 90, 40, 0 ]
+* // WORK => [ ~1.939e292, ~3.137e291, ~2.554e290, ~1.575e286, ~3.137e291, ~2.881e291, ~2.278e292, ~2.278e292, 7.0, ~2.667, 0.0, ~2.881e291 ]
+* // info => 0
+*/
+function dlasq1( N, D, strideD, offsetD, E, strideE, offsetE, WORK, strideWORK, offsetWORK ) { // eslint-disable-line stdlib/jsdoc-doctest-decimal-point
+ return base( N, D, strideD, offsetD, E, strideE, offsetE, WORK, strideWORK, offsetWORK );
+}
+
+
+// EXPORTS //
+
+module.exports = dlasq1;
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/package.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/package.json
new file mode 100644
index 000000000000..77fd1cff82ef
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/package.json
@@ -0,0 +1,72 @@
+{
+ "name": "@stdlib/lapack/base/dlasq1",
+ "version": "0.0.0",
+ "description": "LAPACK routine to computes the singular values of a real `N-by-N` bidiagonal matrix with diagonal `D` and off-diagonal `E`.",
+ "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",
+ "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",
+ "lapack",
+ "svd",
+ "decomposition",
+ "dlasq1",
+ "exchange",
+ "permute",
+ "permutedims",
+ "linear",
+ "algebra",
+ "subroutines",
+ "array",
+ "ndarray",
+ "float64",
+ "double",
+ "float64array"
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/info_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/info_eq_0.json
new file mode 100644
index 000000000000..7d98b2a9a2ad
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/info_eq_0.json
@@ -0,0 +1,45 @@
+{
+ "N": 3,
+
+ "D": [ 100.0, 50.0, 25.0 ],
+ "strideD": 1,
+ "offsetD": 0,
+
+
+ "E": [ 90.0, 40.0, 0.0 ],
+ "strideE": 1,
+ "offsetE": 0,
+
+ "WORK": [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ],
+ "strideWORK": 1,
+ "offsetWORK": 0,
+
+ "info": 0,
+
+ "D_out": [
+ 1.39376808415966082E+002,
+ 5.60643189944822424E+001,
+ 1.59967937870069168E+001
+ ],
+
+ "E_out": [
+ 9.00000000000000000E+001,
+ 4.00000000000000000E+001,
+ 0.00000000000000000E+000
+ ],
+
+ "WORK_out": [
+ 1.93854918696889742E+292,
+ 3.13667047842326929E+291,
+ 2.55365184452744672E+290,
+ 1.57480965220002271E+286,
+ 3.13667047842326929E+291,
+ 2.88130529397052490E+291,
+ 2.27775275325649918E+292,
+ 2.27775275325649874E+292,
+ 7.00000000000000000E+000,
+ 2.66666666666666652E+000,
+ 0.00000000000000000E+000,
+ 2.88132383458501245E+291
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/info_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/info_eq_0.json
new file mode 100644
index 000000000000..7d64b7d50129
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/info_eq_0.json
@@ -0,0 +1,94 @@
+{
+ "N": 3,
+ "D": [
+ 100,
+ 9999,
+ 50,
+ 9999,
+ 25,
+ 9999
+ ],
+ "strideD": 2,
+ "offsetD": 0,
+ "E": [
+ 90,
+ 9999,
+ 40,
+ 9999,
+ 0,
+ 9999
+ ],
+ "strideE": 2,
+ "offsetE": 0,
+ "WORK": [
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999
+ ],
+ "strideWORK": 2,
+ "offsetWORK": 0,
+ "info": 0,
+ "D_out": [
+ 139.37680841596608,
+ 9999,
+ 56.06431899448224,
+ 9999,
+ 15.996793787006917,
+ 9999
+ ],
+ "E_out": [
+ 90,
+ 9999,
+ 40,
+ 9999,
+ 0,
+ 9999
+ ],
+ "WORK_out": [
+ 1.9385491869688974e+292,
+ 9999,
+ 3.1366704784232693e+291,
+ 9999,
+ 2.5536518445274467e+290,
+ 9999,
+ 1.5748096522000227e+286,
+ 9999,
+ 3.1366704784232693e+291,
+ 9999,
+ 2.881305293970525e+291,
+ 9999,
+ 2.277752753256499e+292,
+ 9999,
+ 2.2777527532564987e+292,
+ 9999,
+ 7,
+ 9999,
+ 2.6666666666666665,
+ 9999,
+ 0,
+ 9999,
+ 2.8813238345850124e+291,
+ 9999
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/sigmx_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/sigmx_eq_0.json
new file mode 100644
index 000000000000..0191a798592a
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/large_strides/sigmx_eq_0.json
@@ -0,0 +1,94 @@
+{
+ "N": 3,
+ "D": [
+ 5,
+ 9999,
+ 3,
+ 9999,
+ 7,
+ 9999
+ ],
+ "strideD": 2,
+ "offsetD": 0,
+ "E": [
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999
+ ],
+ "strideE": 2,
+ "offsetE": 0,
+ "WORK": [
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999
+ ],
+ "strideWORK": 2,
+ "offsetWORK": 0,
+ "info": 0,
+ "D_out": [
+ 7,
+ 9999,
+ 5,
+ 9999,
+ 3,
+ 9999
+ ],
+ "E_out": [
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999
+ ],
+ "WORK_out": [
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999,
+ 0,
+ 9999
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/info_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/info_eq_0.json
new file mode 100644
index 000000000000..3c57cb07ab95
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/info_eq_0.json
@@ -0,0 +1,58 @@
+{
+ "N": 3,
+ "D": [
+ 25,
+ 50,
+ 100
+ ],
+ "strideD": -1,
+ "offsetD": 2,
+ "E": [
+ 0,
+ 40,
+ 90
+ ],
+ "strideE": -1,
+ "offsetE": 2,
+ "WORK": [
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ],
+ "strideWORK": -1,
+ "offsetWORK": 11,
+ "info": 0,
+ "D_out": [
+ 15.996793787006917,
+ 56.06431899448224,
+ 139.37680841596608
+ ],
+ "E_out": [
+ 0,
+ 40,
+ 90
+ ],
+ "WORK_out": [
+ 2.8813238345850124e+291,
+ 0,
+ 2.6666666666666665,
+ 7,
+ 2.2777527532564987e+292,
+ 2.277752753256499e+292,
+ 2.881305293970525e+291,
+ 3.1366704784232693e+291,
+ 1.5748096522000227e+286,
+ 2.5536518445274467e+290,
+ 3.1366704784232693e+291,
+ 1.9385491869688974e+292
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/sigmx_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/sigmx_eq_0.json
new file mode 100644
index 000000000000..546f2abbfbdc
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/negative_strides/sigmx_eq_0.json
@@ -0,0 +1,58 @@
+{
+ "N": 3,
+ "D": [
+ 7,
+ 3,
+ 5
+ ],
+ "strideD": -1,
+ "offsetD": 2,
+ "E": [
+ 0,
+ 0,
+ 0
+ ],
+ "strideE": -1,
+ "offsetE": 2,
+ "WORK": [
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ],
+ "strideWORK": -1,
+ "offsetWORK": 11,
+ "info": 0,
+ "D_out": [
+ 3,
+ 5,
+ 7
+ ],
+ "E_out": [
+ 0,
+ 0,
+ 0
+ ],
+ "WORK_out": [
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/info_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/info_eq_0.json
new file mode 100644
index 000000000000..6c17a590cc50
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/info_eq_0.json
@@ -0,0 +1,64 @@
+{
+ "N": 3,
+ "D": [
+ 9999,
+ 100,
+ 50,
+ 25
+ ],
+ "strideD": 1,
+ "offsetD": 1,
+ "E": [
+ 9999,
+ 90,
+ 40,
+ 0
+ ],
+ "strideE": 1,
+ "offsetE": 1,
+ "WORK": [
+ 9999,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ],
+ "strideWORK": 1,
+ "offsetWORK": 1,
+ "info": 0,
+ "D_out": [
+ 9999,
+ 139.37680841596608,
+ 56.06431899448224,
+ 15.996793787006917
+ ],
+ "E_out": [
+ 9999,
+ 90,
+ 40,
+ 0
+ ],
+ "WORK_out": [
+ 9999,
+ 1.9385491869688974e+292,
+ 3.1366704784232693e+291,
+ 2.5536518445274467e+290,
+ 1.5748096522000227e+286,
+ 3.1366704784232693e+291,
+ 2.881305293970525e+291,
+ 2.277752753256499e+292,
+ 2.2777527532564987e+292,
+ 7,
+ 2.6666666666666665,
+ 0,
+ 2.8813238345850124e+291
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/sigmx_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/sigmx_eq_0.json
new file mode 100644
index 000000000000..be5100cb8f18
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/offsets/sigmx_eq_0.json
@@ -0,0 +1,64 @@
+{
+ "N": 3,
+ "D": [
+ 9999,
+ 5,
+ 3,
+ 7
+ ],
+ "strideD": 1,
+ "offsetD": 1,
+ "E": [
+ 9999,
+ 0,
+ 0,
+ 0
+ ],
+ "strideE": 1,
+ "offsetE": 1,
+ "WORK": [
+ 9999,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ],
+ "strideWORK": 1,
+ "offsetWORK": 1,
+ "info": 0,
+ "D_out": [
+ 9999,
+ 7,
+ 5,
+ 3
+ ],
+ "E_out": [
+ 9999,
+ 0,
+ 0,
+ 0
+ ],
+ "WORK_out": [
+ 9999,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0,
+ 0
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/sigmx_eq_0.json b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/sigmx_eq_0.json
new file mode 100644
index 000000000000..4dce29f0c6e9
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/fixtures/sigmx_eq_0.json
@@ -0,0 +1,44 @@
+{
+ "N": 3,
+
+ "D": [ 5.0, 3.0, 7.0 ],
+ "strideD": 1,
+ "offsetD": 0,
+
+ "E": [ 0.0, 0.0, 0.0 ],
+ "strideE": 1,
+ "offsetE": 0,
+
+ "WORK": [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ],
+ "strideWORK": 1,
+ "offsetWORK": 0,
+
+ "info": 0,
+
+ "D_out": [
+ 7.00000000000000000E+000,
+ 5.00000000000000000E+000,
+ 3.00000000000000000E+000
+ ],
+
+ "E_out": [
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000
+ ],
+
+ "WORK_out": [
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000,
+ 0.00000000000000000E+000
+ ]
+}
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.dlasq1.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.dlasq1.js
new file mode 100644
index 000000000000..4e2eae87a498
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.dlasq1.js
@@ -0,0 +1,173 @@
+/*
+* @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';
+
+var tape = require( 'tape' );
+var Float64Array = require( '@stdlib/array/float64' );
+var isAlmostEqual = require( '@stdlib/assert/is-almost-equal-float64array' );
+var dlasq1 = require( './../lib/dlasq1.js' );
+
+
+// FIXTURES //
+
+// All tests have been verified with the LAPACK fortran code.
+var INFO_EQ_0 = require( './fixtures/info_eq_0.json' );
+var SIGMX_EQ_0 = require( './fixtures/sigmx_eq_0.json' );
+
+
+// TESTS //
+
+tape( 'main export is a function', function test( t ) {
+ t.ok( true, __filename );
+ t.strictEqual( typeof dlasq1, 'function', 'main export is a function' );
+ t.end();
+});
+
+tape( 'the function has an arity of 4', function test( t ) {
+ t.strictEqual( dlasq1.length, 4, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=0)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+ E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+ expectedWORK = WORK;
+ expectedD = D;
+ expectedE = E;
+
+ t.strictEqual( dlasq1( 0, D, E, WORK ), 0, 'returns expected value' );
+ t.strictEqual( D, expectedD, 'returns expected value' );
+ t.strictEqual( E, expectedE, 'returns expected value' );
+ t.strictEqual( WORK, expectedWORK, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=1)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ -100.0 ] );
+ E = new Float64Array( [ 90.0 ] );
+ expectedWORK = WORK;
+ expectedD = new Float64Array( [ 100.0 ] );
+ expectedE = E;
+
+ t.strictEqual( dlasq1( 1, D, E, WORK ), 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( E, expectedE, 'returns expected value' );
+ t.strictEqual( WORK, expectedWORK, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=2)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ 100.0, 50.0 ] );
+ E = new Float64Array( [ 90.0, 40.0 ] );
+
+ expectedWORK = WORK;
+ expectedD = new Float64Array( [ 138.9424291276145, 35.98612771774451 ] );
+ expectedE = E;
+
+ info = dlasq1( 2, D, E, WORK );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function computes the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E`', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var data;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ data = INFO_EQ_0;
+
+ D = new Float64Array( data.D );
+ E = new Float64Array( data.E );
+ WORK = new Float64Array( data.WORK );
+
+ expectedD = new Float64Array( data.D_out );
+ expectedE = new Float64Array( data.E_out );
+ expectedWORK = new Float64Array( data.WORK_out );
+
+ info = dlasq1( data.N, D, E, WORK );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function sorts the diagonal elements in descending order when the matrix is already diagonal (sigmx==0)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var data;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ data = SIGMX_EQ_0;
+
+ D = new Float64Array( data.D );
+ E = new Float64Array( data.E );
+ WORK = new Float64Array( data.WORK );
+
+ expectedD = new Float64Array( data.D_out );
+ expectedE = new Float64Array( data.E_out );
+ expectedWORK = new Float64Array( data.WORK_out );
+
+ info = dlasq1( data.N, D, E, WORK );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.js
new file mode 100644
index 000000000000..8eb4bbd1edf3
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.js
@@ -0,0 +1,82 @@
+/**
+* @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 tape = require( 'tape' );
+var proxyquire = require( 'proxyquire' );
+var IS_BROWSER = require( '@stdlib/assert/is-browser' );
+var dlasq1 = require( './../lib' );
+
+
+// VARIABLES //
+
+var opts = {
+ 'skip': IS_BROWSER
+};
+
+
+// TESTS //
+
+tape( 'main export is a function', function test( t ) {
+ t.ok( true, __filename );
+ t.strictEqual( typeof dlasq1, 'function', 'main export is a function' );
+ t.end();
+});
+
+tape( 'attached to the main export is a method providing an ndarray interface', function test( t ) {
+ t.strictEqual( typeof dlasq1.ndarray, 'function', 'method is a function' );
+ t.end();
+});
+
+tape( 'if a native implementation is available, the main export is the native implementation', opts, function test( t ) {
+ var dlasq1 = proxyquire( './../lib', {
+ '@stdlib/utils/try-require': tryRequire
+ });
+
+ t.strictEqual( dlasq1, mock, 'returns expected value' );
+ t.end();
+
+ function tryRequire() {
+ return mock;
+ }
+
+ function mock() {
+ // Mock...
+ }
+});
+
+tape( 'if a native implementation is not available, the main export is a JavaScript implementation', opts, function test( t ) {
+ var dlasq1;
+ var main;
+
+ main = require( './../lib/dlasq1.js' );
+
+ dlasq1 = proxyquire( './../lib', {
+ '@stdlib/utils/try-require': tryRequire
+ });
+
+ t.strictEqual( dlasq1, main, 'returns expected value' );
+ t.end();
+
+ function tryRequire() {
+ return new Error( 'Cannot find module' );
+ }
+});
diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.ndarray.js b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.ndarray.js
new file mode 100644
index 000000000000..e91fb8deb15e
--- /dev/null
+++ b/lib/node_modules/@stdlib/lapack/base/dlasq1/test/test.ndarray.js
@@ -0,0 +1,173 @@
+/*
+* @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';
+
+var tape = require( 'tape' );
+var Float64Array = require( '@stdlib/array/float64' );
+var isAlmostEqual = require( '@stdlib/assert/is-almost-equal-float64array' );
+var dlasq1 = require( './../lib/ndarray.js' );
+
+
+// FIXTURES //
+
+// All tests have been verified with the LAPACK fortran code.
+var INFO_EQ_0 = require( './fixtures/info_eq_0.json' );
+var SIGMX_EQ_0 = require( './fixtures/sigmx_eq_0.json' );
+
+
+// TESTS //
+
+tape( 'main export is a function', function test( t ) {
+ t.ok( true, __filename );
+ t.strictEqual( typeof dlasq1, 'function', 'main export is a function' );
+ t.end();
+});
+
+tape( 'the function has an arity of 10', function test( t ) {
+ t.strictEqual( dlasq1.length, 10, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=0)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ 100.0, 50.0, 25.0 ] );
+ E = new Float64Array( [ 90.0, 40.0, 0.0 ] );
+ expectedWORK = WORK;
+ expectedD = D;
+ expectedE = E;
+
+ t.strictEqual( dlasq1( 0, D, 1, 0, E, 1, 0, WORK, 1, 0 ), 0, 'returns expected value' );
+ t.strictEqual( D, expectedD, 'returns expected value' );
+ t.strictEqual( E, expectedE, 'returns expected value' );
+ t.strictEqual( WORK, expectedWORK, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=1)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ -100.0 ] );
+ E = new Float64Array( [ 90.0 ] );
+ expectedWORK = WORK;
+ expectedD = new Float64Array( [ 100.0 ] );
+ expectedE = E;
+
+ t.strictEqual( dlasq1( 1, D, 1, 0, E, 1, 0, WORK, 1, 0 ), 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( E, expectedE, 'returns expected value' );
+ t.strictEqual( WORK, expectedWORK, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function quick returns (N=2)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ WORK = new Float64Array( 10 );
+ D = new Float64Array( [ 100.0, 50.0 ] );
+ E = new Float64Array( [ 90.0, 40.0 ] );
+
+ expectedWORK = WORK;
+ expectedD = new Float64Array( [ 138.9424291276145, 35.98612771774451 ] );
+ expectedE = E;
+
+ info = dlasq1( 2, D, 1, 0, E, 1, 0, WORK, 1, 0 );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function computes the singular values of a real `N-by-N` bi-diagonal matrix with diagonal `D` and off-diagonal `E`', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var data;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ data = INFO_EQ_0;
+
+ D = new Float64Array( data.D );
+ E = new Float64Array( data.E );
+ WORK = new Float64Array( data.WORK );
+
+ expectedD = new Float64Array( data.D_out );
+ expectedE = new Float64Array( data.E_out );
+ expectedWORK = new Float64Array( data.WORK_out );
+
+ info = dlasq1( data.N, D, data.strideD, data.offsetD, E, data.strideE, data.offsetE, WORK, data.strideWORK, data.offsetWORK );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});
+
+tape( 'the function sorts the diagonal elements in descending order when the matrix is already diagonal (sigmx==0)', function test( t ) {
+ var expectedWORK;
+ var expectedD;
+ var expectedE;
+ var data;
+ var info;
+ var WORK;
+ var D;
+ var E;
+
+ data = SIGMX_EQ_0;
+
+ D = new Float64Array( data.D );
+ E = new Float64Array( data.E );
+ WORK = new Float64Array( data.WORK );
+
+ expectedD = new Float64Array( data.D_out );
+ expectedE = new Float64Array( data.E_out );
+ expectedWORK = new Float64Array( data.WORK_out );
+
+ info = dlasq1( data.N, D, data.strideD, data.offsetD, E, data.strideE, data.offsetE, WORK, data.strideWORK, data.offsetWORK );
+
+ t.strictEqual( info, 0, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( D, expectedD, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( E, expectedE, 1 ), true, 'returns expected value' );
+ t.strictEqual( isAlmostEqual( WORK, expectedWORK, 1 ), true, 'returns expected value' );
+ t.end();
+});