From 7ce9bf5fca3745c91af85ca58daba07c5afa7546 Mon Sep 17 00:00:00 2001 From: Angelika Schwarz Date: Sun, 21 Aug 2022 17:08:34 +0200 Subject: [PATCH] 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 *