From 7ce9bf5fca3745c91af85ca58daba07c5afa7546 Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sun, 21 Aug 2022 17:08:34 +0200 Subject: [PATCH 1/2] Set SCALE early for robust triangular solvers 1) The docs define SCALE as an output argument. Set SCALE *before* the quick return case to have SCALE defined in all cases. This is how the similar routine TRSYL handles the case already. 2) Remove invocations of LABAD in complex routines and be consistent with their real counterparts, which do not call LABAD. --- SRC/clatbs.f | 9 +++------ SRC/clatrs.f | 9 +++------ SRC/dlatbs.f | 2 +- SRC/dlatrs.f | 2 +- SRC/slatbs.f | 2 +- SRC/slatrs.f | 2 +- SRC/zlatbs.f | 9 +++------ SRC/zlatrs.f | 9 +++------ 8 files changed, 16 insertions(+), 28 deletions(-) diff --git a/SRC/clatbs.f b/SRC/clatbs.f index 606f963d38..97abcadce1 100644 --- a/SRC/clatbs.f +++ b/SRC/clatbs.f @@ -278,7 +278,7 @@ SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, $ CDOTU, CLADIV * .. * .. External Subroutines .. - EXTERNAL CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA + EXTERNAL CAXPY, CSSCAL, CTBSV, SSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL @@ -324,17 +324,14 @@ SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * * Determine machine dependent parameters to control overflow. * - SMLNUM = SLAMCH( 'Safe minimum' ) - BIGNUM = ONE / SMLNUM - CALL SLABAD( SMLNUM, BIGNUM ) - SMLNUM = SMLNUM / SLAMCH( 'Precision' ) + SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/clatrs.f b/SRC/clatrs.f index 946ab80689..b27be6d701 100644 --- a/SRC/clatrs.f +++ b/SRC/clatrs.f @@ -274,7 +274,7 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, $ CDOTU, CLADIV * .. * .. External Subroutines .. - EXTERNAL CAXPY, CSSCAL, CTRSV, SLABAD, SSCAL, XERBLA + EXTERNAL CAXPY, CSSCAL, CTRSV, SSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL @@ -318,17 +318,14 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * * Determine machine dependent parameters to control overflow. * - SMLNUM = SLAMCH( 'Safe minimum' ) - BIGNUM = ONE / SMLNUM - CALL SLABAD( SMLNUM, BIGNUM ) - SMLNUM = SMLNUM / SLAMCH( 'Precision' ) + SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/dlatbs.f b/SRC/dlatbs.f index 4b71d53994..6a812743b0 100644 --- a/SRC/dlatbs.f +++ b/SRC/dlatbs.f @@ -310,6 +310,7 @@ SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * @@ -317,7 +318,6 @@ SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/dlatrs.f b/SRC/dlatrs.f index 43f92911d7..f5035e7ce1 100644 --- a/SRC/dlatrs.f +++ b/SRC/dlatrs.f @@ -304,6 +304,7 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * @@ -311,7 +312,6 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/slatbs.f b/SRC/slatbs.f index 617d0b2f50..77940f8cd9 100644 --- a/SRC/slatbs.f +++ b/SRC/slatbs.f @@ -310,6 +310,7 @@ SUBROUTINE SLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * @@ -317,7 +318,6 @@ SUBROUTINE SLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/slatrs.f b/SRC/slatrs.f index 94e0e88bc6..994374d38e 100644 --- a/SRC/slatrs.f +++ b/SRC/slatrs.f @@ -304,6 +304,7 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * @@ -311,7 +312,6 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/zlatbs.f b/SRC/zlatbs.f index b7b2cb8aec..bdffa1ea98 100644 --- a/SRC/zlatbs.f +++ b/SRC/zlatbs.f @@ -278,7 +278,7 @@ SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, $ ZDOTU, ZLADIV * .. * .. External Subroutines .. - EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV, DLABAD + EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN @@ -324,17 +324,14 @@ SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * * Determine machine dependent parameters to control overflow. * - SMLNUM = DLAMCH( 'Safe minimum' ) - BIGNUM = ONE / SMLNUM - CALL DLABAD( SMLNUM, BIGNUM ) - SMLNUM = SMLNUM / DLAMCH( 'Precision' ) + SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * diff --git a/SRC/zlatrs.f b/SRC/zlatrs.f index 91bb688ece..cc977bd7da 100644 --- a/SRC/zlatrs.f +++ b/SRC/zlatrs.f @@ -274,7 +274,7 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, $ ZDOTU, ZLADIV * .. * .. External Subroutines .. - EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV, DLABAD + EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN @@ -318,17 +318,14 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, * * Quick return if possible * + SCALE = ONE IF( N.EQ.0 ) $ RETURN * * Determine machine dependent parameters to control overflow. * - SMLNUM = DLAMCH( 'Safe minimum' ) - BIGNUM = ONE / SMLNUM - CALL DLABAD( SMLNUM, BIGNUM ) - SMLNUM = SMLNUM / DLAMCH( 'Precision' ) + SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) BIGNUM = ONE / SMLNUM - SCALE = ONE * IF( LSAME( NORMIN, 'N' ) ) THEN * From 0a7dec1ff1a4113c7c4195dba3f58d525ed573d0 Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sun, 25 Sep 2022 16:33:46 +0200 Subject: [PATCH 2/2] 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