From 8ad354a807fafe2d299341dd0c54073926d3dbad Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 21:06:15 -0700 Subject: [PATCH] Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal xGBTF2 and xGBTRF form the multipliers of column J as CALL xSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 ) When the pivot is subnormal its reciprocal is not representable, so the whole column of L becomes Inf, or NaN where an entry is zero, the xGER update spreads that through the band, and the factorization completes with INFO = 0. Partial pivoting cannot avoid it: the pivot is the largest entry of the column, and when the matrix is small every column is. A well-conditioned banded system scaled to 2^-1030 is solved by xGESV to 1e-14 and by xGBSV to NaN. The dense routines xGETF2 and xGETRF2 have guarded this since at least LAPACK 3.2 by dividing element-wise when |pivot| < SFMIN. Apply the same test at both banded sites, in all four precisions. SFMIN comes from xLAMCH('S'), computed once after the quick return as in xGETF2. The GB test path gets a matrix type for it: type 9 is the type 1 matrix scaled into the subnormal range, one eighth of the safe minimum, which xLATMS cannot generate because it scales its output to the requested norm. At that scale the matrix carries too few bits to reconstruct, so the type tests what the guarded division promises, that the factor is finite, instead of a residual; the remaining ratios, which estimate a condition number from a subnormal norm, are skipped as for a block size other than the first. On the parent commit the type fails 798 times per precision. The new branch is taken only when |pivot| < SFMIN. Over a sweep of 2124 (precision, n, KL, KU, scale) cases the banded factor and IPIV are bit-identical to the parent commit for every pivot of normal magnitude, and for the subnormal cases the banded solution error now equals the dense one in every case, where before 560 of 708 came back NaN or Inf. The full LAPACK test suite passes: 5447193 LAPACK and 315872 BLAS tests, 0 numerical errors, 0 other errors; the 5292 tests above the parent are the new type. The test files declare the new xSCAL calls EXTERNAL: the extended-API build renames only the routines a file declares, so without the declaration the xlintst*_64 executables failed to link against the 64-bit BLAS. Co-Authored-By: Claude Fable 5.1 --- SRC/cgbtf2.f | 18 +++++++++++++--- SRC/cgbtrf.f | 20 +++++++++++++---- SRC/dgbtf2.f | 18 +++++++++++++--- SRC/dgbtrf.f | 21 +++++++++++++----- SRC/sgbtf2.f | 18 +++++++++++++--- SRC/sgbtrf.f | 21 +++++++++++++----- SRC/zgbtf2.f | 18 +++++++++++++--- SRC/zgbtrf.f | 20 +++++++++++++---- TESTING/LIN/alahd.f | 3 ++- TESTING/LIN/cchkaa.F | 2 +- TESTING/LIN/cchkgb.f | 45 ++++++++++++++++++++++++++++++++------ TESTING/LIN/dchkaa.F | 2 +- TESTING/LIN/dchkgb.f | 45 ++++++++++++++++++++++++++++++++------ TESTING/LIN/schkaa.F | 2 +- TESTING/LIN/schkgb.f | 45 ++++++++++++++++++++++++++++++++------ TESTING/LIN/zchkaa.F | 2 +- TESTING/LIN/zchkgb.f | 51 ++++++++++++++++++++++++++++++++++++-------- TESTING/ctest.in | 2 +- TESTING/dtest.in | 2 +- TESTING/stest.in | 2 +- TESTING/ztest.in | 2 +- 21 files changed, 293 insertions(+), 66 deletions(-) diff --git a/SRC/cgbtf2.f b/SRC/cgbtf2.f index b93a1c398..fdcfdb178 100644 --- a/SRC/cgbtf2.f +++ b/SRC/cgbtf2.f @@ -162,11 +162,13 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) $ ZERO = ( 0.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. + REAL SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + REAL SLAMCH INTEGER ICAMAX - EXTERNAL ICAMAX + EXTERNAL SLAMCH, ICAMAX * .. * .. External Subroutines .. EXTERNAL CGERU, CSCAL, CSWAP, XERBLA @@ -205,6 +207,10 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -248,8 +254,14 @@ SUBROUTINE CGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), - $ 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL CSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), + $ 1 ) + ELSE + DO 35 I = 1, KM + AB( KV+1+I, J ) = AB( KV+1+I, J ) / AB( KV+1, J ) + 35 CONTINUE + END IF * * Update trailing submatrix within the band. * diff --git a/SRC/cgbtrf.f b/SRC/cgbtrf.f index 4f0da00cc..eb00b7fb1 100644 --- a/SRC/cgbtrf.f +++ b/SRC/cgbtrf.f @@ -166,14 +166,16 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW COMPLEX TEMP + REAL SFMIN * .. * .. Local Arrays .. COMPLEX WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + REAL SLAMCH INTEGER ICAMAX, ILAENV - EXTERNAL ICAMAX, ILAENV + EXTERNAL SLAMCH, ICAMAX, ILAENV * .. * .. External Subroutines .. EXTERNAL CCOPY, CGBTF2, CGEMM, CGERU, CLASWP, @@ -214,6 +216,10 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'CGBTRF', ' ', M, N, KL, KU ) @@ -325,9 +331,15 @@ SUBROUTINE CGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL CSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, - $ JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL CSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column diff --git a/SRC/dgbtf2.f b/SRC/dgbtf2.f index 8ff949b96..6b1991e76 100644 --- a/SRC/dgbtf2.f +++ b/SRC/dgbtf2.f @@ -161,11 +161,13 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. + DOUBLE PRECISION SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IDAMAX - EXTERNAL IDAMAX + EXTERNAL DLAMCH, IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA @@ -204,6 +206,10 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -248,8 +254,14 @@ SUBROUTINE DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), - $ 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), + $ 1 ) + ELSE + DO 35 I = 1, KM + AB( KV+1+I, J ) = AB( KV+1+I, J ) / AB( KV+1, J ) + 35 CONTINUE + END IF * * Update trailing submatrix within the band. * diff --git a/SRC/dgbtrf.f b/SRC/dgbtrf.f index 6dd8d9008..8bee3120c 100644 --- a/SRC/dgbtrf.f +++ b/SRC/dgbtrf.f @@ -164,15 +164,16 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW - DOUBLE PRECISION TEMP + DOUBLE PRECISION SFMIN, TEMP * .. * .. Local Arrays .. DOUBLE PRECISION WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IDAMAX, ILAENV - EXTERNAL IDAMAX, ILAENV + EXTERNAL DLAMCH, IDAMAX, ILAENV * .. * .. External Subroutines .. EXTERNAL DCOPY, DGBTF2, DGEMM, DGER, DLASWP, @@ -213,6 +214,10 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU ) @@ -324,9 +329,15 @@ SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, - $ JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL DSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column diff --git a/SRC/sgbtf2.f b/SRC/sgbtf2.f index d8e3918f8..53026dd26 100644 --- a/SRC/sgbtf2.f +++ b/SRC/sgbtf2.f @@ -161,11 +161,13 @@ SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. * .. Local Scalars .. + REAL SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + REAL SLAMCH INTEGER ISAMAX - EXTERNAL ISAMAX + EXTERNAL SLAMCH, ISAMAX * .. * .. External Subroutines .. EXTERNAL SGER, SSCAL, SSWAP, XERBLA @@ -204,6 +206,10 @@ SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -248,8 +254,14 @@ SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL SSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), - $ 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL SSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), + $ 1 ) + ELSE + DO 35 I = 1, KM + AB( KV+1+I, J ) = AB( KV+1+I, J ) / AB( KV+1, J ) + 35 CONTINUE + END IF * * Update trailing submatrix within the band. * diff --git a/SRC/sgbtrf.f b/SRC/sgbtrf.f index 322f4d6db..c9ac5327b 100644 --- a/SRC/sgbtrf.f +++ b/SRC/sgbtrf.f @@ -164,15 +164,16 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * .. Local Scalars .. INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW - REAL TEMP + REAL SFMIN, TEMP * .. * .. Local Arrays .. REAL WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + REAL SLAMCH INTEGER ILAENV, ISAMAX - EXTERNAL ILAENV, ISAMAX + EXTERNAL SLAMCH, ILAENV, ISAMAX * .. * .. External Subroutines .. EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, @@ -213,6 +214,10 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = SLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU ) @@ -324,9 +329,15 @@ SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, - $ JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL SSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column diff --git a/SRC/zgbtf2.f b/SRC/zgbtf2.f index 90f83e19a..b8f200910 100644 --- a/SRC/zgbtf2.f +++ b/SRC/zgbtf2.f @@ -162,11 +162,13 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) $ ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Local Scalars .. + DOUBLE PRECISION SFMIN INTEGER I, J, JP, JU, KM, KV * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER IZAMAX - EXTERNAL IZAMAX + EXTERNAL DLAMCH, IZAMAX * .. * .. External Subroutines .. EXTERNAL XERBLA, ZGERU, ZSCAL, ZSWAP @@ -205,6 +207,10 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Gaussian elimination with partial pivoting * * Set fill-in elements in columns KU+2 to KV to zero. @@ -248,8 +254,14 @@ SUBROUTINE ZGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers. * - CALL ZSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), - $ 1 ) + IF( ABS( AB( KV+1, J ) ).GE.SFMIN ) THEN + CALL ZSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), + $ 1 ) + ELSE + DO 35 I = 1, KM + AB( KV+1+I, J ) = AB( KV+1+I, J ) / AB( KV+1, J ) + 35 CONTINUE + END IF * * Update trailing submatrix within the band. * diff --git a/SRC/zgbtrf.f b/SRC/zgbtrf.f index a642ab094..9c8f55f7f 100644 --- a/SRC/zgbtrf.f +++ b/SRC/zgbtrf.f @@ -166,14 +166,16 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP, $ JU, K2, KM, KV, NB, NW COMPLEX*16 TEMP + DOUBLE PRECISION SFMIN * .. * .. Local Arrays .. COMPLEX*16 WORK13( LDWORK, NBMAX ), $ WORK31( LDWORK, NBMAX ) * .. * .. External Functions .. + DOUBLE PRECISION DLAMCH INTEGER ILAENV, IZAMAX - EXTERNAL ILAENV, IZAMAX + EXTERNAL DLAMCH, ILAENV, IZAMAX * .. * .. External Subroutines .. EXTERNAL XERBLA, ZCOPY, ZGBTF2, ZGEMM, ZGERU, @@ -214,6 +216,10 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * +* Compute machine safe minimum +* + SFMIN = DLAMCH('S') +* * Determine the block size for this environment * NB = ILAENV( 1, 'ZGBTRF', ' ', M, N, KL, KU ) @@ -325,9 +331,15 @@ SUBROUTINE ZGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO ) * * Compute multipliers * - CALL ZSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, - $ JJ ), - $ 1 ) + IF( ABS( AB( KV+1, JJ ) ).GE.SFMIN ) THEN + CALL ZSCAL( KM, ONE / AB( KV+1, JJ ), + $ AB( KV+2, JJ ), 1 ) + ELSE + DO 75 I = 1, KM + AB( KV+1+I, JJ ) = AB( KV+1+I, JJ ) / + $ AB( KV+1, JJ ) + 75 CONTINUE + END IF * * Update trailing submatrix within the band and within * the current block. JM is the index of the last column diff --git a/TESTING/LIN/alahd.f b/TESTING/LIN/alahd.f index b04a3f796..0039efdb3 100644 --- a/TESTING/LIN/alahd.f +++ b/TESTING/LIN/alahd.f @@ -831,7 +831,8 @@ SUBROUTINE ALAHD( IOUNIT, PATH ) $ '2. First column zero', 15X, '6. Random, CNDNUM = .01/EPS', $ / 4X, '3. Last column zero', 16X, $ '7. Scaled near underflow', / 4X, - $ '4. Last n/2 columns zero', 11X, '8. Scaled near overflow' ) + $ '4. Last n/2 columns zero', 11X, '8. Scaled near overflow', + $ / 39X, '9. Scaled into the subnormal range' ) * * GT matrix types * diff --git a/TESTING/LIN/cchkaa.F b/TESTING/LIN/cchkaa.F index fa4000050..3430b0475 100644 --- a/TESTING/LIN/cchkaa.F +++ b/TESTING/LIN/cchkaa.F @@ -490,7 +490,7 @@ PROGRAM CCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/cchkgb.f b/TESTING/LIN/cchkgb.f index 17e4a495e..c93abfa0e 100644 --- a/TESTING/LIN/cchkgb.f +++ b/TESTING/LIN/cchkgb.f @@ -214,7 +214,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -228,6 +228,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + REAL ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -236,14 +237,17 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH REAL CLANGB, CLANGE, SGET06 EXTERNAL CLANGB, CLANGE, SGET06 + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, CCOPY, CERRGE, CGBCON, $ CGBRFS, CGBT01, CGBT02, CGBT05, CGBTRF, CGBTRS, $ CGET04, CLACPY, CLARHS, CLASET, CLATB4, CLATMS, - $ XLAENV + $ CSSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC CMPLX, MAX, MIN @@ -401,6 +405,17 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. CLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = SLAMCH( 'Safe minimum' ) / 8 + CALL CSSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -478,13 +493,31 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL CGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = CLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( SISNAN( ANRMF ) .OR. + $ ANRMF.GT.SLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL CGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ SISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -496,7 +529,7 @@ SUBROUTINE CCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = CLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/TESTING/LIN/dchkaa.F b/TESTING/LIN/dchkaa.F index 856185c65..606c4f290 100644 --- a/TESTING/LIN/dchkaa.F +++ b/TESTING/LIN/dchkaa.F @@ -486,7 +486,7 @@ PROGRAM DCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/dchkgb.f b/TESTING/LIN/dchkgb.f index 80e2b9b8f..7bb238a6c 100644 --- a/TESTING/LIN/dchkgb.f +++ b/TESTING/LIN/dchkgb.f @@ -213,7 +213,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -227,6 +227,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN DOUBLE PRECISION AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + DOUBLE PRECISION ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -235,14 +236,17 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH DOUBLE PRECISION DGET06, DLANGB, DLANGE EXTERNAL DGET06, DLANGB, DLANGE + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, DCOPY, DERRGE, DGBCON, $ DGBRFS, DGBT01, DGBT02, DGBT05, DGBTRF, DGBTRS, $ DGET04, DLACPY, DLARHS, DLASET, DLATB4, DLATMS, - $ XLAENV + $ DSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -401,6 +405,17 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. DLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = DLAMCH( 'Safe minimum' ) / 8 + CALL DSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -478,13 +493,31 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL DGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = DLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( DISNAN( ANRMF ) .OR. + $ ANRMF.GT.DLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL DGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ DISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -496,7 +529,7 @@ SUBROUTINE DCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = DLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/TESTING/LIN/schkaa.F b/TESTING/LIN/schkaa.F index 25c455d4f..50630d0e1 100644 --- a/TESTING/LIN/schkaa.F +++ b/TESTING/LIN/schkaa.F @@ -482,7 +482,7 @@ PROGRAM SCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/schkgb.f b/TESTING/LIN/schkgb.f index 5711438d0..b02ff97d2 100644 --- a/TESTING/LIN/schkgb.f +++ b/TESTING/LIN/schkgb.f @@ -213,7 +213,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL ONE, ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -227,6 +227,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN REAL AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + REAL ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -235,14 +236,17 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, REAL RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL SISNAN + REAL SLAMCH REAL SGET06, SLANGB, SLANGE EXTERNAL SGET06, SLANGB, SLANGE + EXTERNAL SISNAN, SLAMCH * .. * .. External Subroutines .. EXTERNAL ALAERH, ALAHD, ALASUM, SCOPY, SERRGE, SGBCON, $ SGBRFS, SGBT01, SGBT02, SGBT05, SGBTRF, SGBTRS, $ SGET04, SLACPY, SLARHS, SLASET, SLATB4, SLATMS, - $ XLAENV + $ SSCAL, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -401,6 +405,17 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. SLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = SLAMCH( 'Safe minimum' ) / 8 + CALL SSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -478,13 +493,31 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = SLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( SISNAN( ANRMF ) .OR. + $ ANRMF.GT.SLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL SGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ SISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -496,7 +529,7 @@ SUBROUTINE SCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = SLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/TESTING/LIN/zchkaa.F b/TESTING/LIN/zchkaa.F index a86793907..04fce98a3 100644 --- a/TESTING/LIN/zchkaa.F +++ b/TESTING/LIN/zchkaa.F @@ -490,7 +490,7 @@ PROGRAM ZCHKAA * LA = ( 2*KDMAX+1 )*NMAX LAFAC = ( 3*KDMAX+1 )*NMAX - NTYPES = 8 + NTYPES = 9 CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT ) * IF( TSTCHK ) THEN diff --git a/TESTING/LIN/zchkgb.f b/TESTING/LIN/zchkgb.f index bd54ac5e8..b4a288362 100644 --- a/TESTING/LIN/zchkgb.f +++ b/TESTING/LIN/zchkgb.f @@ -214,7 +214,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) INTEGER NTYPES, NTESTS - PARAMETER ( NTYPES = 8, NTESTS = 7 ) + PARAMETER ( NTYPES = 9, NTESTS = 7 ) INTEGER NBW, NTRAN PARAMETER ( NBW = 4, NTRAN = 3 ) * .. @@ -228,6 +228,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NIMAT, NKL, NKU, NRHS, NRUN DOUBLE PRECISION AINVNM, ANORM, ANORMI, ANORMO, CNDNUM, RCOND, $ RCONDC, RCONDI, RCONDO + DOUBLE PRECISION ANRMF, SUBNRM * .. * .. Local Arrays .. CHARACTER TRANSS( NTRAN ) @@ -236,14 +237,17 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, DOUBLE PRECISION RESULT( NTESTS ) * .. * .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLAMCH DOUBLE PRECISION DGET06, ZLANGB, ZLANGE EXTERNAL DGET06, ZLANGB, ZLANGE + EXTERNAL DISNAN, DLAMCH * .. * .. External Subroutines .. - EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZCOPY, ZERRGE, - $ ZGBCON, ZGBRFS, ZGBT01, ZGBT02, ZGBT05, ZGBTRF, - $ ZGBTRS, ZGET04, ZLACPY, ZLARHS, ZLASET, ZLATB4, - $ ZLATMS + EXTERNAL ALAERH, ALAHD, ALASUM, XLAENV, ZCOPY, ZDSCAL, + $ ZERRGE, ZGBCON, ZGBRFS, ZGBT01, ZGBT02, ZGBT05, + $ ZGBTRF, ZGBTRS, ZGET04, ZLACPY, ZLARHS, ZLASET, + $ ZLATB4, ZLATMS * .. * .. Intrinsic Functions .. INTRINSIC DCMPLX, MAX, MIN @@ -401,6 +405,17 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, $ NERRS, NOUT ) GO TO 120 END IF +* +* Type 9: scale the matrix into the subnormal +* range, where the reciprocal of a pivot +* overflows. ZLATMS cannot generate such a +* matrix, because it scales its output by the +* requested norm. +* + IF( IMAT.EQ.9 ) THEN + SUBNRM = DLAMCH( 'Safe minimum' ) / 8 + CALL ZDSCAL( LDA*N, SUBNRM, A, 1 ) + END IF ELSE IF( IZERO.GT.0 ) THEN * * Use the same matrix for types 3 and 4 as for @@ -478,13 +493,31 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Reconstruct matrix from factors and compute * residual. * - CALL ZGBT01( M, N, KL, KU, A, LDA, AFAC, LDAFAC, - $ IWORK, WORK, RESULT( 1 ) ) + IF( IMAT.EQ.9 ) THEN +* +* The subnormal matrix carries no accuracy to +* reconstruct, so test what the guarded pivot +* division promises: a finite factor. +* + ANRMF = ZLANGB( 'M', N, KL, KL+KU, AFAC, + $ LDAFAC, RWORK ) + IF( DISNAN( ANRMF ) .OR. + $ ANRMF.GT.DLAMCH( 'Overflow' ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = ZERO + END IF + ELSE + CALL ZGBT01( M, N, KL, KU, A, LDA, AFAC, + $ LDAFAC, IWORK, WORK, + $ RESULT( 1 ) ) + END IF * * Print information about the tests so far that * did not pass the threshold. * - IF( RESULT( 1 ).GE.THRESH ) THEN + IF( RESULT( 1 ).GE.THRESH .OR. + $ DISNAN( RESULT( 1 ) ) ) THEN IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 ) $ CALL ALAHD( NOUT, PATH ) WRITE( NOUT, FMT = 9997 )M, N, KL, KU, NB, @@ -496,7 +529,7 @@ SUBROUTINE ZCHKGB( DOTYPE, NM, MVAL, NN, NVAL, NNB, NBVAL, NNS, * Skip the remaining tests if this is not the * first block size or if M .ne. N. * - IF( INB.GT.1 .OR. M.NE.N ) + IF( INB.GT.1 .OR. M.NE.N .OR. IMAT.EQ.9 ) $ GO TO 110 * ANORMO = ZLANGB( 'O', N, KL, KU, A, LDA, RWORK ) diff --git a/TESTING/ctest.in b/TESTING/ctest.in index 4e30224d7..d740caddf 100644 --- a/TESTING/ctest.in +++ b/TESTING/ctest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits CGE 11 List types on next line if 0 < NTYPES < 11 -CGB 8 List types on next line if 0 < NTYPES < 8 +CGB 9 List types on next line if 0 < NTYPES < 9 CGT 12 List types on next line if 0 < NTYPES < 12 CPO 9 List types on next line if 0 < NTYPES < 9 CPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/TESTING/dtest.in b/TESTING/dtest.in index cde62db50..40db81d78 100644 --- a/TESTING/dtest.in +++ b/TESTING/dtest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits DGE 11 List types on next line if 0 < NTYPES < 11 -DGB 8 List types on next line if 0 < NTYPES < 8 +DGB 9 List types on next line if 0 < NTYPES < 9 DGT 12 List types on next line if 0 < NTYPES < 12 DPO 9 List types on next line if 0 < NTYPES < 9 DPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/TESTING/stest.in b/TESTING/stest.in index abfd639fd..1bc94b348 100644 --- a/TESTING/stest.in +++ b/TESTING/stest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits SGE 11 List types on next line if 0 < NTYPES < 11 -SGB 8 List types on next line if 0 < NTYPES < 8 +SGB 9 List types on next line if 0 < NTYPES < 9 SGT 12 List types on next line if 0 < NTYPES < 12 SPO 9 List types on next line if 0 < NTYPES < 9 SPS 9 List types on next line if 0 < NTYPES < 9 diff --git a/TESTING/ztest.in b/TESTING/ztest.in index bf4c9d100..eeb1258b7 100644 --- a/TESTING/ztest.in +++ b/TESTING/ztest.in @@ -15,7 +15,7 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits ZGE 11 List types on next line if 0 < NTYPES < 11 -ZGB 8 List types on next line if 0 < NTYPES < 8 +ZGB 9 List types on next line if 0 < NTYPES < 9 ZGT 12 List types on next line if 0 < NTYPES < 12 ZPO 9 List types on next line if 0 < NTYPES < 9 ZPS 9 List types on next line if 0 < NTYPES < 9