Skip to content

Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite factorizations #1377

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 4 commits into
base: master
Choose a base branch
from
Open

Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite factorizations #1377

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7bd0c3d
Avoid overflow in the 2x2 pivot inverse of the symmetric-indefinite f…
rmlarsen Sep 6, 2026
273ab29
LAPACK: Guard complex 2x2 pivot division
rmlarsen Sep 7, 2026
36b5001
Revert "LAPACK: Guard complex 2x2 pivot division"
rmlarsen Sep 7, 2026
6b2ba55
Add a subnormal 2x2 pivot block to the symmetric-indefinite tests
rmlarsen Sep 8, 2026
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
12 changes: 6 additions & 6 deletions SRC/chetf2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -403,11 +403,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
D11 = REAL( A( K, K ) ) / D
TT = ONE / ( D11*D22-ONE )
D12 = A( K-1, K ) / D
D = TT / D
*
DO 40 J = K - 2, 1, -1
WKM1 = D*( D11*A( J, K-1 )-CONJG( D12 )*A( J, K ) )
WK = D*( D22*A( J, K )-D12*A( J, K-1 ) )
WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )*
$ A( J, K ) ) / D )
WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D )
DO 30 I = J, 1, -1
A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
$ A( I, K-1 )*CONJG( WKM1 )
Expand Down Expand Up @@ -593,11 +593,11 @@ SUBROUTINE CHETF2( UPLO, N, A, LDA, IPIV, INFO )
D22 = REAL( A( K, K ) ) / D
TT = ONE / ( D11*D22-ONE )
D21 = A( K+1, K ) / D
D = TT / D
*
DO 80 J = K + 2, N
WK = D*( D11*A( J, K )-D21*A( J, K+1 ) )
WKP1 = D*( D22*A( J, K+1 )-CONJG( D21 )*A( J, K ) )
WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D )
WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )*
$ A( J, K ) ) / D )
DO 70 I = J, N
A( I, J ) = A( I, J ) - A( I, K )*CONJG( WK ) -
$ A( I, K+1 )*CONJG( WKP1 )
Expand Down
32 changes: 16 additions & 16 deletions SRC/chetf2_rk.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -617,24 +617,24 @@ SUBROUTINE CHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
*
DO 30 J = K - 2, 1, -1
*
* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J
*
WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )*
$ A( J, K ) )
WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) )
WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )*
$ A( J, K ) ) / D )
WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D )
*
* Perform a rank-2 update of A(1:k-2,1:k-2)
*
DO 20 I = J, 1, -1
A( I, J ) = A( I, J ) -
$ ( A( I, K ) / D )*CONJG( WK ) -
$ ( A( I, K-1 ) / D )*CONJG( WKM1 )
$ A( I, K )*CONJG( WK ) -
$ A( I, K-1 )*CONJG( WKM1 )
20 CONTINUE
*
* Store U(k) and U(k-1) in cols k and k-1 for row J
*
A( J, K ) = WK / D
A( J, K-1 ) = WKM1 / D
A( J, K ) = WK
A( J, K-1 ) = WKM1
* (*) Make sure that diagonal element of pivot is real
A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
*
Expand Down Expand Up @@ -978,24 +978,24 @@ SUBROUTINE CHETF2_RK( UPLO, N, A, LDA, E, IPIV, INFO )
*
DO 60 J = K + 2, N
*
* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J
*
WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) )
WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )*
$ A( J, K ) )
WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D )
WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )*
$ A( J, K ) ) / D )
*
* Perform a rank-2 update of A(k+2:n,k+2:n)
*
DO 50 I = J, N
A( I, J ) = A( I, J ) -
$ ( A( I, K ) / D )*CONJG( WK ) -
$ ( A( I, K+1 ) / D )*CONJG( WKP1 )
$ A( I, K )*CONJG( WK ) -
$ A( I, K+1 )*CONJG( WKP1 )
50 CONTINUE
*
* Store L(k) and L(k+1) in cols k and k+1 for row J
*
A( J, K ) = WK / D
A( J, K+1 ) = WKP1 / D
A( J, K ) = WK
A( J, K+1 ) = WKP1
* (*) Make sure that diagonal element of pivot is real
A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
*
Expand Down
32 changes: 16 additions & 16 deletions SRC/chetf2_rook.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -536,24 +536,24 @@ SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
*
DO 30 J = K - 2, 1, -1
*
* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J
*
WKM1 = TT*( D11*A( J, K-1 )-CONJG( D12 )*
$ A( J, K ) )
WK = TT*( D22*A( J, K )-D12*A( J, K-1 ) )
WKM1 = TT*( ( D11*A( J, K-1 )-CONJG( D12 )*
$ A( J, K ) ) / D )
WK = TT*( ( D22*A( J, K )-D12*A( J, K-1 ) ) / D )
*
* Perform a rank-2 update of A(1:k-2,1:k-2)
*
DO 20 I = J, 1, -1
A( I, J ) = A( I, J ) -
$ ( A( I, K ) / D )*CONJG( WK ) -
$ ( A( I, K-1 ) / D )*CONJG( WKM1 )
$ A( I, K )*CONJG( WK ) -
$ A( I, K-1 )*CONJG( WKM1 )
20 CONTINUE
*
* Store U(k) and U(k-1) in cols k and k-1 for row J
*
A( J, K ) = WK / D
A( J, K-1 ) = WKM1 / D
A( J, K ) = WK
A( J, K-1 ) = WKM1
* (*) Make sure that diagonal element of pivot is real
A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
*
Expand Down Expand Up @@ -857,24 +857,24 @@ SUBROUTINE CHETF2_ROOK( UPLO, N, A, LDA, IPIV, INFO )
*
DO 60 J = K + 2, N
*
* Compute D21 * ( W(k)W(k+1) ) * inv(D(k)) for row J
* Compute ( W(k)W(k+1) ) * inv(D(k)) for row J
*
WK = TT*( D11*A( J, K )-D21*A( J, K+1 ) )
WKP1 = TT*( D22*A( J, K+1 )-CONJG( D21 )*
$ A( J, K ) )
WK = TT*( ( D11*A( J, K )-D21*A( J, K+1 ) ) / D )
WKP1 = TT*( ( D22*A( J, K+1 )-CONJG( D21 )*
$ A( J, K ) ) / D )
*
* Perform a rank-2 update of A(k+2:n,k+2:n)
*
DO 50 I = J, N
A( I, J ) = A( I, J ) -
$ ( A( I, K ) / D )*CONJG( WK ) -
$ ( A( I, K+1 ) / D )*CONJG( WKP1 )
$ A( I, K )*CONJG( WK ) -
$ A( I, K+1 )*CONJG( WKP1 )
50 CONTINUE
*
* Store L(k) and L(k+1) in cols k and k+1 for row J
*
A( J, K ) = WK / D
A( J, K+1 ) = WKP1 / D
A( J, K ) = WK
A( J, K+1 ) = WKP1
* (*) Make sure that diagonal element of pivot is real
A( J, J ) = CMPLX( REAL( A( J, J ) ), ZERO )
*
Expand Down
19 changes: 9 additions & 10 deletions SRC/chptrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -389,13 +389,12 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO )
D11 = REAL( AP( K+( K-1 )*K / 2 ) ) / D
TT = ONE / ( D11*D22-ONE )
D12 = AP( K-1+( K-1 )*K / 2 ) / D
D = TT / D
*
DO 50 J = K - 2, 1, -1
WKM1 = D*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
$ CONJG( D12 )*AP( J+( K-1 )*K / 2 ) )
WK = D*( D22*AP( J+( K-1 )*K / 2 )-D12*
$ AP( J+( K-2 )*( K-1 ) / 2 ) )
WKM1 = TT*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
$ CONJG( D12 )*AP( J+( K-1 )*K / 2 ) ) / D )
WK = TT*( ( D22*AP( J+( K-1 )*K / 2 )-D12*
$ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D )
DO 40 I = J, 1, -1
AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
$ AP( I+( K-1 )*K / 2 )*CONJG( WK ) -
Expand Down Expand Up @@ -603,13 +602,13 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO )
D22 = REAL( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D
TT = ONE / ( D11*D22-ONE )
D21 = AP( K+1+( K-1 )*( 2*N-K ) / 2 ) / D
D = TT / D
*
DO 100 J = K + 2, N
WK = D*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21*
$ AP( J+K*( 2*N-K-1 ) / 2 ) )
WKP1 = D*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
$ CONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 2 ) )
WK = TT*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-D21*
$ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D )
WKP1 = TT*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
$ CONJG( D21 )*AP( J+( K-1 )*( 2*N-K ) / 2 )
$ ) / D )
DO 90 I = J, N
AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
$ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
Expand Down
36 changes: 18 additions & 18 deletions SRC/clahef.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -481,17 +481,17 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
* ( ( -1 ) ( D22 ) )
*
* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ),
* ( ( -1 ) ( D22 ) )
*
* = ( conj(D21)*( D11 ) D21*( -1 ) )
* ( ( -1 ) ( D22 ) ),
*
* where D11 = d22/d21,
* D22 = d11/conj(d21),
* D21 = T/d21,
* T = 1/(D22*D11-1).
*
* T/d21 is not formed, since it overflows when d21 is
* subnormal: each entry of the product is divided by
* d21 or conj(d21) and then scaled by T.
*
* (NOTE: No need to check for division by ZERO,
* since that was ensured earlier in pivot search:
* (a) d21 != 0, since in 2x2 pivot case(4)
Expand All @@ -503,16 +503,16 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
D11 = W( K, KW ) / CONJG( D21 )
D22 = W( K-1, KW-1 ) / D21
T = ONE / ( REAL( D11*D22 )-ONE )
D21 = T / D21
*
* Update elements in columns A(k-1) and A(k) as
* dot products of rows of ( W(kw-1) W(kw) ) and columns
* of D**(-1)
*
DO 20 J = 1, K - 2
A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
A( J, K ) = CONJG( D21 )*
$ ( D22*W( J, KW )-W( J, KW-1 ) )
A( J, K-1 ) = T*( ( D11*W( J, KW-1 )-W( J, KW ) ) /
$ D21 )
A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) /
$ CONJG( D21 ) )
20 CONTINUE
END IF
*
Expand Down Expand Up @@ -828,17 +828,17 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* = (1/|d21|**2) * T * ( d21*( D11 ) conj(d21)*( -1 ) ) =
* ( ( -1 ) ( D22 ) )
*
* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ) =
* = ( (T/conj(d21))*( D11 ) (T/d21)*( -1 ) ),
* ( ( -1 ) ( D22 ) )
*
* = ( conj(D21)*( D11 ) D21*( -1 ) )
* ( ( -1 ) ( D22 ) )
*
* where D11 = d22/d21,
* D22 = d11/conj(d21),
* D21 = T/d21,
* T = 1/(D22*D11-1).
*
* T/d21 is not formed, since it overflows when d21 is
* subnormal: each entry of the product is divided by
* d21 or conj(d21) and then scaled by T.
*
* (NOTE: No need to check for division by ZERO,
* since that was ensured earlier in pivot search:
* (a) d21 != 0, since in 2x2 pivot case(4)
Expand All @@ -850,16 +850,16 @@ SUBROUTINE CLAHEF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
D11 = W( K+1, K+1 ) / D21
D22 = W( K, K ) / CONJG( D21 )
T = ONE / ( REAL( D11*D22 )-ONE )
D21 = T / D21
*
* Update elements in columns A(k) and A(k+1) as
* dot products of rows of ( W(k) W(k+1) ) and columns
* of D**(-1)
*
DO 80 J = K + 2, N
A( J, K ) = CONJG( D21 )*
$ ( D11*W( J, K )-W( J, K+1 ) )
A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) /
$ CONJG( D21 ) )
A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) /
$ D21 )
80 CONTINUE
END IF
*
Expand Down
24 changes: 14 additions & 10 deletions SRC/clasyf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,8 +432,9 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* = 1/d21 * T * ( ( D11 ) ( -1 ) )
* ( ( -1 ) ( D22 ) )
*
* = D21 * ( ( D11 ) ( -1 ) )
* ( ( -1 ) ( D22 ) )
* T/d21 is not formed, since it overflows when d21 is
* subnormal: each entry of the product is divided by
* d21 and then scaled by T.
*
D21 = W( K-1, KW )
D11 = W( K, KW ) / D21
Expand All @@ -444,10 +445,11 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* dot products of rows of ( W(kw-1) W(kw) ) and columns
* of D**(-1)
*
D21 = T / D21
DO 20 J = 1, K - 2
A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
A( J, K-1 ) = T*( (D11*W( J, KW-1 )-W( J, KW ) ) /
$ D21 )
A( J, K ) = T*( ( D22*W( J, KW )-W( J, KW-1 ) ) /
$ D21 )
20 CONTINUE
END IF
*
Expand Down Expand Up @@ -712,22 +714,24 @@ SUBROUTINE CLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
* = 1/d21 * T * ( ( D11 ) ( -1 ) )
* ( ( -1 ) ( D22 ) )
*
* = D21 * ( ( D11 ) ( -1 ) )
* ( ( -1 ) ( D22 ) )
* T/d21 is not formed, since it overflows when d21 is
* subnormal: each entry of the product is divided by
* d21 and then scaled by T.
*
D21 = W( K+1, K )
D11 = W( K+1, K+1 ) / D21
D22 = W( K, K ) / D21
T = CONE / ( D11*D22-CONE )
D21 = T / D21
*
* Update elements in columns A(k) and A(k+1) as
* dot products of rows of ( W(k) W(k+1) ) and columns
* of D**(-1)
*
DO 80 J = K + 2, N
A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
A( J, K ) = T*( ( D11*W( J, K )-W( J, K+1 ) ) /
$ D21 )
A( J, K+1 ) = T*( ( D22*W( J, K+1 )-W( J, K ) ) /
$ D21 )
80 CONTINUE
END IF
*
Expand Down
18 changes: 8 additions & 10 deletions SRC/csptrf.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -375,13 +375,12 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO )
D22 = AP( K-1+( K-2 )*( K-1 ) / 2 ) / D12
D11 = AP( K+( K-1 )*K / 2 ) / D12
T = CONE / ( D11*D22-CONE )
D12 = T / D12
*
DO 50 J = K - 2, 1, -1
WKM1 = D12*( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
$ AP( J+( K-1 )*K / 2 ) )
WK = D12*( D22*AP( J+( K-1 )*K / 2 )-
$ AP( J+( K-2 )*( K-1 ) / 2 ) )
WKM1 = T*( ( D11*AP( J+( K-2 )*( K-1 ) / 2 )-
$ AP( J+( K-1 )*K / 2 ) ) / D12 )
WK = T*( ( D22*AP( J+( K-1 )*K / 2 )-
$ AP( J+( K-2 )*( K-1 ) / 2 ) ) / D12 )
DO 40 I = J, 1, -1
AP( I+( J-1 )*J / 2 ) = AP( I+( J-1 )*J / 2 ) -
$ AP( I+( K-1 )*K / 2 )*WK -
Expand Down Expand Up @@ -576,13 +575,12 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO )
D11 = AP( K+1+K*( 2*N-K-1 ) / 2 ) / D21
D22 = AP( K+( K-1 )*( 2*N-K ) / 2 ) / D21
T = CONE / ( D11*D22-CONE )
D21 = T / D21
*
DO 100 J = K + 2, N
WK = D21*( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
$ AP( J+K*( 2*N-K-1 ) / 2 ) )
WKP1 = D21*( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
$ AP( J+( K-1 )*( 2*N-K ) / 2 ) )
WK = T*( ( D11*AP( J+( K-1 )*( 2*N-K ) / 2 )-
$ AP( J+K*( 2*N-K-1 ) / 2 ) ) / D21 )
WKP1 = T*( ( D22*AP( J+K*( 2*N-K-1 ) / 2 )-
$ AP( J+( K-1 )*( 2*N-K ) / 2 ) ) / D21 )
DO 90 I = J, N
AP( I+( J-1 )*( 2*N-J ) / 2 ) = AP( I+( J-1 )*
$ ( 2*N-J ) / 2 ) - AP( I+( K-1 )*( 2*N-K ) /
Expand Down
Loading