Skip to content

Commit 0a7dec1

Browse files
angschangsch
committed
Avoid NaN generation in LATRS
Fixes #714
1 parent 7ce9bf5 commit 0a7dec1

File tree

4 files changed

+262
-12
lines changed

4 files changed

+262
-12
lines changed

‎SRC/clatrs.f‎

Lines changed: 68 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -357,8 +357,74 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
357357
IF( TMAX.LE.BIGNUM*HALF ) THEN
358358
TSCAL = ONE
359359
ELSE
360-
TSCAL = HALF / ( SMLNUM*TMAX )
361-
CALL SSCAL( N, TSCAL, CNORM, 1 )
360+
*
361+
* Avoid NaN generation if entries in CNORM exceed the
362+
* overflow threshold
363+
*
364+
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
365+
* Case 1: All entries in CNORM are valid floating-point numbers
366+
TSCAL = HALF / ( SMLNUM*TMAX )
367+
CALL SSCAL( N, TSCAL, CNORM, 1 )
368+
ELSE
369+
* Case 2: At least one column norm of A cannot be
370+
* represented as a floating-point number. Find the
371+
* maximum offdiagonal absolute value
372+
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
373+
* not +/- Infinity, use this value as TSCAL.
374+
TMAX = ZERO
375+
IF( UPPER ) THEN
376+
*
377+
* A is upper triangular.
378+
*
379+
DO J = 2, N
380+
DO I = 1, J - 1
381+
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
382+
$ ABS( AIMAG(A ( I, J ) ) ) )
383+
END DO
384+
END DO
385+
ELSE
386+
*
387+
* A is lower triangular.
388+
*
389+
DO J = 1, N - 1
390+
DO I = J + 1, N
391+
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
392+
$ ABS( AIMAG(A ( I, J ) ) ) )
393+
END DO
394+
END DO
395+
END IF
396+
*
397+
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
398+
TSCAL = ONE / ( SMLNUM*TMAX )
399+
DO J = 1, N
400+
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
401+
CNORM( J ) = CNORM( J )*TSCAL
402+
ELSE
403+
* Recompute the 1-norm of each column without
404+
* introducing Infinity in the summation.
405+
TSCAL = TWO * TSCAL
406+
CNORM( J ) = ZERO
407+
IF( UPPER ) THEN
408+
DO I = 1, J - 1
409+
CNORM( J ) = CNORM( J ) +
410+
$ TSCAL * CABS2( A( I, J ) )
411+
END DO
412+
ELSE
413+
DO I = J + 1, N
414+
CNORM( J ) = CNORM( J ) +
415+
$ TSCAL * CABS2( A( I, J ) )
416+
END DO
417+
END IF
418+
TSCAL = TSCAL * HALF
419+
END IF
420+
END DO
421+
ELSE
422+
* At least one entry of A is not a valid floating-point
423+
* entry. Rely on TRSV to propagate Inf and NaN.
424+
CALL CTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
425+
RETURN
426+
END IF
427+
END IF
362428
END IF
363429
*
364430
* Compute a bound on the computed solution vector to see if the

‎SRC/dlatrs.f‎

Lines changed: 63 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -264,8 +264,8 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
264264
* .. External Functions ..
265265
LOGICAL LSAME
266266
INTEGER IDAMAX
267-
DOUBLE PRECISION DASUM, DDOT, DLAMCH
268-
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
267+
DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
268+
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
269269
* ..
270270
* .. External Subroutines ..
271271
EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
@@ -343,8 +343,67 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
343343
IF( TMAX.LE.BIGNUM ) THEN
344344
TSCAL = ONE
345345
ELSE
346-
TSCAL = ONE / ( SMLNUM*TMAX )
347-
CALL DSCAL( N, TSCAL, CNORM, 1 )
346+
*
347+
* Avoid NaN generation if entries in CNORM exceed the
348+
* overflow threshold
349+
*
350+
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
351+
* Case 1: All entries in CNORM are valid floating-point numbers
352+
TSCAL = ONE / ( SMLNUM*TMAX )
353+
CALL DSCAL( N, TSCAL, CNORM, 1 )
354+
ELSE
355+
* Case 2: At least one column norm of A cannot be represented
356+
* as floating-point number. Find the offdiagonal entry A( I, J )
357+
* with the largest absolute value. If this entry is not +/- Infinity,
358+
* use this value as TSCAL.
359+
TMAX = ZERO
360+
IF( UPPER ) THEN
361+
*
362+
* A is upper triangular.
363+
*
364+
DO J = 2, N
365+
TMAX = MAX( DLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
366+
$ TMAX )
367+
END DO
368+
ELSE
369+
*
370+
* A is lower triangular.
371+
*
372+
DO J = 1, N - 1
373+
TMAX = MAX( DLANGE( 'M', N-J, 1, A( J+1, J ), 1,
374+
$ SUMJ ), TMAX )
375+
END DO
376+
END IF
377+
*
378+
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
379+
TSCAL = ONE / ( SMLNUM*TMAX )
380+
DO J = 1, N
381+
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
382+
CNORM( J ) = CNORM( J )*TSCAL
383+
ELSE
384+
* Recompute the 1-norm without introducing Infinity
385+
* in the summation
386+
CNORM( J ) = ZERO
387+
IF( UPPER ) THEN
388+
DO I = 1, J - 1
389+
CNORM( J ) = CNORM( J ) +
390+
$ TSCAL * ABS( A( I, J ) )
391+
END DO
392+
ELSE
393+
DO I = J + 1, N
394+
CNORM( J ) = CNORM( J ) +
395+
$ TSCAL * ABS( A( I, J ) )
396+
END DO
397+
END IF
398+
END IF
399+
END DO
400+
ELSE
401+
* At least one entry of A is not a valid floating-point entry.
402+
* Rely on TRSV to propagate Inf and NaN.
403+
CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
404+
RETURN
405+
END IF
406+
END IF
348407
END IF
349408
*
350409
* Compute a bound on the computed solution vector to see if the

‎SRC/slatrs.f‎

Lines changed: 63 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -264,8 +264,8 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
264264
* .. External Functions ..
265265
LOGICAL LSAME
266266
INTEGER ISAMAX
267-
REAL SASUM, SDOT, SLAMCH
268-
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH
267+
REAL SASUM, SDOT, SLAMCH, SLANGE
268+
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
269269
* ..
270270
* .. External Subroutines ..
271271
EXTERNAL SAXPY, SSCAL, STRSV, XERBLA
@@ -343,8 +343,67 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
343343
IF( TMAX.LE.BIGNUM ) THEN
344344
TSCAL = ONE
345345
ELSE
346-
TSCAL = ONE / ( SMLNUM*TMAX )
347-
CALL SSCAL( N, TSCAL, CNORM, 1 )
346+
*
347+
* Avoid NaN generation if entries in CNORM exceed the
348+
* overflow threshold
349+
*
350+
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
351+
* Case 1: All entries in CNORM are valid floating-point numbers
352+
TSCAL = ONE / ( SMLNUM*TMAX )
353+
CALL SSCAL( N, TSCAL, CNORM, 1 )
354+
ELSE
355+
* Case 2: At least one column norm of A cannot be represented
356+
* as floating-point number. Find the offdiagonal entry A( I, J )
357+
* with the largest absolute value. If this entry is not +/- Infinity,
358+
* use this value as TSCAL.
359+
TMAX = ZERO
360+
IF( UPPER ) THEN
361+
*
362+
* A is upper triangular.
363+
*
364+
DO J = 2, N
365+
TMAX = MAX( SLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
366+
$ TMAX )
367+
END DO
368+
ELSE
369+
*
370+
* A is lower triangular.
371+
*
372+
DO J = 1, N - 1
373+
TMAX = MAX( SLANGE( 'M', N-J, 1, A( J+1, J ), 1,
374+
$ SUMJ ), TMAX )
375+
END DO
376+
END IF
377+
*
378+
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
379+
TSCAL = ONE / ( SMLNUM*TMAX )
380+
DO J = 1, N
381+
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
382+
CNORM( J ) = CNORM( J )*TSCAL
383+
ELSE
384+
* Recompute the 1-norm without introducing Infinity
385+
* in the summation
386+
CNORM( J ) = ZERO
387+
IF( UPPER ) THEN
388+
DO I = 1, J - 1
389+
CNORM( J ) = CNORM( J ) +
390+
$ TSCAL * ABS( A( I, J ) )
391+
END DO
392+
ELSE
393+
DO I = J + 1, N
394+
CNORM( J ) = CNORM( J ) +
395+
$ TSCAL * ABS( A( I, J ) )
396+
END DO
397+
END IF
398+
END IF
399+
END DO
400+
ELSE
401+
* At least one entry of A is not a valid floating-point entry.
402+
* Rely on TRSV to propagate Inf and NaN.
403+
CALL STRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
404+
RETURN
405+
END IF
406+
END IF
348407
END IF
349408
*
350409
* Compute a bound on the computed solution vector to see if the

