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