Skip to content

Commit

Permalink
Avoid NaN generation in LATRS
Browse files Browse the repository at this point in the history 
Fixes Reference-LAPACK#714
  • Loading branch information
@angsch
angsch committed Sep 25, 2022
1 parent 7ce9bf5 commit 0a7dec1
Show file tree
Hide file tree
Showing 4 changed files with 262 additions and 12 deletions.
70 changes: 68 additions & 2 deletions SRC/clatrs.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -357,8 +357,74 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
IF( TMAX.LE.BIGNUM*HALF ) THEN
TSCAL = ONE
ELSE
TSCAL = HALF / ( SMLNUM*TMAX )
CALL SSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = HALF / ( SMLNUM*TMAX )
CALL SSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be
* represented as a floating-point number. Find the
* maximum offdiagonal absolute value
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
* not +/- Infinity, use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
DO I = 1, J - 1
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
$ ABS( AIMAG(A ( I, J ) ) ) )
END DO
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
DO I = J + 1, N
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
$ ABS( AIMAG(A ( I, J ) ) ) )
END DO
END DO
END IF
*
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm of each column without
* introducing Infinity in the summation.
TSCAL = TWO * TSCAL
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * CABS2( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * CABS2( A( I, J ) )
END DO
END IF
TSCAL = TSCAL * HALF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point
* entry. Rely on TRSV to propagate Inf and NaN.
CALL CTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the
Expand Down
67 changes: 63 additions & 4 deletions SRC/dlatrs.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,8 +264,8 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DASUM, DDOT, DLAMCH
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
Expand Down Expand Up @@ -343,8 +343,67 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
IF( TMAX.LE.BIGNUM ) THEN
TSCAL = ONE
ELSE
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = ONE / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be represented
* as floating-point number. Find the offdiagonal entry A( I, J )
* with the largest absolute value. If this entry is not +/- Infinity,
* use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
TMAX = MAX( DLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
$ TMAX )
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
TMAX = MAX( DLANGE( 'M', N-J, 1, A( J+1, J ), 1,
$ SUMJ ), TMAX )
END DO
END IF
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm without introducing Infinity
* in the summation
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
END IF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point entry.
* Rely on TRSV to propagate Inf and NaN.
CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the
Expand Down
67 changes: 63 additions & 4 deletions SRC/slatrs.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,8 +264,8 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
* .. External Functions ..
LOGICAL LSAME
INTEGER ISAMAX
REAL SASUM, SDOT, SLAMCH
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH
REAL SASUM, SDOT, SLAMCH, SLANGE
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
* ..
* .. External Subroutines ..
EXTERNAL SAXPY, SSCAL, STRSV, XERBLA
Expand Down Expand Up @@ -343,8 +343,67 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
IF( TMAX.LE.BIGNUM ) THEN
TSCAL = ONE
ELSE
TSCAL = ONE / ( SMLNUM*TMAX )
CALL SSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = ONE / ( SMLNUM*TMAX )
CALL SSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be represented
* as floating-point number. Find the offdiagonal entry A( I, J )
* with the largest absolute value. If this entry is not +/- Infinity,
* use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
TMAX = MAX( SLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
$ TMAX )
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
TMAX = MAX( SLANGE( 'M', N-J, 1, A( J+1, J ), 1,
$ SUMJ ), TMAX )
END DO
END IF
*
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm without introducing Infinity
* in the summation
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * ABS( A( I, J ) )
END DO
END IF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point entry.
* Rely on TRSV to propagate Inf and NaN.
CALL STRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the
Expand Down
70 changes: 68 additions & 2 deletions SRC/zlatrs.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -357,8 +357,74 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
IF( TMAX.LE.BIGNUM*HALF ) THEN
TSCAL = ONE
ELSE
TSCAL = HALF / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
*
* Avoid NaN generation if entries in CNORM exceed the
* overflow threshold
*
IF ( TMAX.LE.DLAMCH('Overflow') ) THEN
* Case 1: All entries in CNORM are valid floating-point numbers
TSCAL = HALF / ( SMLNUM*TMAX )
CALL DSCAL( N, TSCAL, CNORM, 1 )
ELSE
* Case 2: At least one column norm of A cannot be
* represented as a floating-point number. Find the
* maximum offdiagonal absolute value
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
* not +/- Infinity, use this value as TSCAL.
TMAX = ZERO
IF( UPPER ) THEN
*
* A is upper triangular.
*
DO J = 2, N
DO I = 1, J - 1
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
$ ABS( DIMAG(A ( I, J ) ) ) )
END DO
END DO
ELSE
*
* A is lower triangular.
*
DO J = 1, N - 1
DO I = J + 1, N
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
$ ABS( DIMAG(A ( I, J ) ) ) )
END DO
END DO
END IF
*
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
TSCAL = ONE / ( SMLNUM*TMAX )
DO J = 1, N
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
CNORM( J ) = CNORM( J )*TSCAL
ELSE
* Recompute the 1-norm of each column without
* introducing Infinity in the summation.
TSCAL = TWO * TSCAL
CNORM( J ) = ZERO
IF( UPPER ) THEN
DO I = 1, J - 1
CNORM( J ) = CNORM( J ) +
$ TSCAL * CABS2( A( I, J ) )
END DO
ELSE
DO I = J + 1, N
CNORM( J ) = CNORM( J ) +
$ TSCAL * CABS2( A( I, J ) )
END DO
END IF
TSCAL = TSCAL * HALF
END IF
END DO
ELSE
* At least one entry of A is not a valid floating-point
* entry. Rely on TRSV to propagate Inf and NaN.
CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
RETURN
END IF
END IF
END IF
*
* Compute a bound on the computed solution vector to see if the
Expand Down

0 comments on commit 0a7dec1

Please sign in to comment.