‎SRC/zlatrs.f‎

Lines changed: 68 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -357,8 +357,74 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
357357
IF( TMAX.LE.BIGNUM*HALF ) THEN
358358
TSCAL = ONE
359359
ELSE
360-
TSCAL = HALF / ( SMLNUM*TMAX )
361-
CALL DSCAL( N, TSCAL, CNORM, 1 )
360+
*
361+
* Avoid NaN generation if entries in CNORM exceed the
362+
* overflow threshold
363+
*
364+
IF ( TMAX.LE.DLAMCH('Overflow') ) THEN
365+
* Case 1: All entries in CNORM are valid floating-point numbers
366+
TSCAL = HALF / ( SMLNUM*TMAX )
367+
CALL DSCAL( N, TSCAL, CNORM, 1 )
368+
ELSE
369+
* Case 2: At least one column norm of A cannot be
370+
* represented as a floating-point number. Find the
371+
* maximum offdiagonal absolute value
372+
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
373+
* not +/- Infinity, use this value as TSCAL.
374+
TMAX = ZERO
375+
IF( UPPER ) THEN
376+
*
377+
* A is upper triangular.
378+
*
379+
DO J = 2, N
380+
DO I = 1, J - 1
381+
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
382+
$ ABS( DIMAG(A ( I, J ) ) ) )
383+
END DO
384+
END DO
385+
ELSE
386+
*
387+
* A is lower triangular.
388+
*
389+
DO J = 1, N - 1
390+
DO I = J + 1, N
391+
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
392+
$ ABS( DIMAG(A ( I, J ) ) ) )
393+
END DO
394+
END DO
395+
END IF
396+
*
397+
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
398+
TSCAL = ONE / ( SMLNUM*TMAX )
399+
DO J = 1, N
400+
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
401+
CNORM( J ) = CNORM( J )*TSCAL
402+
ELSE
403+
* Recompute the 1-norm of each column without
404+
* introducing Infinity in the summation.
405+
TSCAL = TWO * TSCAL
406+
CNORM( J ) = ZERO
407+
IF( UPPER ) THEN
408+
DO I = 1, J - 1
409+
CNORM( J ) = CNORM( J ) +
410+
$ TSCAL * CABS2( A( I, J ) )
411+
END DO
412+
ELSE
413+
DO I = J + 1, N
414+
CNORM( J ) = CNORM( J ) +
415+
$ TSCAL * CABS2( A( I, J ) )
416+
END DO
417+
END IF
418+
TSCAL = TSCAL * HALF
419+
END IF
420+
END DO
421+
ELSE
422+
* At least one entry of A is not a valid floating-point
423+
* entry. Rely on TRSV to propagate Inf and NaN.
424+
CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
425+
RETURN
426+
END IF
427+
END IF
362428
END IF
363429
*
364430
* Compute a bound on the computed solution vector to see if the

0 commit comments

Comments
 (0)