Skip to content

Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal #1379

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rmlarsen wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Divide by a subnormal pivot in the banded LU factorizations instead of scaling by its reciprocal #1379

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 15 additions & 3 deletions SRC/cgbtf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
*
Expand Down
20 changes: 16 additions & 4 deletions SRC/cgbtrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand Down
18 changes: 15 additions & 3 deletions SRC/dgbtf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
*
Expand Down
21 changes: 16 additions & 5 deletions SRC/dgbtrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand Down
18 changes: 15 additions & 3 deletions SRC/sgbtf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
*
Expand Down
21 changes: 16 additions & 5 deletions SRC/sgbtrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand Down
18 changes: 15 additions & 3 deletions SRC/zgbtf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
*
Expand Down
20 changes: 16 additions & 4 deletions SRC/zgbtrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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 )
Expand Down Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion TESTING/LIN/alahd.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
*
Expand Down
2 changes: 1 addition & 1 deletion TESTING/LIN/cchkaa.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading