1616* limitations under the License.
1717*/
1818
19- /* eslint-disable max-len, max-params, max-statements, max-lines-per-function */
19+ /* eslint-disable max-len, max-params, max-statements */
2020
2121'use strict' ;
2222
@@ -28,6 +28,7 @@ var dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray;
2828var idamax = require ( '@stdlib/blas/base/idamax' ) . ndarray ;
2929var abs = require ( '@stdlib/math/base/special/abs' ) ;
3030var isnan = require ( '@stdlib/assert/is-nan' ) ;
31+ var dfill = require ( '@stdlib/blas/ext/base/dfill' ) . ndarray ;
3132var max = require ( '@stdlib/math/base/special/maxn' ) ;
3233var min = require ( '@stdlib/math/base/special/minn' ) ;
3334var dscal = require ( '@stdlib/blas/base/dscal' ) . ndarray ;
@@ -37,6 +38,10 @@ var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
3738
3839var sclfac = 2.0 ;
3940var factor = 0.95 ;
41+ var sfmin1 = dlamch ( 'S' ) / dlamch ( 'P' ) ;
42+ var sfmax1 = 1.0 / sfmin1 ;
43+ var sfmin2 = sfmin1 * sclfac ;
44+ var sfmax2 = 1.0 / sfmin2 ;
4045
4146
4247// MAIN //
@@ -74,7 +79,7 @@ var factor = 0.95;
7479* @param {Float64Array } scale - array containing permutation and scaling information
7580* @param {integer } strideScale - stride of `scale`
7681* @param {NonNegativeInteger } offsetScale - starting index for `scale`
77- * @throws {RangeError } should not return NaN
82+ * @throws {RangeError } the input matrix cannot be balanced. Encountered a NaN during computation.
7883* @returns {integer } status code
7984*
8085* @example
@@ -93,10 +98,6 @@ var factor = 0.95;
9398function dgebal ( job , N , A , strideA1 , strideA2 , offsetA , out , strideOut , offsetOut , scale , strideScale , offsetScale ) {
9499 var canSwap ;
95100 var noconv ;
96- var sfmin1 ;
97- var sfmin2 ;
98- var sfmax1 ;
99- var sfmax2 ;
100101 var ica ;
101102 var ira ;
102103 var ia1 ;
@@ -124,11 +125,7 @@ function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetO
124125 }
125126
126127 if ( job === 'none' ) {
127- is = offsetScale ;
128- for ( i = 0 ; i < N ; i ++ ) {
129- scale [ is ] = 1.0 ;
130- is += strideScale ;
131- }
128+ dfill ( N , 1.0 , scale , strideScale , offsetScale ) ;
132129
133130 out [ offsetOut ] = 0.0 ; // ilo
134131 out [ offsetOut + strideOut ] = N - 1 ; // ihi
@@ -227,11 +224,7 @@ function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetO
227224 }
228225
229226 // Initialize scale for non-permuted submatrix
230- is = offsetScale + ( k * strideScale ) ;
231- for ( i = k ; i <= l ; i ++ ) {
232- scale [ is ] = 1.0 ;
233- is += strideScale ;
234- }
227+ dfill ( N - k , 1.0 , scale , strideScale , offsetScale + ( k * strideScale ) ) ;
235228
236229 if ( job === 'permute' ) {
237230 out [ offsetOut ] = k ; // ilo
@@ -240,11 +233,6 @@ function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetO
240233 }
241234
242235 // Balance the submatrix in rows K to L, iterative loop for norm reduction (job = 'B')
243- sfmin1 = dlamch ( 'S' ) / dlamch ( 'P' ) ;
244- sfmax1 = 1.0 / sfmin1 ;
245- sfmin2 = sfmin1 * sclfac ;
246- sfmax2 = 1.0 / sfmin2 ;
247-
248236 noconv = true ;
249237 while ( noconv ) {
250238 is = offsetScale + ( k * strideScale ) ; // Follows scale
@@ -270,7 +258,7 @@ function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetO
270258 }
271259
272260 if ( isnan ( c ) || isnan ( r ) || isnan ( ca ) || isnan ( ra ) ) {
273- throw new RangeError ( 'should not return NaN' ) ;
261+ throw new RangeError ( 'invalid argument. The input matrix cannot be balanced. Encountered a NaN during computation. ' ) ;
274262 }
275263
276264 g = r / sclfac ;
@@ -306,24 +294,20 @@ function dgebal( job, N, A, strideA1, strideA2, offsetA, out, strideOut, offsetO
306294 continue ;
307295 }
308296
309- if ( f < 1.0 && scale [ is ] < 1.0 ) {
310- if ( f * scale [ is ] <= sfmin1 ) {
311- ia1 += strideA2 ;
312- ia2 += strideA1 ;
313- is += strideScale ;
314- ia3 += strideA2 ;
315- continue ;
316- }
297+ if ( f < 1.0 && scale [ is ] < 1.0 && f * scale [ is ] <= sfmin1 ) {
298+ ia1 += strideA2 ;
299+ ia2 += strideA1 ;
300+ is += strideScale ;
301+ ia3 += strideA2 ;
302+ continue ;
317303 }
318304
319- if ( f > 1.0 && scale [ is ] > 1.0 ) {
320- if ( scale [ is ] >= sfmax1 / f ) {
321- ia1 += strideA2 ;
322- ia2 += strideA1 ;
323- is += strideScale ;
324- ia3 += strideA2 ;
325- continue ;
326- }
305+ if ( f > 1.0 && scale [ is ] > 1.0 && scale [ is ] >= sfmax1 / f ) {
306+ ia1 += strideA2 ;
307+ ia2 += strideA1 ;
308+ is += strideScale ;
309+ ia3 += strideA2 ;
310+ continue ;
327311 }
328312
329313 g = 1.0 / f ;
0 commit comments