Set SCALE early for robust triangular solvers

angsch · angsch · commit 7ce9bf5fca37 · 2022-09-19T17:20:47.000+02:00
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.
diff --git 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
@@ -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
@@ -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
 *
diff --git a/SRC/dlatrs.f b/SRC/dlatrs.f
@@ -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
 *
diff --git a/SRC/slatbs.f b/SRC/slatbs.f
@@ -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
 *
diff --git a/SRC/slatrs.f b/SRC/slatrs.f
@@ -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
 *
diff --git 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
@@ -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
 *

Original file line number	Diff line number	Diff line change
`@@ -278,7 +278,7 @@ SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`278`	`278`	`$ CDOTU, CLADIV`
`279`	`279`	`* ..`
`280`	`280`	`* .. External Subroutines ..`
`281`		`- EXTERNAL CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA`
	`281`	`+ EXTERNAL CAXPY, CSSCAL, CTBSV, SSCAL, XERBLA`
`282`	`282`	`* ..`
`283`	`283`	`* .. Intrinsic Functions ..`
`284`	`284`	`INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL`
`@@ -324,17 +324,14 @@ SUBROUTINE CLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`324`	`324`	`*`
`325`	`325`	`* Quick return if possible`
`326`	`326`	`*`
	`327`	`+ SCALE = ONE`
`327`	`328`	`IF( N.EQ.0 )`
`328`	`329`	`$ RETURN`
`329`	`330`	`*`
`330`	`331`	`* Determine machine dependent parameters to control overflow.`
`331`	`332`	`*`
`332`		`- SMLNUM = SLAMCH( 'Safe minimum' )`
`333`		`- BIGNUM = ONE / SMLNUM`
`334`		`- CALL SLABAD( SMLNUM, BIGNUM )`
`335`		`- SMLNUM = SMLNUM / SLAMCH( 'Precision' )`
	`333`	`+ SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )`
`336`	`334`	`BIGNUM = ONE / SMLNUM`
`337`		`- SCALE = ONE`
`338`	`335`	`*`
`339`	`336`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`340`	`337`	`*`
Original file line number	Diff line number	Diff line change
`@@ -274,7 +274,7 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`274`	`274`	`$ CDOTU, CLADIV`
`275`	`275`	`* ..`
`276`	`276`	`* .. External Subroutines ..`
`277`		`- EXTERNAL CAXPY, CSSCAL, CTRSV, SLABAD, SSCAL, XERBLA`
	`277`	`+ EXTERNAL CAXPY, CSSCAL, CTRSV, SSCAL, XERBLA`
`278`	`278`	`* ..`
`279`	`279`	`* .. Intrinsic Functions ..`
`280`	`280`	`INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL`
`@@ -318,17 +318,14 @@ SUBROUTINE CLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`318`	`318`	`*`
`319`	`319`	`* Quick return if possible`
`320`	`320`	`*`
	`321`	`+ SCALE = ONE`
`321`	`322`	`IF( N.EQ.0 )`
`322`	`323`	`$ RETURN`
`323`	`324`	`*`
`324`	`325`	`* Determine machine dependent parameters to control overflow.`
`325`	`326`	`*`
`326`		`- SMLNUM = SLAMCH( 'Safe minimum' )`
`327`		`- BIGNUM = ONE / SMLNUM`
`328`		`- CALL SLABAD( SMLNUM, BIGNUM )`
`329`		`- SMLNUM = SMLNUM / SLAMCH( 'Precision' )`
	`327`	`+ SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )`
`330`	`328`	`BIGNUM = ONE / SMLNUM`
`331`		`- SCALE = ONE`
`332`	`329`	`*`
`333`	`330`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`334`	`331`	`*`
Original file line number	Diff line number	Diff line change
`@@ -310,14 +310,14 @@ SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`310`	`310`	`*`
`311`	`311`	`* Quick return if possible`
`312`	`312`	`*`
	`313`	`+ SCALE = ONE`
`313`	`314`	`IF( N.EQ.0 )`
`314`	`315`	`$ RETURN`
`315`	`316`	`*`
`316`	`317`	`* Determine machine dependent parameters to control overflow.`
`317`	`318`	`*`
`318`	`319`	`SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )`
`319`	`320`	`BIGNUM = ONE / SMLNUM`
`320`		`- SCALE = ONE`
`321`	`321`	`*`
`322`	`322`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`323`	`323`	`*`
Original file line number	Diff line number	Diff line change
`@@ -304,14 +304,14 @@ SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`304`	`304`	`*`
`305`	`305`	`* Quick return if possible`
`306`	`306`	`*`
	`307`	`+ SCALE = ONE`
