Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 | 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 10x 10x 10x 10x 10x 10x 10x 10x 10x 10x 10x 2x 2x 10x 2x 2x 2x 10x 2x 2x 2x 2x 2x 2x 2x 2x 4x 4x 4x 4x 4x 10x 8x 8x 8x 8x 8x 4x 4x 4x 10x 2x 2x 2x 2x 2x 10x 6x 6x 6x 2x 2x 2x 2x 2x 2x 2x 2x 10x 10x 10x 10x 2x 2x 2x 2x 2x 2x 2x 2x 2x 6x 6x 6x 6x 2x 10x 2x 10x 3x 3x 3x 3x 3x | /**
* @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 => <Float64Array>[ ~139.377, ~56.064, ~15.997 ]
* // E => <Float64Array>[ 90, 40, 0 ]
* // WORK => <Float64Array>[ ~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;
|