From 4d899ac502e20977b2f31ceb844296899bc46579 Mon Sep 17 00:00:00 2001 From: Prajjwal Bajpai Date: Mon, 6 Jul 2026 23:20:16 +0530 Subject: [PATCH] feat: add `lapack/base/dlasq2` --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown_pkg_readmes status: na - task: lint_markdown_docs status: na - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../@stdlib/lapack/base/dlasq2/lib/base.js | 514 ++++++++++++++++++ .../@stdlib/lapack/base/dlasq2/lib/dlasq3.js | 368 +++++++++++++ .../@stdlib/lapack/base/dlasq2/lib/dlasq4.js | 401 ++++++++++++++ .../@stdlib/lapack/base/dlasq2/lib/dlasq5.js | 331 +++++++++++ .../@stdlib/lapack/base/dlasq2/lib/dlasq6.js | 220 ++++++++ .../@stdlib/lapack/base/dlasq2/lib/dlasrt.js | 250 +++++++++ 6 files changed, 2084 insertions(+) create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/base.js create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq3.js create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq4.js create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq5.js create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq6.js create mode 100644 lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasrt.js diff --git a/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/base.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/base.js new file mode 100644 index 000000000000..2908773f5e26 --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/base.js @@ -0,0 +1,514 @@ +/** +* @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 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 +*/ +function dlasq2( N, Z, strideZ, offsetZ ) { + var deemin; + var oldemn; + var safmin; + 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( 15 ); + + 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 ) ) + 7; + 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 * i0 ) + 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 ] = dmin; + out[ 1 ] = sigma; + out[ 2 ] = desig; + out[ 3 ] = qmax; + out[ 4 ] = nfail; + out[ 5 ] = iter; + out[ 6 ] = ndiv; + out[ 7 ] = ttype; + out[ 8 ] = dmin1; + out[ 9 ] = dmin2; + out[ 10 ] = dn; + out[ 11 ] = dn1; + out[ 12 ] = dn2; + out[ 13 ] = g; + out[ 14 ] = tau; + + // Call dlasq3 + dlasq3( i0, n0, Z, strideZ, offsetZ, pp, IEEE, out, 1, 0 ); + + dmin = out[ 0 ]; + sigma = out[ 1 ]; + desig = out[ 2 ]; + qmax = out[ 3 ]; + nfail = out[ 4 ]; + iter = out[ 5 ]; + ndiv = out[ 6 ]; + ttype = out[ 7 ]; + dmin1 = out[ 8 ]; + dmin2 = out[ 9 ]; + dn = out[ 10 ]; + dn1 = out[ 11 ]; + dn2 = out[ 12 ]; + g = out[ 13 ]; + tau = out[ 14 ]; + + 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 * n0 ) * strideZ ) ]; + emin = Z[ offsetZ + ( ( ( 4 * n0 ) + 2 ) * strideZ ) ]; + oldemn = Z[ offsetZ + ( ( ( 4 * n0 ) + 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/dlasq2/lib/dlasq3.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq3.js new file mode 100644 index 000000000000..6298fd895f5b --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq3.js @@ -0,0 +1,368 @@ +/** +* @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 {integer} N0 - last index +* @param {Float64Array} Z - qd array +* @param {integer} strideZ - stride length for `Z` +* @param {NonNegativeInteger} offsetZ - starting index of `Z` +* @param {integer} PP - ping-pong flag (0 or 1) +* @param {boolean} IEEE - IEEE arithmetic flag +* @param {Float64Array} out - output array containing `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( 15 ); +* 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, 3, Z, 1, 0, 0, true, out, 1, 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 ] +* // out => [ ~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 ] +*/ +function dlasq3( I0, N0, Z, strideZ, offsetZ, PP, 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 nn; + var g; + var s; + var t; + + // Read input values from output array + idx = offsetOut; + 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 ) { + 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 ] = 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/dlasq2/lib/dlasq4.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq4.js new file mode 100644 index 000000000000..53117a0a6373 --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq4.js @@ -0,0 +1,401 @@ +/** +* @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' ); + + +// 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 - 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 CNST1; + var CNST2; + var CNST3; + var ttype; + var idx2; + var gap1; + var gap2; + var gam; + var idx; + 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 ]; + + CNST1 = 0.563; + CNST2 = 1.010; + CNST3 = 1.050; + + // 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 ) ); + + 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; + for ( idx = ( np - offsetZ ) / strideZ; idx >= ( 4*I0 ) + PP + 2; idx -= 4 ) { + if ( b2 === 0 ) { + break; + } + b1 = b2; + if ( Z[ i4 ] > Z[ i4 - ( 2*strideZ ) ] ) { + idx2 = offsetOut + strideOut; + out[ idx2 ] = ttype; + idx2 += strideOut; + out[ idx2 ] = G; + return out; + } + b2 *= Z[ i4 ] / Z[ i4 - ( 2*strideZ ) ]; + a2 += b2; + if ( ( 100*max( b2, b1 ) ) < a2 || CNST1 < a2 ) { + break; + } + i4 -= 4*strideZ; + } + 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 - ( strideZ*17 ); + for ( idx = ( ( nn - offsetZ ) / strideZ ) - 17; idx >= ( 4*I0 ) + PP + 2; idx -= 4 ) { + if ( b2 === 0 ) { + break; + } + b1 = b2; + if ( Z[ i4 ] > Z[ i4 - ( 2*strideZ ) ] ) { + idx2 = offsetOut + strideOut; + out[ idx2 ] = ttype; + idx2 += strideOut; + out[ idx2 ] = G; + return out; + } + b2 *= ( Z[ i4 ] / Z[ i4 - ( 2*strideZ ) ] ); + a2 += b2; + if ( ( 100*max( b2, b1 ) ) < a2 || CNST1 < a2 ) { + break; + } + i4 -= 4*strideZ; + } + 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 ) ); + for ( idx = 2; idx >= 0; idx-- ) { + a2 = b1; + if ( Z[ i4 ] > Z[ i4 - ( 2*strideZ ) ] ) { + idx2 = offsetOut + strideOut; + out[ idx2 ] = ttype; + idx2 += strideOut; + out[ idx2 ] = G; + return out; + } + b1 *= Z[ i4 ] / Z[ i4 - ( 2*strideZ ) ]; + b2 += b1; + + if ( ( 100*max( b1, a2 ) ) < b2 ) { + break; + } + i4 -= 4*strideZ; + } + } + + 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 ) ); + for ( idx = 1; idx >= 0; idx-- ) { + if ( Z[ i4 ] > Z[ i4 - ( 2*strideZ ) ] ) { + idx2 = offsetOut + strideOut; + out[ idx2 ] = ttype; + idx2 += strideOut; + out[ idx2 ] = G; + return out; + } + b1 *= ( Z[ i4 ] / Z[ i4 - ( 2*strideZ ) ] ); + b2 += b1; + if ( ( 100*b1 ) < b2 ) { + break; + } + i4 -= 4*strideZ; + } + } + + 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/dlasq2/lib/dlasq5.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq5.js new file mode 100644 index 000000000000..8c938fe35377 --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/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/dlasq2/lib/dlasq6.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq6.js new file mode 100644 index 000000000000..8c925518502d --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasq6.js @@ -0,0 +1,220 @@ +/** +* @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 iz1; + var iz2; + var iz3; + var iz4; + 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 ) { + iz1 = j4 + strideZ; + iz2 = j4 - strideZ; + iz3 = j4 - ( 2*strideZ ); + + Z[ iz3 ] = d + Z[ iz2 ]; + if ( Z[ iz3 ] === 0 ) { + Z[ j4 ] = 0; + d = Z[ iz1 ]; + dmin = d; + emin = 0; + } else if ( ( SFMIN*Z[ iz1 ] < Z[ iz3 ] ) && + ( SFMIN*Z[ iz3 ] < Z[ iz1 ] ) ) { + temp = Z[ iz1 ] / Z[ iz3 ]; + Z[ j4 ] = Z[ iz2 ]*temp; + d *= temp; + } else { + Z[ j4 ] = Z[ iz1 ]*( Z[ iz2 ] / Z[ iz3 ] ); + d = Z[ iz1 ]*( d / Z[ iz3 ] ); + } + dmin = min( dmin, d ); + emin = min( emin, Z[ j4 ] ); + j4 += 4*strideZ; + } + } else { + for ( idx = ( 4*I0 ) + 3; idx <= ( 4*N0 ) - 9; idx += 4 ) { + iz1 = j4 + ( 2*strideZ ); + iz2 = j4 - strideZ; + iz3 = j4 - ( 2*strideZ ); + iz4 = j4 - ( 3*strideZ ); + + Z[ iz4 ] = d + Z[ j4 ]; + if ( Z[ iz4 ] === 0 ) { + Z[ iz2 ] = 0; + d = Z[ iz1 ]; + dmin = d; + emin = 0; + } else if ( ( SFMIN*Z[ iz1 ] < Z[ iz4 ] ) && + ( SFMIN*Z[ iz4 ] < Z[ iz1 ] ) ) { + temp = Z[ iz1 ] / Z[ iz4 ]; + Z[ j4 - strideZ ] = Z[ j4 ]*temp; + d *= temp; + } else { + Z[ iz2 ] = Z[ iz1 ]*( Z[ j4 ] / Z[ iz4 ] ); + d = Z[ iz1 ]*( d / Z[ iz4 ] ); + } + dmin = min( dmin, d ); + emin = min( emin, Z[ iz2 ] ); + j4 += 4*strideZ; + } + } + + // Unroll last two steps. + dnm2 = d; + dmin2 = dmin; + j4 = offsetZ + ( strideZ*( ( 4*N0 ) - 5 - PP ) ); + j4p2 = j4 + ( strideZ*( ( 2*PP ) - 1 ) ); + + iz1 = j4 - ( 2*strideZ ); + iz2 = j4p2 + ( 2*strideZ ); + + Z[ j4 - ( 2*strideZ ) ] = dnm2 + Z[ j4p2 ]; + if ( Z[ iz1 ] === 0 ) { + Z[ j4 ] = 0; + dnm1 = Z[ iz2 ]; + dmin = dnm1; + emin = 0; + } else if ( ( SFMIN*Z[ iz2 ] < Z[ iz1 ] ) && + ( SFMIN*Z[ iz1 ] < Z[ iz2 ] ) ) { + temp = Z[ iz2 ] / Z[ iz1 ]; + Z[ j4 ] = Z[ j4p2 ]*temp; + dnm1 = dnm2*temp; + } else { + Z[ j4 ] = Z[ iz2 ]*( Z[ j4p2 ] / Z[ iz1 ] ); + dnm1 = Z[ iz2 ]*( dnm2 / Z[ iz1 ] ); + } + dmin = min( dmin, dnm1 ); + + dmin1 = dmin; + j4 += 4*strideZ; + j4p2 = j4 + ( strideZ*( ( 2*PP ) - 1 ) ); + + iz1 = j4 - ( 2*strideZ ); + iz2 = j4p2 + ( 2*strideZ ); + + Z[ iz1 ] = dnm1 + Z[ j4p2 ]; + if ( Z[ iz1 ] === 0 ) { + Z[ j4 ] = 0; + dn = Z[ iz2 ]; + dmin = dn; + emin = 0; + } else if ( ( SFMIN*Z[ iz2 ] < Z[ iz1 ] ) && + ( SFMIN*Z[ iz1 ] < Z[ iz2 ] ) ) { + temp = Z[ iz2 ] / Z[ iz1 ]; + Z[ j4 ] = Z[ j4p2 ]*temp; + dn = dnm1*temp; + } else { + Z[ j4 ] = Z[ iz2 ]*( Z[ j4p2 ] / Z[ iz1 ] ); + dn = Z[ iz2 ]*( dnm1 / Z[ iz1 ] ); + } + 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/dlasq2/lib/dlasrt.js b/lib/node_modules/@stdlib/lapack/base/dlasq2/lib/dlasrt.js new file mode 100644 index 000000000000..01ad9afb80a9 --- /dev/null +++ b/lib/node_modules/@stdlib/lapack/base/dlasq2/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;