`307`	`308`	`IF( N.EQ.0 )`
`308`	`309`	`$ RETURN`
`309`	`310`	`*`
`310`	`311`	`* Determine machine dependent parameters to control overflow.`
`311`	`312`	`*`
`312`	`313`	`SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )`
`313`	`314`	`BIGNUM = ONE / SMLNUM`
`314`		`- SCALE = ONE`
`315`	`315`	`*`
`316`	`316`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`317`	`317`	`*`
Original file line number	Diff line number	Diff line change
`@@ -310,14 +310,14 @@ SUBROUTINE SLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`310`	`310`	`*`
`311`	`311`	`* Quick return if possible`
`312`	`312`	`*`
	`313`	`+ SCALE = ONE`
`313`	`314`	`IF( N.EQ.0 )`
`314`	`315`	`$ RETURN`
`315`	`316`	`*`
`316`	`317`	`* Determine machine dependent parameters to control overflow.`
`317`	`318`	`*`
`318`	`319`	`SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )`
`319`	`320`	`BIGNUM = ONE / SMLNUM`
`320`		`- SCALE = ONE`
`321`	`321`	`*`
`322`	`322`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`323`	`323`	`*`
Original file line number	Diff line number	Diff line change
`@@ -304,14 +304,14 @@ SUBROUTINE SLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`304`	`304`	`*`
`305`	`305`	`* Quick return if possible`
`306`	`306`	`*`
	`307`	`+ SCALE = ONE`
`307`	`308`	`IF( N.EQ.0 )`
`308`	`309`	`$ RETURN`
`309`	`310`	`*`
`310`	`311`	`* Determine machine dependent parameters to control overflow.`
`311`	`312`	`*`
`312`	`313`	`SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )`
`313`	`314`	`BIGNUM = ONE / SMLNUM`
`314`		`- SCALE = ONE`
`315`	`315`	`*`
`316`	`316`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`317`	`317`	`*`
Original file line number	Diff line number	Diff line change
`@@ -278,7 +278,7 @@ SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`278`	`278`	`$ ZDOTU, ZLADIV`
`279`	`279`	`* ..`
`280`	`280`	`* .. External Subroutines ..`
`281`		`- EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV, DLABAD`
	`281`	`+ EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV`
`282`	`282`	`* ..`
`283`	`283`	`* .. Intrinsic Functions ..`
`284`	`284`	`INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN`
`@@ -324,17 +324,14 @@ SUBROUTINE ZLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X,`
`324`	`324`	`*`
`325`	`325`	`* Quick return if possible`
`326`	`326`	`*`
	`327`	`+ SCALE = ONE`
`327`	`328`	`IF( N.EQ.0 )`
`328`	`329`	`$ RETURN`
`329`	`330`	`*`
`330`	`331`	`* Determine machine dependent parameters to control overflow.`
`331`	`332`	`*`
`332`		`- SMLNUM = DLAMCH( 'Safe minimum' )`
`333`		`- BIGNUM = ONE / SMLNUM`
`334`		`- CALL DLABAD( SMLNUM, BIGNUM )`
`335`		`- SMLNUM = SMLNUM / DLAMCH( 'Precision' )`
	`333`	`+ SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )`
`336`	`334`	`BIGNUM = ONE / SMLNUM`
`337`		`- SCALE = ONE`
`338`	`335`	`*`
`339`	`336`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`340`	`337`	`*`
Original file line number	Diff line number	Diff line change
`@@ -274,7 +274,7 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`274`	`274`	`$ ZDOTU, ZLADIV`
`275`	`275`	`* ..`
`276`	`276`	`* .. External Subroutines ..`
`277`		`- EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV, DLABAD`
	`277`	`+ EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV`
`278`	`278`	`* ..`
`279`	`279`	`* .. Intrinsic Functions ..`
`280`	`280`	`INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN`
`@@ -318,17 +318,14 @@ SUBROUTINE ZLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,`
`318`	`318`	`*`
`319`	`319`	`* Quick return if possible`
`320`	`320`	`*`
	`321`	`+ SCALE = ONE`
`321`	`322`	`IF( N.EQ.0 )`
`322`	`323`	`$ RETURN`
`323`	`324`	`*`
`324`	`325`	`* Determine machine dependent parameters to control overflow.`
`325`	`326`	`*`
`326`		`- SMLNUM = DLAMCH( 'Safe minimum' )`
`327`		`- BIGNUM = ONE / SMLNUM`
`328`		`- CALL DLABAD( SMLNUM, BIGNUM )`
`329`		`- SMLNUM = SMLNUM / DLAMCH( 'Precision' )`
	`327`	`+ SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )`
`330`	`328`	`BIGNUM = ONE / SMLNUM`
`331`		`- SCALE = ONE`
`332`	`329`	`*`
`333`	`330`	`IF( LSAME( NORMIN, 'N' ) ) THEN`
`334`	`331`	`*`