Skip to content

Set SCALE early for robust triangular solvers #712

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Sep 29, 2022
Merged

Set SCALE early for robust triangular solvers #712

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7ce9bf5
Set SCALE early for robust triangular solvers
angsch Aug 21, 2022
0a7dec1
Avoid NaN generation in LATRS
angsch Sep 25, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 3 additions & 6 deletions SRC/clatbs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
*
Expand Down
79 changes: 71 additions & 8 deletions SRC/clatrs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
*
Expand Down Expand Up @@ -360,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
2 changes: 1 addition & 1 deletion SRC/dlatbs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,14 +310,14 @@ SUBROUTINE DLATBS( 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' ) / DLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
Expand Down
69 changes: 64 additions & 5 deletions SRC/dlatrs.f
View file
Open in desktop
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 @@ -304,14 +304,14 @@ SUBROUTINE DLATRS( 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' ) / DLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
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
2 changes: 1 addition & 1 deletion SRC/slatbs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,14 +310,14 @@ SUBROUTINE SLATBS( 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' ) / SLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
Expand Down
69 changes: 64 additions & 5 deletions SRC/slatrs.f
View file
Open in desktop
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 @@ -304,14 +304,14 @@ SUBROUTINE SLATRS( 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' ) / SLAMCH( 'Precision' )
BIGNUM = ONE / SMLNUM
SCALE = ONE
*
IF( LSAME( NORMIN, 'N' ) ) THEN
*
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
9 changes: 3 additions & 6 deletions SRC/zlatbs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
*
Expand Down
Loading