|
| 1 | +/** |
| 2 | +* @license Apache-2.0 |
| 3 | +* |
| 4 | +* Copyright (c) 2025 The Stdlib Authors. |
| 5 | +* |
| 6 | +* Licensed under the Apache License, Version 2.0 (the "License"); |
| 7 | +* you may not use this file except in compliance with the License. |
| 8 | +* You may obtain a copy of the License at |
| 9 | +* |
| 10 | +* http://www.apache.org/licenses/LICENSE-2.0 |
| 11 | +* |
| 12 | +* Unless required by applicable law or agreed to in writing, software |
| 13 | +* distributed under the License is distributed on an "AS IS" BASIS, |
| 14 | +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 15 | +* See the License for the specific language governing permissions and |
| 16 | +* limitations under the License. |
| 17 | +*/ |
| 18 | + |
| 19 | +/* eslint-disable max-len, max-params, max-statements */ |
| 20 | + |
| 21 | +'use strict'; |
| 22 | + |
| 23 | +// MODULES // |
| 24 | + |
| 25 | +var dswap = require( '@stdlib/blas/base/dswap' ).ndarray; |
| 26 | +var dlamch = require( '@stdlib/lapack/base/dlamch' ); |
| 27 | +var dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray; |
| 28 | +var idamax = require( '@stdlib/blas/base/idamax' ).ndarray; |
| 29 | +var abs = require( '@stdlib/math/base/special/abs' ); |
| 30 | +var isnan = require( '@stdlib/assert/is-nan' ); |
| 31 | +var max = require( '@stdlib/math/base/special/maxn' ); |
| 32 | +var min = require( '@stdlib/math/base/special/minn' ); |
| 33 | +var dscal = require( '@stdlib/blas/base/dscal' ).ndarray; |
| 34 | + |
| 35 | + |
| 36 | +// MAIN // |
| 37 | + |
| 38 | +/** |
| 39 | +* Balances a general real matrix `A`. |
| 40 | +* |
| 41 | +* ## Notes |
| 42 | +* |
| 43 | +* The job parameter can be one of the following: |
| 44 | +* |
| 45 | +* - 'N': none, return immediately |
| 46 | +* - 'P': permute only |
| 47 | +* - 'S': scale only |
| 48 | +* - 'B': both permute and scale |
| 49 | +* - The matrix `A` is overwritten by the balanced matrix. |
| 50 | +* |
| 51 | +* @private |
| 52 | +* @param {string} job - indicates the operations to be performed |
| 53 | +* @param {NonNegativeInteger} N - number of rows/columns in matrix `A` |
| 54 | +* @param {Float64Array} A - input matrix to be balanced |
| 55 | +* @param {integer} strideA1 - stride of the first dimension of `A` |
| 56 | +* @param {integer} strideA2 - stride of the second dimension of `A` |
| 57 | +* @param {NonNegativeInteger} offsetA - starting index for `A` |
| 58 | +* @param {Float64Array} out - stores the first and last row/column of the balanced submatrix |
| 59 | +* @param {integer} strideOut - stride of `out` |
| 60 | +* @param {NonNegativeInteger} offsetOut - starting index for `out` |
| 61 | +* @param {Float64Array} scale - array containing permutation and scaling information |
| 62 | +* @param {integer} strideScale - stride of `scale` |
| 63 | +* @param {NonNegativeInteger} offsetScale - starting index for `scale` |
| 64 | +* @returns {integer} status code |
| 65 | +* |
| 66 | +* @example |
| 67 | +* var Float64Array = require( '@stdlib/array/float64' ); |
| 68 | +* |
| 69 | +* var A = new Float64Array( [ 1.0, 100.0, 0.0, 2.0, 200.0, 0.0, 0.0, 0.0, 3.0 ] ); |
| 70 | +* var out = new Float64Array( 2 ); |
| 71 | +* var scale = new Float64Array( 3 ); |
| 72 | +* |
| 73 | +* dgebal( 'B', 3, A, 3, 1, 0, out, 1, 0, scale, 1, 0 ); |
| 74 | +* // A => <Float64Array>[ 1, 12.5, 0, 16, 200, 0, 0, 0, 3 ] |
| 75 | +* // out => <Float64Array>[ 0, 1 ] |
| 76 | +* // scale => <Float64Array>[ 8, 1, 2 ] |
| 77 | +*/ |
| 78 | +function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetOut, scale, strideScale, offsetScale ) { |
| 79 | + var canSwap; |
| 80 | + var noconv; |
| 81 | + var sfmin1; |
| 82 | + var sfmin2; |
| 83 | + var sfmax1; |
| 84 | + var sfmax2; |
| 85 | + var sclfac; |
| 86 | + var factor; |
| 87 | + var ica; |
| 88 | + var ira; |
| 89 | + var ca; |
| 90 | + var ra; |
| 91 | + var is; |
| 92 | + var c; |
| 93 | + var r; |
| 94 | + var k; |
| 95 | + var l; |
| 96 | + var i; |
| 97 | + var j; |
| 98 | + var g; |
| 99 | + var f; |
| 100 | + var s; |
| 101 | + |
| 102 | + sclfac = 2.0; |
| 103 | + factor = 0.95; |
| 104 | + |
| 105 | + // Quick return if possible |
| 106 | + if ( N === 0 ) { |
| 107 | + out[ offsetOut ] = 0.0; // ilo |
| 108 | + out[ offsetOut + strideOut ] = -1.0; // ihi (invalid) |
| 109 | + return 0; |
| 110 | + } |
| 111 | + |
| 112 | + if ( job === 'N' ) { |
| 113 | + is = offsetScale; |
| 114 | + for ( i = 0; i < N; i++ ) { |
| 115 | + scale[ is ] = 1.0; |
| 116 | + is += strideScale; |
| 117 | + } |
| 118 | + |
| 119 | + out[ offsetOut ] = 0.0; // ilo |
| 120 | + out[ offsetOut + strideOut ] = N - 1; // ihi |
| 121 | + return 0; |
| 122 | + } |
| 123 | + |
| 124 | + // Permutation to isolate eigenvalues if possible |
| 125 | + k = 0; |
| 126 | + l = N - 1; |
| 127 | + |
| 128 | + if ( job !== 'S' ) { |
| 129 | + // Row and column exchange |
| 130 | + noconv = true; |
| 131 | + while ( noconv ) { |
| 132 | + // Search for rows isolating an eigenvalue and push them down |
| 133 | + noconv = false; |
| 134 | + for ( i = l; i >= 0; i-- ) { |
| 135 | + canSwap = true; |
| 136 | + for ( j = 0; j <= l; j++ ) { |
| 137 | + if ( i !== j && A[ offsetA + (i*strideA1) + (j*strideA2) ] !== 0.0 ) { |
| 138 | + canSwap = false; |
| 139 | + break; |
| 140 | + } |
| 141 | + } |
| 142 | + |
| 143 | + if ( canSwap ) { |
| 144 | + scale[ offsetScale + (l*strideScale) ] = i; |
| 145 | + if ( i !== l ) { |
| 146 | + dswap( l+1, A, strideA1, offsetA + (i*strideA2), A, strideA1, offsetA + (l*strideA2) ); |
| 147 | + dswap( N - k, A, strideA2, offsetA + (i*strideA1) + (k*strideA2), A, strideA2, offsetA + (l*strideA1) + (k*strideA2) ); |
| 148 | + } |
| 149 | + noconv = true; |
| 150 | + |
| 151 | + if ( l === 0.0 ) { |
| 152 | + out[ offsetOut ] = 0.0; // ilo |
| 153 | + out[ offsetOut + strideOut ] = 0.0; // ihi |
| 154 | + return 0; |
| 155 | + } |
| 156 | + l -= 1; |
| 157 | + } |
| 158 | + } |
| 159 | + } |
| 160 | + |
| 161 | + noconv = true; |
| 162 | + while ( noconv ) { |
| 163 | + // Search for columns isolating an eigenvalue and push them left |
| 164 | + noconv = false; |
| 165 | + for ( j = k; j <= l; j++ ) { |
| 166 | + canSwap = true; |
| 167 | + for ( i = k; i <= l; i++ ) { |
| 168 | + if ( i !== j && A[ offsetA + (i*strideA1) + (j*strideA2) ] !== 0.0 ) { |
| 169 | + canSwap = false; |
| 170 | + break; |
| 171 | + } |
| 172 | + } |
| 173 | + |
| 174 | + if ( canSwap ) { |
| 175 | + scale[ offsetScale + (k*strideScale) ] = j; |
| 176 | + if ( j !== k ) { |
| 177 | + dswap( l+1, A, strideA1, offsetA + (j*strideA2), A, strideA1, offsetA + (k*strideA2) ); |
| 178 | + dswap( N-k, A, strideA2, offsetA + (j*strideA1), A, strideA2, offsetA + (k*strideA1) ); |
| 179 | + } |
| 180 | + noconv = true; |
| 181 | + k += 1; |
| 182 | + } |
| 183 | + } |
| 184 | + } |
| 185 | + } |
| 186 | + |
| 187 | + // Initialize `scale` for non-permuted submatrix |
| 188 | + is = offsetScale; |
| 189 | + for ( i = k; i <= l; i++ ) { |
| 190 | + scale[ is ] = 1.0; |
| 191 | + is += strideScale; |
| 192 | + } |
| 193 | + |
| 194 | + if ( job === 'P' ) { |
| 195 | + out[ offsetOut ] = k; // ilo |
| 196 | + out[ offsetOut + strideOut ] = l; // ihi |
| 197 | + return 0; |
| 198 | + } |
| 199 | + |
| 200 | + // Balance the submatrix in rows K to L, iterative loop for norm reduction |
| 201 | + sfmin1 = dlamch( 'S' ) / dlamch( 'P' ); |
| 202 | + sfmax1 = 1.0 / sfmin1; |
| 203 | + sfmin2 = sfmin1 * sclfac; |
| 204 | + sfmax2 = 1.0 / sfmin2; |
| 205 | + |
| 206 | + noconv = true; |
| 207 | + while ( noconv ) { |
| 208 | + noconv = false; |
| 209 | + for ( i = k; i <= l; i++ ) { |
| 210 | + c = dnrm2( l-k+1, A, strideA1, offsetA + (k*strideA1) + (i*strideA2) ); |
| 211 | + r = dnrm2( l-k+1, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) ); |
| 212 | + ica = idamax( l+1, A, strideA1, offsetA + (i*strideA2) ); |
| 213 | + ca = abs( A[ offsetA + (ica*strideA1) + (i*strideA2) ] ); |
| 214 | + ira = idamax( N-k+1, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) ); |
| 215 | + ra = abs( A[ offsetA + (i*strideA1) + ((ira+k)*strideA2) ] ); |
| 216 | + |
| 217 | + if ( c === 0.0 || r === 0.0 ) { |
| 218 | + continue; |
| 219 | + } |
| 220 | + |
| 221 | + if ( isnan( c ) || isnan( r ) || isnan( ca ) || isnan( ra ) ) { |
| 222 | + return -3; |
| 223 | + } |
| 224 | + |
| 225 | + g = r / sclfac; |
| 226 | + f = 1.0; |
| 227 | + s = c + r; |
| 228 | + |
| 229 | + while ( c < g && max( f, c, ca ) < sfmax2 && min( r, g, ra ) > sfmin2 ) { |
| 230 | + f *= sclfac; |
| 231 | + c *= sclfac; |
| 232 | + ca *= sclfac; |
| 233 | + r /= sclfac; |
| 234 | + g /= sclfac; |
| 235 | + ra /= sclfac; |
| 236 | + } |
| 237 | + |
| 238 | + g = c / sclfac; |
| 239 | + |
| 240 | + while ( g >= r && max( r, ra ) < sfmax2 && min( f, c, g, ca ) > sfmin2 ) { |
| 241 | + f /= sclfac; |
| 242 | + c /= sclfac; |
| 243 | + g /= sclfac; |
| 244 | + ca /= sclfac; |
| 245 | + r *= sclfac; |
| 246 | + ra *= sclfac; |
| 247 | + } |
| 248 | + |
| 249 | + // Now balance |
| 250 | + if ( ( c + r ) >= factor * s ) { |
| 251 | + continue; |
| 252 | + } |
| 253 | + |
| 254 | + if ( f < 1.0 && scale[ offsetScale + (i*strideScale) ] < 1.0 ) { |
| 255 | + if ( f * scale[ offsetScale + (i*strideScale) ] <= sfmin1 ) { |
| 256 | + continue; |
| 257 | + } |
| 258 | + } |
| 259 | + |
| 260 | + if ( f > 1.0 && scale[ offsetScale + (i*strideScale) ] > 1.0 ) { |
| 261 | + if ( scale[ offsetScale + (i*strideScale) ] >= sfmax1 / f ) { |
| 262 | + continue; |
| 263 | + } |
| 264 | + } |
| 265 | + |
| 266 | + g = 1.0 / f; |
| 267 | + scale[ offsetScale + (i*strideScale) ] *= f; |
| 268 | + noconv = true; |
| 269 | + |
| 270 | + dscal( N-k, g, A, strideA2, offsetA + (i*strideA1) + (k*strideA2) ); |
| 271 | + dscal( l+1, f, A, strideA1, offsetA + (i*strideA2) ); |
| 272 | + } |
| 273 | + } |
| 274 | + |
| 275 | + out[ offsetOut ] = k; // ilo |
| 276 | + out[ offsetOut + strideOut ] = l; // ihi |
| 277 | + return 0; |
| 278 | +} |
| 279 | + |
| 280 | + |
| 281 | +// EXPORTS // |
| 282 | + |
| 283 | +module.exports = dgebal; |
0 commit comments