From 0a7dec1ff1a4113c7c4195dba3f58d525ed573d0 Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sun, 25 Sep 2022 16:33:46 +0200 Subject: [PATCH] Avoid NaN generation in LATRS Fixes #714 --- SRC/clatrs.f | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++-- SRC/dlatrs.f | 67 ++++++++++++++++++++++++++++++++++++++++++++++--- SRC/slatrs.f | 67 ++++++++++++++++++++++++++++++++++++++++++++++--- SRC/zlatrs.f | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++-- 4 files changed, 262 insertions(+), 12 deletions(-) diff --git a/SRC/clatrs.f b/SRC/clatrs.f index b27be6d701..91334b7066 100644 --- a/SRC/clatrs.f +++ b/SRC/clatrs.f @@ -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 diff --git a/SRC/dlatrs.f b/SRC/dlatrs.f index f5035e7ce1..be156bee20 100644 --- a/SRC/dlatrs.f +++ b/SRC/dlatrs.f @@ -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 @@ -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 diff --git a/SRC/slatrs.f b/SRC/slatrs.f index 994374d38e..0761d656f2 100644 --- a/SRC/slatrs.f +++ b/SRC/slatrs.f @@ -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 @@ -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 diff --git a/SRC/zlatrs.f b/SRC/zlatrs.f index cc977bd7da..2276ace875 100644 --- a/SRC/zlatrs.f +++ b/SRC/zlatrs.f @@ -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