From 97cfda26ea07b83a106077866b2f170e22bcc1f6 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 20:55:21 -0700 Subject: [PATCH 1/3] Scale down in xLARFG and xLARFGP when |BETA| exceeds half the overflow threshold The Householder generators rescale their input only at the small end, when |BETA| < SAFMIN. At the large end they form ALPHA-BETA (xLARFG) and ALPHA+BETA (xLARFGP), sums of two like-signed terms each bounded by |BETA|, which overflow whenever |BETA| > OVFL/2 although BETA itself and the reflector are representable. xLARFG then returns TAU = Inf and v = 0; xLARFGP gets ALPHA = Inf, TAU = 0, takes its flush branch and returns H = I with the tail of the column not annihilated. Every QR-type factorization built on them (xGEQRF, xGEQRFP, xGELQF, xGEQLF, xGERQF, xGEQRT, xGEQR, xGEQP3, xGEQP3RK, xTZRZF) returns a wrong R or Inf/NaN with INFO = 0 for a column whose leading entry exceeds about 9e307 (1.7e38 in single precision). The least-squares drivers pre-scale A and are not affected. When |BETA| > OVFL/2, scale X and ALPHA by SAFMIN, recompute XNORM and BETA, and multiply BETA back by 1/SAFMIN on exit, mirroring the existing small-end branch. The factors are powers of two, so TAU and v are the same reflector. Since |ALPHA| <= |BETA|, the sums cannot overflow for |BETA| <= OVFL/2, and the test is placed there so that no input the old code handled takes the new branch. The overflow threshold is fetched from xLAMCH( 'O' ) only when |BETA|*SAFMIN > 1, which keeps the extra cost on the common path to one multiply and one compare. The QR, LQ, QL and RQ test paths get a matrix type for it: type 9 is the random matrix of type 4 with the entry the first reflector works on raised to three quarters of the overflow threshold, A( 1, 1 ) for QR and LQ and A( M, N ) for QL and RQ. The generator cannot produce such a matrix, because it scales what it generates by the requested norm. The existing test ratios detect the failure, but only once the comparison with the threshold is guarded with xISNAN: a NaN is not .GE. THRESH, so without the guard every ratio of the new type passes on the parent commit. With it, the parent fails 5670 ratios per precision and this branch none. Over a sweep of 856704 (precision, generator, n, INCX, exponent of ALPHA, exponent of X, sign, pattern) cases the outputs are bit-identical to the parent commit whenever the true |BETA| <= OVFL/2; the parent returns Inf, NaN or a reflector with residual above tolerance in 12622 cases with a representable BETA, this branch in none, and each reflector was verified against its defining relation in quadruple precision. The full LAPACK test suite passes: 5506497 LAPACK and 315872 BLAS tests, 0 numerical errors, 0 other errors; the 64596 tests above the parent are the four new matrix types. Co-Authored-By: Claude Fable 5.1 --- SRC/clarfg.f | 24 ++++++++++++++++++++---- SRC/clarfgp.f | 24 ++++++++++++++++++++---- SRC/dlarfg.f | 23 +++++++++++++++++++---- SRC/dlarfgp.f | 23 +++++++++++++++++++---- SRC/slarfg.f | 23 +++++++++++++++++++---- SRC/slarfgp.f | 23 +++++++++++++++++++---- SRC/zlarfg.f | 24 ++++++++++++++++++++---- SRC/zlarfgp.f | 24 ++++++++++++++++++++---- TESTING/LIN/alahd.f | 3 ++- TESTING/LIN/cchkaa.F | 8 ++++---- TESTING/LIN/cchklq.f | 23 +++++++++++++++++++---- TESTING/LIN/cchkql.f | 24 ++++++++++++++++++++---- TESTING/LIN/cchkqr.f | 21 +++++++++++++++++---- TESTING/LIN/cchkrq.f | 24 ++++++++++++++++++++---- TESTING/LIN/dchkaa.F | 8 ++++---- TESTING/LIN/dchklq.f | 23 +++++++++++++++++++---- TESTING/LIN/dchkql.f | 24 ++++++++++++++++++++---- TESTING/LIN/dchkqr.f | 21 +++++++++++++++++---- TESTING/LIN/dchkrq.f | 24 ++++++++++++++++++++---- TESTING/LIN/schkaa.F | 8 ++++---- TESTING/LIN/schklq.f | 23 +++++++++++++++++++---- TESTING/LIN/schkql.f | 24 ++++++++++++++++++++---- TESTING/LIN/schkqr.f | 21 +++++++++++++++++---- TESTING/LIN/schkrq.f | 24 ++++++++++++++++++++---- TESTING/LIN/zchkaa.F | 8 ++++---- TESTING/LIN/zchklq.f | 23 +++++++++++++++++++---- TESTING/LIN/zchkql.f | 24 ++++++++++++++++++++---- TESTING/LIN/zchkqr.f | 21 +++++++++++++++++---- TESTING/LIN/zchkrq.f | 24 ++++++++++++++++++++---- TESTING/ctest.in | 8 ++++---- TESTING/dtest.in | 8 ++++---- TESTING/stest.in | 8 ++++---- TESTING/ztest.in | 8 ++++---- 33 files changed, 494 insertions(+), 129 deletions(-) diff --git a/SRC/clarfg.f b/SRC/clarfg.f index d9e6cf4bc..b77be5446 100644 --- a/SRC/clarfg.f +++ b/SRC/clarfg.f @@ -118,8 +118,9 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL ONE, ZERO - PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL ONE, ZERO, HALF + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN + INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN * .. * .. External Subroutines .. EXTERNAL CSCAL, CSSCAL @@ -179,16 +180,31 @@ SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) XNORM = SCNRM2( N-1, X, INCX ) ALPHA = CMPLX( ALPHR, ALPHI ) BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL CSSCAL( N-1, SAFMIN, X, INCX ) + ALPHI = ALPHI*SAFMIN + ALPHR = ALPHR*SAFMIN + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA ) CALL CSCAL( N-1, ALPHA, X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/clarfgp.f b/SRC/clarfgp.f index ca24dc8ec..e0b4d6d68 100644 --- a/SRC/clarfgp.f +++ b/SRC/clarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL TWO, ONE, ZERO - PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) $ CLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN + INTRINSIC ABS, AIMAG, CMPLX, HUGE, REAL, SIGN * .. * .. External Subroutines .. EXTERNAL CSCAL, CSSCAL @@ -193,6 +194,18 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) XNORM = SCNRM2( N-1, X, INCX ) ALPHA = CMPLX( ALPHR, ALPHI ) BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL CSSCAL( N-1, SMLNUM, X, INCX ) + ALPHI = ALPHI*SMLNUM + ALPHR = ALPHR*SMLNUM + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -245,11 +258,14 @@ SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/dlarfg.f b/SRC/dlarfg.f index 38f65af0f..03c489b89 100644 --- a/SRC/dlarfg.f +++ b/SRC/dlarfg.f @@ -118,8 +118,9 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE, ZERO, HALF + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -130,7 +131,7 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL @@ -173,15 +174,29 @@ SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) * XNORM = DNRM2( N-1, X, INCX ) BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + RSAFMN = ONE / SAFMIN + KNT = -1 + CALL DSCAL( N-1, SAFMIN, X, INCX ) + ALPHA = ALPHA*SAFMIN + XNORM = DNRM2( N-1, X, INCX ) + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) END IF TAU = ( BETA-ALPHA ) / BETA CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/dlarfgp.f b/SRC/dlarfgp.f index fe6613997..c74daa312 100644 --- a/SRC/dlarfgp.f +++ b/SRC/dlarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION TWO, ONE, ZERO - PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -128,7 +129,7 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL @@ -185,6 +186,17 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * XNORM = DNRM2( N-1, X, INCX ) BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + BIGNUM = ONE / SMLNUM + KNT = -1 + CALL DSCAL( N-1, SMLNUM, X, INCX ) + ALPHA = ALPHA*SMLNUM + XNORM = DNRM2( N-1, X, INCX ) + BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -229,11 +241,14 @@ SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/slarfg.f b/SRC/slarfg.f index a0893439d..2e7163c81 100644 --- a/SRC/slarfg.f +++ b/SRC/slarfg.f @@ -118,8 +118,9 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL ONE, ZERO - PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL ONE, ZERO, HALF + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -130,7 +131,7 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL SLAMCH, SLAPY2, SNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL SSCAL @@ -173,15 +174,29 @@ SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) * XNORM = SNRM2( N-1, X, INCX ) BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + RSAFMN = ONE / SAFMIN + KNT = -1 + CALL SSCAL( N-1, SAFMIN, X, INCX ) + ALPHA = ALPHA*SAFMIN + XNORM = SNRM2( N-1, X, INCX ) + BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) END IF TAU = ( BETA-ALPHA ) / BETA CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/slarfgp.f b/SRC/slarfgp.f index 0a3d72313..d23a74b95 100644 --- a/SRC/slarfgp.f +++ b/SRC/slarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - REAL TWO, ONE, ZERO - PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) + REAL TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0, + $ HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -128,7 +129,7 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) EXTERNAL SLAMCH, SLAPY2, SNRM2 * .. * .. Intrinsic Functions .. - INTRINSIC ABS, SIGN + INTRINSIC ABS, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL SSCAL @@ -185,6 +186,17 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * XNORM = SNRM2( N-1, X, INCX ) BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + BIGNUM = ONE / SMLNUM + KNT = -1 + CALL SSCAL( N-1, SMLNUM, X, INCX ) + ALPHA = ALPHA*SMLNUM + XNORM = SNRM2( N-1, X, INCX ) + BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -229,11 +241,14 @@ SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/SRC/zlarfg.f b/SRC/zlarfg.f index 576d9c64d..ece8ed5c9 100644 --- a/SRC/zlarfg.f +++ b/SRC/zlarfg.f @@ -118,8 +118,9 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE, ZERO, HALF + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN + INTRINSIC ABS, DBLE, DCMPLX, DIMAG, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZSCAL @@ -179,16 +180,31 @@ SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) XNORM = DZNRM2( N-1, X, INCX ) ALPHA = DCMPLX( ALPHR, ALPHI ) BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA-BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL ZDSCAL( N-1, SAFMIN, X, INCX ) + ALPHI = ALPHI*SAFMIN + ALPHR = ALPHR*SAFMIN + XNORM = DZNRM2( N-1, X, INCX ) + ALPHA = DCMPLX( ALPHR, ALPHI ) + BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) CALL ZSCAL( N-1, ALPHA, X, INCX ) * -* If ALPHA is subnormal, it may lose relative accuracy +* Undo the scaling. If ALPHA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*RSAFMN ALPHA = BETA END IF * diff --git a/SRC/zlarfgp.f b/SRC/zlarfgp.f index 29a6b824f..87d49259e 100644 --- a/SRC/zlarfgp.f +++ b/SRC/zlarfgp.f @@ -116,8 +116,9 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION TWO, ONE, ZERO - PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) + DOUBLE PRECISION TWO, ONE, ZERO, HALF + PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0, + $ HALF = 0.5D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT @@ -131,7 +132,7 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) $ ZLADIV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN + INTRINSIC ABS, DBLE, DCMPLX, DIMAG, HUGE, SIGN * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZSCAL @@ -193,6 +194,18 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) XNORM = DZNRM2( N-1, X, INCX ) ALPHA = DCMPLX( ALPHR, ALPHI ) BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + ELSE IF( ABS( BETA ).GT.HALF*HUGE( ZERO ) ) THEN +* +* |ALPHA| <= |BETA|, so ALPHA+BETA can overflow only when +* |BETA| > HUGE/2; scale X down and recompute them. +* + KNT = -1 + CALL ZDSCAL( N-1, SMLNUM, X, INCX ) + ALPHI = ALPHI*SMLNUM + ALPHR = ALPHR*SMLNUM + XNORM = DZNRM2( N-1, X, INCX ) + ALPHA = DCMPLX( ALPHR, ALPHI ) + BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA @@ -245,11 +258,14 @@ SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) * END IF * -* If BETA is subnormal, it may lose relative accuracy +* Undo the scaling. If BETA is subnormal, it may lose relative +* accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE + IF( KNT.LT.0 ) + $ BETA = BETA*BIGNUM ALPHA = BETA END IF * diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index b04a3f796..85d216b3a 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -945,7 +945,8 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) $ '2. Upper triangular', 16X, '6. Random, CNDNUM = 0.1/EPS', $ / 4X, '3. Lower triangular', 16X, $ '7. Scaled near underflow', / 4X, '4. Random, CNDNUM = 2', - $ 14X, '8. Scaled near overflow' ) + $ 14X, '8. Scaled near overflow', / 39X, + $ '9. Leading entry near overflow' ) * * QP matrix types * diff --git a/TESTING/LIN/cchkaa.F b/TESTING/LIN/cchkaa.F index ac181df8c..a46fe949c 100644 --- a/TESTING/LIN/cchkaa.F +++ b/TESTING/LIN/cchkaa.F @@ -1038,7 +1038,7 @@ PROGRAM CCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1055,7 +1055,7 @@ PROGRAM CCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1072,7 +1072,7 @@ PROGRAM CCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1089,7 +1089,7 @@ PROGRAM CCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/cchklq.f b/TESTING/LIN/cchklq.f index 175fd1a84..e9648ff25 100644 --- a/TESTING/LIN/cchklq.f +++ b/TESTING/LIN/cchklq.f @@ -219,9 +219,10 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRLQ, CGELQF, CGELS, $ CGET02, CLACPY, CLARHS, CLATB4, CLATMS, CLQT01, @@ -314,6 +320,14 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CLQT01; other values are * used in the calls of CLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE CCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkql.f b/TESTING/LIN/cchkql.f index 8b029263b..db3c2bf94 100644 --- a/TESTING/LIN/cchkql.f +++ b/TESTING/LIN/cchkql.f @@ -219,9 +219,10 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRQL, CGEQLS, CGET02, $ CLACPY, CLARHS, CLATB4, CLATMS, CQLT01, CQLT02, @@ -314,6 +320,15 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CQLT01; other values are * used in the calls of CQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE CCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkqr.f b/TESTING/LIN/cchkqr.f index 2f3e64000..556ade52d 100644 --- a/TESTING/LIN/cchkqr.f +++ b/TESTING/LIN/cchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH LOGICAL CGENND EXTERNAL CGENND + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRQR, CGELS, CGET02, @@ -323,6 +327,14 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CQRT01; other values are * used in the calls of CQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE CCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/cchkrq.f b/TESTING/LIN/cchkrq.f index 4d77d192f..4fb7f6c5a 100644 --- a/TESTING/LIN/cchkrq.f +++ b/TESTING/LIN/cchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CERRRQ, CGERQS, CGET02, $ CLACPY, CLARHS, CLATB4, CLATMS, CRQT01, CRQT02, @@ -319,6 +325,15 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of CRQT01; other values are * used in the calls of CRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE CCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkaa.F b/TESTING/LIN/dchkaa.F index d0f7dbcb4..1411d5633 100644 --- a/TESTING/LIN/dchkaa.F +++ b/TESTING/LIN/dchkaa.F @@ -880,7 +880,7 @@ PROGRAM DCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -897,7 +897,7 @@ PROGRAM DCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -914,7 +914,7 @@ PROGRAM DCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -931,7 +931,7 @@ PROGRAM DCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/dchklq.f b/TESTING/LIN/dchklq.f index 8fa7c8d6f..d297af90d 100644 --- a/TESTING/LIN/dchklq.f +++ b/TESTING/LIN/dchklq.f @@ -219,9 +219,10 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRLQ, DGELQF, DGELS, $ DGET02, DLACPY, DLARHS, DLATB4, DLATMS, DLQT01, @@ -314,6 +320,14 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DLQT01; other values are * used in the calls of DLQT02, and must not exceed MINMN. @@ -433,7 +447,8 @@ SUBROUTINE DCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkql.f b/TESTING/LIN/dchkql.f index 24837f1c4..0f907b827 100644 --- a/TESTING/LIN/dchkql.f +++ b/TESTING/LIN/dchkql.f @@ -219,9 +219,10 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRQL, DGEQLS, DGET02, $ DLACPY, DLARHS, DLATB4, DLATMS, DQLT01, DQLT02, @@ -314,6 +320,15 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DQLT01; other values are * used in the calls of DQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE DCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkqr.f b/TESTING/LIN/dchkqr.f index 013f7676a..157e9b5cb 100644 --- a/TESTING/LIN/dchkqr.f +++ b/TESTING/LIN/dchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH LOGICAL DGENND EXTERNAL DGENND + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRQR, DGELS, DGET02, @@ -323,6 +327,14 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DQRT01; other values are * used in the calls of DQRT02, and must not exceed MINMN. @@ -435,7 +447,8 @@ SUBROUTINE DCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/dchkrq.f b/TESTING/LIN/dchkrq.f index e85ac113b..3d2d663d5 100644 --- a/TESTING/LIN/dchkrq.f +++ b/TESTING/LIN/dchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DERRRQ, DGERQS, DGET02, $ DLACPY, DLARHS, DLATB4, DLATMS, DRQT01, DRQT02, @@ -319,6 +325,15 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of DRQT01; other values are * used in the calls of DRQT02, and must not exceed MINMN. @@ -416,7 +431,8 @@ SUBROUTINE DCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkaa.F b/TESTING/LIN/schkaa.F index dd3f8de4b..d7a0390f5 100644 --- a/TESTING/LIN/schkaa.F +++ b/TESTING/LIN/schkaa.F @@ -874,7 +874,7 @@ PROGRAM SCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -891,7 +891,7 @@ PROGRAM SCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -908,7 +908,7 @@ PROGRAM SCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -925,7 +925,7 @@ PROGRAM SCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/schklq.f b/TESTING/LIN/schklq.f index 5dce71ecd..ae2a2ca23 100644 --- a/TESTING/LIN/schklq.f +++ b/TESTING/LIN/schklq.f @@ -219,9 +219,10 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRLQ, SGELQF, SGELS, $ SGET02, SLACPY, SLARHS, SLATB4, SLATMS, SLQT01, @@ -314,6 +320,14 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SLQT01; other values are * used in the calls of SLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE SCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkql.f b/TESTING/LIN/schkql.f index 20726e188..7ea38fb0b 100644 --- a/TESTING/LIN/schkql.f +++ b/TESTING/LIN/schkql.f @@ -219,9 +219,10 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRQL, SGEQLS, SGET02, $ SLACPY, SLARHS, SLATB4, SLATMS, SQLT01, SQLT02, @@ -314,6 +320,15 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SQLT01; other values are * used in the calls of SQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE SCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkqr.f b/TESTING/LIN/schkqr.f index 49d775c17..a680e7b84 100644 --- a/TESTING/LIN/schkqr.f +++ b/TESTING/LIN/schkqr.f @@ -224,9 +224,10 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH LOGICAL SGENND EXTERNAL SGENND + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRQR, SGELS, SGET02, @@ -323,6 +327,14 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SQRT01; other values are * used in the calls of SQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE SCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/schkrq.f b/TESTING/LIN/schkrq.f index 1d919f19c..5bd14fdc8 100644 --- a/TESTING/LIN/schkrq.f +++ b/TESTING/LIN/schkrq.f @@ -224,9 +224,10 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - REAL ZERO - PARAMETER ( ZERO = 0.0E0 ) + PARAMETER ( NTYPES = 9 ) + REAL ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, + $ QUARTER = 0.25E0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) REAL RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH + EXTERNAL SISNAN, SLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SERRRQ, SGERQS, SGET02, $ SLACPY, SLARHS, SLATB4, SLATMS, SRQT01, SRQT02, @@ -319,6 +325,15 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ SLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of SRQT01; other values are * used in the calls of SRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE SCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ SISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkaa.F b/TESTING/LIN/zchkaa.F index fb1b5cd59..b64e3f31d 100644 --- a/TESTING/LIN/zchkaa.F +++ b/TESTING/LIN/zchkaa.F @@ -1037,7 +1037,7 @@ PROGRAM ZCHKAA * * QR: QR factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1054,7 +1054,7 @@ PROGRAM ZCHKAA * * LQ: LQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1071,7 +1071,7 @@ PROGRAM ZCHKAA * * QL: QL factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN @@ -1088,7 +1088,7 @@ PROGRAM ZCHKAA * * RQ: RQ factorization * - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/zchklq.f b/TESTING/LIN/zchklq.f index d4106fa4c..1305ba8ca 100644 --- a/TESTING/LIN/zchklq.f +++ b/TESTING/LIN/zchklq.f @@ -219,9 +219,10 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, ZERRLQ, ZGELQF, ZGELS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZLQT01, @@ -314,6 +320,14 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZLQT01; other values are * used in the calls of ZLQT02, and must not exceed MINMN. @@ -423,7 +437,8 @@ SUBROUTINE ZCHKLQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkql.f b/TESTING/LIN/zchkql.f index cafb2c623..3f8b005d4 100644 --- a/TESTING/LIN/zchkql.f +++ b/TESTING/LIN/zchkql.f @@ -219,9 +219,10 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -235,6 +236,11 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQL, ZGEQLS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZQLT01, @@ -314,6 +320,15 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZQLT01; other values are * used in the calls of ZQLT02, and must not exceed MINMN. @@ -410,7 +425,8 @@ SUBROUTINE ZCHKQL( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkqr.f b/TESTING/LIN/zchkqr.f index 740abf4df..6cb7d2c39 100644 --- a/TESTING/LIN/zchkqr.f +++ b/TESTING/LIN/zchkqr.f @@ -224,9 +224,10 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 9 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -241,8 +242,11 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH LOGICAL ZGENND EXTERNAL ZGENND + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRQR, ZGELS, @@ -323,6 +327,14 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( 1 ) = ( ONE - QUARTER )*DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZQRT01; other values are * used in the calls of ZQRT02, and must not exceed MINMN. @@ -434,7 +446,8 @@ SUBROUTINE ZCHKQR( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NTESTS - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/LIN/zchkrq.f b/TESTING/LIN/zchkrq.f index 6bcb3d6b9..7e2ccc80a 100644 --- a/TESTING/LIN/zchkrq.f +++ b/TESTING/LIN/zchkrq.f @@ -224,9 +224,10 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER NTESTS PARAMETER ( NTESTS = 7 ) INTEGER NTYPES - PARAMETER ( NTYPES = 8 ) - DOUBLE PRECISION ZERO - PARAMETER ( ZERO = 0.0D0 ) + PARAMETER ( NTYPES = 9 ) + DOUBLE PRECISION ZERO, ONE, QUARTER + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, + $ QUARTER = 0.25D0 ) * .. * .. Local Scalars .. CHARACTER DIST, TYPE @@ -240,6 +241,11 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, INTEGER ISEED( 4 ), ISEEDY( 4 ), KVAL( 4 ) DOUBLE PRECISION RESULT( NTESTS ) * .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH + EXTERNAL DISNAN, DLAMCH +* .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZERRRQ, ZGERQS, $ ZGET02, ZLACPY, ZLARHS, ZLATB4, ZLATMS, ZRQT01, @@ -319,6 +325,15 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, GO TO 50 END IF * +* Type 9: make the entry the first reflector works on +* large enough that its sum with the norm of the column +* overflows. The generator cannot produce such a matrix, +* because it scales the matrix by its norm. +* + IF( IMAT.EQ.9 .AND. MINMN.GT.0 ) + $ A( ( N-1 )*LDA+M ) = ( ONE-QUARTER )* + $ DLAMCH( 'Overflow' ) +* * Set some values for K: the first value must be MINMN, * corresponding to the call of ZRQT01; other values are * used in the calls of ZRQT02, and must not exceed MINMN. @@ -415,7 +430,8 @@ SUBROUTINE ZCHKRQ( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NXVAL, * pass the threshold. * DO 20 I = 1, NT - IF( RESULT( I ).GE.THRESH ) THEN + IF( RESULT( I ).GE.THRESH .OR. + $ DISNAN( RESULT( I ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9999 )M, N, K, NB, NX, diff --git a/TESTING/ctest.in b/TESTING/ctest.in index 4e30224d7..a16b56295 100644 --- a/TESTING/ctest.in +++ b/TESTING/ctest.in @@ -37,10 +37,10 @@ CSP 11 List types on next line if 0 < NTYPES < 11 CTR 18 List types on next line if 0 < NTYPES < 18 CTP 18 List types on next line if 0 < NTYPES < 18 CTB 17 List types on next line if 0 < NTYPES < 17 -CQR 8 List types on next line if 0 < NTYPES < 8 -CRQ 8 List types on next line if 0 < NTYPES < 8 -CLQ 8 List types on next line if 0 < NTYPES < 8 -CQL 8 List types on next line if 0 < NTYPES < 8 +CQR 9 List types on next line if 0 < NTYPES < 9 +CRQ 9 List types on next line if 0 < NTYPES < 9 +CLQ 9 List types on next line if 0 < NTYPES < 9 +CQL 9 List types on next line if 0 < NTYPES < 9 CQP 6 List types on next line if 0 < NTYPES < 6 CQK 19 List types on next line if 0 < NTYPES < 19 CCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/dtest.in b/TESTING/dtest.in index cde62db50..2c5ce7659 100644 --- a/TESTING/dtest.in +++ b/TESTING/dtest.in @@ -31,10 +31,10 @@ DSP 10 List types on next line if 0 < NTYPES < 10 DTR 18 List types on next line if 0 < NTYPES < 18 DTP 18 List types on next line if 0 < NTYPES < 18 DTB 17 List types on next line if 0 < NTYPES < 17 -DQR 8 List types on next line if 0 < NTYPES < 8 -DRQ 8 List types on next line if 0 < NTYPES < 8 -DLQ 8 List types on next line if 0 < NTYPES < 8 -DQL 8 List types on next line if 0 < NTYPES < 8 +DQR 9 List types on next line if 0 < NTYPES < 9 +DRQ 9 List types on next line if 0 < NTYPES < 9 +DLQ 9 List types on next line if 0 < NTYPES < 9 +DQL 9 List types on next line if 0 < NTYPES < 9 DQP 6 List types on next line if 0 < NTYPES < 6 DQK 19 LIst types on next line if 0 < NTYPES < 19 DCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/stest.in b/TESTING/stest.in index abfd639fd..b6a859ef1 100644 --- a/TESTING/stest.in +++ b/TESTING/stest.in @@ -31,10 +31,10 @@ SSP 10 List types on next line if 0 < NTYPES < 10 STR 18 List types on next line if 0 < NTYPES < 18 STP 18 List types on next line if 0 < NTYPES < 18 STB 17 List types on next line if 0 < NTYPES < 17 -SQR 8 List types on next line if 0 < NTYPES < 8 -SRQ 8 List types on next line if 0 < NTYPES < 8 -SLQ 8 List types on next line if 0 < NTYPES < 8 -SQL 8 List types on next line if 0 < NTYPES < 8 +SQR 9 List types on next line if 0 < NTYPES < 9 +SRQ 9 List types on next line if 0 < NTYPES < 9 +SLQ 9 List types on next line if 0 < NTYPES < 9 +SQL 9 List types on next line if 0 < NTYPES < 9 SQP 6 List types on next line if 0 < NTYPES < 6 SQK 19 List types on next line if 0 < NTYPES < 19 SCX 19 LIst types on next line if 0 < NTYPES < 19 diff --git a/TESTING/ztest.in b/TESTING/ztest.in index bf4c9d100..f123e6630 100644 --- a/TESTING/ztest.in +++ b/TESTING/ztest.in @@ -37,10 +37,10 @@ ZSP 11 List types on next line if 0 < NTYPES < 11 ZTR 18 List types on next line if 0 < NTYPES < 18 ZTP 18 List types on next line if 0 < NTYPES < 18 ZTB 17 List types on next line if 0 < NTYPES < 17 -ZQR 8 List types on next line if 0 < NTYPES < 8 -ZRQ 8 List types on next line if 0 < NTYPES < 8 -ZLQ 8 List types on next line if 0 < NTYPES < 8 -ZQL 8 List types on next line if 0 < NTYPES < 8 +ZQR 9 List types on next line if 0 < NTYPES < 9 +ZRQ 9 List types on next line if 0 < NTYPES < 9 +ZLQ 9 List types on next line if 0 < NTYPES < 9 +ZQL 9 List types on next line if 0 < NTYPES < 9 ZQP 6 List types on next line if 0 < NTYPES < 6 ZQK 19 List types on next line if 0 < NTYPES < 19 ZCX 19 List types on next line if 0 < NTYPES < 19 From f9d174cc80335246e198c990a12432b4e9cd5d57 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Sat, 12 Sep 2026 12:20:52 -0700 Subject: [PATCH 2/3] TESTING: Scale large right hand sides in xGEQLS The top-end QL cases can overflow the tau*w update while applying Q to B, even when the reflectors and the final solution are representable. Scale large right hand sides before applying Q and restore their scale after the triangular solve. Use component magnitudes for complex B so the scaling decision itself cannot overflow. This preserves the extreme-value factorization cases and fixes their least-squares residual failures. All four linear test suites pass with Flang 20; GNU 13 passes with both 32-bit and 64-bit integer APIs. --- TESTING/LIN/cgeqls.f | 36 +++++++++++++++++++++++++++++++++--- TESTING/LIN/dgeqls.f | 28 ++++++++++++++++++++++++++-- TESTING/LIN/sgeqls.f | 28 ++++++++++++++++++++++++++-- TESTING/LIN/zgeqls.f | 36 +++++++++++++++++++++++++++++++++--- 4 files changed, 118 insertions(+), 10 deletions(-) diff --git a/TESTING/LIN/cgeqls.f b/TESTING/LIN/cgeqls.f index 232b5dd23..da9e95c58 100644 --- a/TESTING/LIN/cgeqls.f +++ b/TESTING/LIN/cgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by CGEQLF. +*> computed by CGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,11 +140,19 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, COMPLEX ONE PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) * .. +* .. Local Scalars .. + REAL BIGNUM, BNRM + INTEGER I, J +* .. +* .. External Functions .. + REAL SLAMCH + EXTERNAL SLAMCH +* .. * .. External Subroutines .. - EXTERNAL CTRSM, CUNMQL, XERBLA + EXTERNAL CLASCL, CTRSM, CUNMQL, XERBLA * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, AIMAG, MAX, REAL * .. * .. Executable Statements .. * @@ -174,6 +183,23 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* Use component magnitudes since ABS of a finite complex entry +* can exceed the overflow threshold. +* + BNRM = 0.0E+0 + DO 20 J = 1, NRHS + DO 10 I = 1, M + BNRM = MAX( BNRM, ABS( REAL( B( I, J ) ) ), + $ ABS( AIMAG( B( I, J ) ) ) ) + 10 CONTINUE + 20 CONTINUE + BIGNUM = SLAMCH( 'Precision' ) / SLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL CUNMQL( 'Left', 'Conjugate transpose', M, NRHS, N, A, LDA, @@ -183,6 +209,10 @@ SUBROUTINE CGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL CTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/dgeqls.f b/TESTING/LIN/dgeqls.f index 749acb1ea..2061e14ec 100644 --- a/TESTING/LIN/dgeqls.f +++ b/TESTING/LIN/dgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by DGEQLF. +*> computed by DGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,8 +140,18 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. +* .. Local Scalars .. + DOUBLE PRECISION BIGNUM, BNRM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL DLAMCH, DLANGE +* .. * .. External Subroutines .. - EXTERNAL DORMQL, DTRSM, XERBLA + EXTERNAL DLASCL, DORMQL, DTRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -174,6 +185,15 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* + BNRM = DLANGE( 'M', M, NRHS, B, LDB, RWORK ) + BIGNUM = DLAMCH( 'Precision' ) / DLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL DORMQL( 'Left', 'Transpose', M, NRHS, N, A, LDA, TAU, B, LDB, @@ -183,6 +203,10 @@ SUBROUTINE DGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/sgeqls.f b/TESTING/LIN/sgeqls.f index e89fb37b7..1baf89239 100644 --- a/TESTING/LIN/sgeqls.f +++ b/TESTING/LIN/sgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by SGEQLF. +*> computed by SGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,8 +140,18 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, REAL ONE PARAMETER ( ONE = 1.0E+0 ) * .. +* .. Local Scalars .. + REAL BIGNUM, BNRM +* .. +* .. Local Arrays .. + REAL RWORK( 1 ) +* .. +* .. External Functions .. + REAL SLAMCH, SLANGE + EXTERNAL SLAMCH, SLANGE +* .. * .. External Subroutines .. - EXTERNAL SORMQL, STRSM, XERBLA + EXTERNAL SLASCL, SORMQL, STRSM, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -174,6 +185,15 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* + BNRM = SLANGE( 'M', M, NRHS, B, LDB, RWORK ) + BIGNUM = SLAMCH( 'Precision' ) / SLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL SORMQL( 'Left', 'Transpose', M, NRHS, N, A, LDA, TAU, B, LDB, @@ -183,6 +203,10 @@ SUBROUTINE SGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL STRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * diff --git a/TESTING/LIN/zgeqls.f b/TESTING/LIN/zgeqls.f index 5d0eb55e8..aa92bb044 100644 --- a/TESTING/LIN/zgeqls.f +++ b/TESTING/LIN/zgeqls.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the QL factorization *> A = Q*L -*> computed by ZGEQLF. +*> computed by ZGEQLF. Large right hand sides are scaled before +*> applying Q, and the solution is returned at the original scale. *> \endverbatim * * Arguments: @@ -139,11 +140,19 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) * .. +* .. Local Scalars .. + DOUBLE PRECISION BIGNUM, BNRM + INTEGER I, J +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH +* .. * .. External Subroutines .. - EXTERNAL XERBLA, ZTRSM, ZUNMQL + EXTERNAL XERBLA, ZLASCL, ZTRSM, ZUNMQL * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC ABS, DBLE, DIMAG, MAX * .. * .. Executable Statements .. * @@ -174,6 +183,23 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, IF( N.EQ.0 .OR. NRHS.EQ.0 .OR. M.EQ.0 ) $ RETURN * +* Scale large right hand sides before applying Q: tau*w can +* overflow even when Q' * B and the solution are representable. +* Use component magnitudes since ABS of a finite complex entry +* can exceed the overflow threshold. +* + BNRM = 0.0D+0 + DO 20 J = 1, NRHS + DO 10 I = 1, M + BNRM = MAX( BNRM, ABS( DBLE( B( I, J ) ) ), + $ ABS( DIMAG( B( I, J ) ) ) ) + 10 CONTINUE + 20 CONTINUE + BIGNUM = DLAMCH( 'Precision' ) / DLAMCH( 'Safe minimum' ) + IF( BNRM.GT.BIGNUM ) + $ CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, + $ INFO ) +* * B := Q' * B * CALL ZUNMQL( 'Left', 'Conjugate transpose', M, NRHS, N, A, LDA, @@ -183,6 +209,10 @@ SUBROUTINE ZGEQLS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N, NRHS, $ ONE, A( M-N+1, 1 ), LDA, B( M-N+1, 1 ), LDB ) +* + IF( BNRM.GT.BIGNUM ) + $ CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, M, NRHS, B, LDB, + $ INFO ) * RETURN * From 5aa9f46d95e0c6414288b47b204436f0066ee7c7 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Tue, 15 Sep 2026 20:13:45 -0700 Subject: [PATCH 3/3] TESTING: Protect complex RQ solves and reject nonfinite residuals --- TESTING/LIN/CMakeLists.txt | 12 +++++- TESTING/LIN/cgerqs.f | 30 ++++++++++---- TESTING/LIN/cget02.f | 13 ++++-- TESTING/LIN/ctest_rq_helpers.f90 | 68 ++++++++++++++++++++++++++++++++ TESTING/LIN/zgerqs.f | 30 ++++++++++---- TESTING/LIN/zget02.f | 13 ++++-- TESTING/LIN/ztest_rq_helpers.f90 | 68 ++++++++++++++++++++++++++++++++ 7 files changed, 211 insertions(+), 23 deletions(-) create mode 100644 TESTING/LIN/ctest_rq_helpers.f90 create mode 100644 TESTING/LIN/ztest_rq_helpers.f90 diff --git a/TESTING/LIN/CMakeLists.txt b/TESTING/LIN/CMakeLists.txt index 2313fa0c4..ec9bb134a 100644 --- a/TESTING/LIN/CMakeLists.txt +++ b/TESTING/LIN/CMakeLists.txt @@ -257,7 +257,7 @@ set(ZLINTSTRFP chkxer.f xerbla.f alaerh.f aladhd.f alahd.f alasvm.f) function(add_lin_executable name) - cmake_parse_arguments(PARSE_ARGV 1 LIN "" "" "SOURCES;DEFAULT_API;EXT_API") + cmake_parse_arguments(PARSE_ARGV 1 LIN "TEST" "" "SOURCES;DEFAULT_API;EXT_API") set(sources ${LIN_SOURCES} ${LIN_DEFAULT_API}) set(ext_sources ${LIN_SOURCES} ${LIN_EXT_API}) @@ -265,6 +265,9 @@ function(add_lin_executable name) add_executable(${name} ${sources}) target_link_libraries(${name} PRIVATE ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) lapack_add_coverage(${name}) + if(LIN_TEST) + add_test(NAME LAPACK-${name} COMMAND $) + endif() endif() if(BUILD_INDEX64_EXT_API) @@ -276,6 +279,9 @@ function(add_lin_executable name) target_compile_options(${name}_64 PRIVATE ${FOPT_ILP64}) target_link_libraries(${name}_64 PRIVATE ${TMGLIB} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}) lapack_add_coverage(${name}_64) + if(LIN_TEST) + add_test(NAME LAPACK-${name}_64 COMMAND $) + endif() # Add depedency to the global codegen target. Since we reuse the same # generated source file for multiple tests, the generation could be @@ -310,6 +316,8 @@ if(BUILD_COMPLEX) DEFAULT_API ${CXLINTST} EXT_API ${CLINTST_NO_XBLAS}) add_lin_executable(xlintstrfc SOURCES ${CLINTSTRFP}) + add_lin_executable(xrq_helpers_c TEST + SOURCES ctest_rq_helpers.f90 cgerqs.f cget02.f) endif() if(BUILD_COMPLEX16) @@ -318,6 +326,8 @@ if(BUILD_COMPLEX16) DEFAULT_API ${ZXLINTST} EXT_API ${ZLINTST_NO_XBLAS}) add_lin_executable(xlintstrfz SOURCES ${ZLINTSTRFP}) + add_lin_executable(xrq_helpers_z TEST + SOURCES ztest_rq_helpers.f90 zgerqs.f zget02.f) endif() if(BUILD_COMPLEX AND BUILD_COMPLEX16) diff --git a/TESTING/LIN/cgerqs.f b/TESTING/LIN/cgerqs.f index ac28a3643..0a48f6252 100644 --- a/TESTING/LIN/cgerqs.f +++ b/TESTING/LIN/cgerqs.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the RQ factorization *> A = R*Q -*> computed by CGERQF. +*> computed by CGERQF. The triangular solve uses scaled +*> complex division to avoid intermediate overflow. *> \endverbatim * * Arguments: @@ -136,12 +137,18 @@ SUBROUTINE CGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * ===================================================================== * * .. Parameters .. - COMPLEX CZERO, CONE - PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), - $ CONE = ( 1.0E+0, 0.0E+0 ) ) + COMPLEX CZERO + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J +* .. +* .. External Functions .. + COMPLEX CLADIV + EXTERNAL CLADIV * .. * .. External Subroutines .. - EXTERNAL CLASET, CTRSM, CUNMRQ, XERBLA + EXTERNAL CAXPY, CLASET, CUNMRQ, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -177,8 +184,17 @@ SUBROUTINE CGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * * Solve R*X = B(n-m+1:n,:) * - CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', M, NRHS, - $ CONE, A( 1, N-M+1 ), LDA, B( N-M+1, 1 ), LDB ) +* LADIV avoids overflow in intrinsic complex division. Scaling +* all of B instead could erase small solution components. +* + DO 20 J = 1, NRHS + DO 10 I = M, 1, -1 + B( N-M+I, J ) = CLADIV( B( N-M+I, J ), + $ A( I, N-M+I ) ) + CALL CAXPY( I-1, -B( N-M+I, J ), A( 1, N-M+I ), + $ 1, B( N-M+1, J ), 1 ) + 10 CONTINUE + 20 CONTINUE * * Set B(1:n-m,:) to zero * diff --git a/TESTING/LIN/cget02.f b/TESTING/LIN/cget02.f index cb92c3553..3e4c7e16a 100644 --- a/TESTING/LIN/cget02.f +++ b/TESTING/LIN/cget02.f @@ -168,7 +168,7 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, EXTERNAL CGEMM * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC HUGE, MAX * .. * .. Executable Statements .. * @@ -187,7 +187,7 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, N2 = N END IF * -* Exit with RESID = 1/EPS if ANORM = 0. +* Exit with RESID = 1/EPS if ANORM is zero or nonfinite. * EPS = SLAMCH( 'Epsilon' ) IF( LSAME( TRANS, 'N' ) ) THEN @@ -195,7 +195,8 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, ELSE ANORM = CLANGE( 'I', M, N, A, LDA, RWORK ) END IF - IF( ANORM.LE.ZERO ) THEN + IF( .NOT.( ANORM.GT.ZERO .AND. + $ ANORM.LE.HUGE( ANORM ) ) ) THEN RESID = ONE / EPS RETURN END IF @@ -212,8 +213,12 @@ SUBROUTINE CGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, DO 10 J = 1, NRHS BNORM = SCASUM( N1, B( 1, J ), 1 ) XNORM = SCASUM( N2, X( 1, J ), 1 ) - IF( XNORM.LE.ZERO ) THEN +* Reject nonfinite norms before MAX can hide a NaN ratio. + IF( .NOT.( XNORM.GT.ZERO .AND. + $ XNORM.LE.HUGE( XNORM ) .AND. + $ BNORM.LE.HUGE( BNORM ) ) ) THEN RESID = ONE / EPS + RETURN ELSE RESID = MAX( RESID, ( ( BNORM/ANORM )/XNORM )/EPS ) END IF diff --git a/TESTING/LIN/ctest_rq_helpers.f90 b/TESTING/LIN/ctest_rq_helpers.f90 new file mode 100644 index 000000000..ca77488cb --- /dev/null +++ b/TESTING/LIN/ctest_rq_helpers.f90 @@ -0,0 +1,68 @@ +! SPDX-FileCopyrightText: The LAPACK Authors +! SPDX-License-Identifier: BSD-3-Clause +! +! Regression tests for the RQ solve and residual test helpers. +program ctest_rq_helpers + use, intrinsic :: ieee_arithmetic + implicit none + integer, parameter :: wp = kind(1.0) + complex(wp) :: a(2,3), b(3,2), x(3,2), tau(2), work(128) + real(wp) :: large, small, resid, rwork(3), bad(2) + integer :: info, i, j, k + external :: cgerqs, cget02 + character :: trans(3) = ['N', 'T', 'C'] + + ! R has both ends of the normal range in one right-hand side. + ! Uniform scaling of B would lose the first solution component. + large = 0.75_wp*huge(1.0_wp) + small = tiny(1.0_wp) + a = 0 + a(1,2) = small + a(2,3) = large + tau = 0 + b = 0 + b(2,1) = small + b(3,1) = cmplx(-0.95_wp*large, 0.95_wp*large, wp) + b(2,2) = 2*small + b(3,2) = 0.5_wp*large + x = 0 + x(2,1) = 1 + x(3,1) = cmplx(-0.95_wp, 0.95_wp, wp) + x(2,2) = 2 + x(3,2) = 0.5_wp + call cgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 1 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 2 + + ! Check ordinary magnitudes as well. + a(1,2) = 2 + a(2,3) = 4 + b = 0 + b(2,:) = 2*x(2,:) + b(3,:) = 4*x(3,:) + call cgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 3 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 4 + + bad = [ieee_value(0.0_wp, ieee_quiet_nan), & + ieee_value(0.0_wp, ieee_positive_inf)] + do k = 1, 3 + do j = 1, 3 + do i = 1, 2 + a = 1 + x = 1 + b = 1 + if (j == 1) a(1,1) = bad(i) + if (j == 2) x(1,1) = bad(i) + if (j == 3) b(1,1) = bad(i) + call cget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (.not. (resid > 30)) stop 5 + end do + end do + a = 1 + x = 1 + b = 1 + call cget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (resid /= 0) stop 6 + end do +end program diff --git a/TESTING/LIN/zgerqs.f b/TESTING/LIN/zgerqs.f index ce3021284..2ff05346e 100644 --- a/TESTING/LIN/zgerqs.f +++ b/TESTING/LIN/zgerqs.f @@ -29,7 +29,8 @@ *> min || A*X - B || *> using the RQ factorization *> A = R*Q -*> computed by ZGERQF. +*> computed by ZGERQF. The triangular solve uses scaled +*> complex division to avoid intermediate overflow. *> \endverbatim * * Arguments: @@ -136,12 +137,18 @@ SUBROUTINE ZGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * ===================================================================== * * .. Parameters .. - COMPLEX*16 CZERO, CONE - PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), - $ CONE = ( 1.0D+0, 0.0D+0 ) ) + COMPLEX*16 CZERO + PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J +* .. +* .. External Functions .. + COMPLEX*16 ZLADIV + EXTERNAL ZLADIV * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLASET, ZTRSM, ZUNMRQ + EXTERNAL XERBLA, ZAXPY, ZLASET, ZUNMRQ * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -177,8 +184,17 @@ SUBROUTINE ZGERQS( M, N, NRHS, A, LDA, TAU, B, LDB, WORK, LWORK, * * Solve R*X = B(n-m+1:n,:) * - CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', M, NRHS, - $ CONE, A( 1, N-M+1 ), LDA, B( N-M+1, 1 ), LDB ) +* LADIV avoids overflow in intrinsic complex division. Scaling +* all of B instead could erase small solution components. +* + DO 20 J = 1, NRHS + DO 10 I = M, 1, -1 + B( N-M+I, J ) = ZLADIV( B( N-M+I, J ), + $ A( I, N-M+I ) ) + CALL ZAXPY( I-1, -B( N-M+I, J ), A( 1, N-M+I ), + $ 1, B( N-M+1, J ), 1 ) + 10 CONTINUE + 20 CONTINUE * * Set B(1:n-m,:) to zero * diff --git a/TESTING/LIN/zget02.f b/TESTING/LIN/zget02.f index 8f45177b5..ffe41aba0 100644 --- a/TESTING/LIN/zget02.f +++ b/TESTING/LIN/zget02.f @@ -168,7 +168,7 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, EXTERNAL ZGEMM * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC HUGE, MAX * .. * .. Executable Statements .. * @@ -187,7 +187,7 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, N2 = N END IF * -* Exit with RESID = 1/EPS if ANORM = 0. +* Exit with RESID = 1/EPS if ANORM is zero or nonfinite. * EPS = DLAMCH( 'Epsilon' ) IF( LSAME( TRANS, 'N' ) ) THEN @@ -195,7 +195,8 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, ELSE ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK ) END IF - IF( ANORM.LE.ZERO ) THEN + IF( .NOT.( ANORM.GT.ZERO .AND. + $ ANORM.LE.HUGE( ANORM ) ) ) THEN RESID = ONE / EPS RETURN END IF @@ -212,8 +213,12 @@ SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, DO 10 J = 1, NRHS BNORM = DZASUM( N1, B( 1, J ), 1 ) XNORM = DZASUM( N2, X( 1, J ), 1 ) - IF( XNORM.LE.ZERO ) THEN +* Reject nonfinite norms before MAX can hide a NaN ratio. + IF( .NOT.( XNORM.GT.ZERO .AND. + $ XNORM.LE.HUGE( XNORM ) .AND. + $ BNORM.LE.HUGE( BNORM ) ) ) THEN RESID = ONE / EPS + RETURN ELSE RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) END IF diff --git a/TESTING/LIN/ztest_rq_helpers.f90 b/TESTING/LIN/ztest_rq_helpers.f90 new file mode 100644 index 000000000..2030c8adb --- /dev/null +++ b/TESTING/LIN/ztest_rq_helpers.f90 @@ -0,0 +1,68 @@ +! SPDX-FileCopyrightText: The LAPACK Authors +! SPDX-License-Identifier: BSD-3-Clause +! +! Regression tests for the RQ solve and residual test helpers. +program ztest_rq_helpers + use, intrinsic :: ieee_arithmetic + implicit none + integer, parameter :: wp = kind(1.0d0) + complex(wp) :: a(2,3), b(3,2), x(3,2), tau(2), work(128) + real(wp) :: large, small, resid, rwork(3), bad(2) + integer :: info, i, j, k + external :: zgerqs, zget02 + character :: trans(3) = ['N', 'T', 'C'] + + ! R has both ends of the normal range in one right-hand side. + ! Uniform scaling of B would lose the first solution component. + large = 0.75_wp*huge(1.0_wp) + small = tiny(1.0_wp) + a = 0 + a(1,2) = small + a(2,3) = large + tau = 0 + b = 0 + b(2,1) = small + b(3,1) = cmplx(-0.95_wp*large, 0.95_wp*large, wp) + b(2,2) = 2*small + b(3,2) = 0.5_wp*large + x = 0 + x(2,1) = 1 + x(3,1) = cmplx(-0.95_wp, 0.95_wp, wp) + x(2,2) = 2 + x(3,2) = 0.5_wp + call zgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 1 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 2 + + ! Check ordinary magnitudes as well. + a(1,2) = 2 + a(2,3) = 4 + b = 0 + b(2,:) = 2*x(2,:) + b(3,:) = 4*x(3,:) + call zgerqs(2, 3, 2, a, 2, tau, b, 3, work, 128, info) + if (info /= 0) stop 3 + if (.not. all(abs(b-x) <= 8*epsilon(1.0_wp))) stop 4 + + bad = [ieee_value(0.0_wp, ieee_quiet_nan), & + ieee_value(0.0_wp, ieee_positive_inf)] + do k = 1, 3 + do j = 1, 3 + do i = 1, 2 + a = 1 + x = 1 + b = 1 + if (j == 1) a(1,1) = bad(i) + if (j == 2) x(1,1) = bad(i) + if (j == 3) b(1,1) = bad(i) + call zget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (.not. (resid > 30)) stop 5 + end do + end do + a = 1 + x = 1 + b = 1 + call zget02(trans(k), 1, 1, 1, a, 2, x, 3, b, 3, rwork, resid) + if (resid /= 0) stop 6 + end do +end program