Skip to content

*tgex2: loosen swap acceptance threshold to fix false rejections #1299

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
jschueller wants to merge 2 commits into
base: master
Choose a base branch
from
Open

*tgex2: loosen swap acceptance threshold to fix false rejections #1299

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
52f50ab
*tgex2: check raw Sylvester off-diagonal of both A and B in weak test
jschueller Jun 13, 2026
3aa1272
tgex2: use LANGE
jschueller Jul 23, 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
47 changes: 28 additions & 19 deletions SRC/dtgex2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,11 +246,11 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
PARAMETER ( WANDS = .TRUE. )
* ..
* .. Local Scalars ..
LOGICAL STRONG, WEAK
LOGICAL RAWACC, STRONG, WEAK
INTEGER I, IDUM, LINFO, M
DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORMA, DNORMB, DSCALE,
$ DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM,
$ THRESHA, THRESHB
DOUBLE PRECISION BQRA21, BQRB21, BRQA21, DDUM, DNORMA, DNORMB,
$ DSCALE, DSUM, EPS, F, G, SA, SB, SCALE,
$ SMLNUM, THRESHA, THRESHB
* ..
* .. Local Arrays ..
INTEGER IWORK( LDST + 2 )
Expand All @@ -261,8 +261,8 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
$ TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DGEMM, DGEQR2, DGERQ2, DLACPY, DLAGV2,
Expand Down Expand Up @@ -322,6 +322,8 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
* on 04/01/10.
* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
* Jim Demmel and Guillaume Revy. See forum post 1783.
* The weak test now checks both A and B off-diagonal blocks from the
* raw Sylvester result.
*
THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM )
THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM )
Expand Down Expand Up @@ -504,6 +506,16 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
CALL DGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO,
$ T,
$ LDST )
*
* Compute F-norm(S21) and F-norm(T21) for the raw Sylvester
* result before the QR/RQ cleanup.
*
BQRA21 = DLANGE( 'F', N1, N2, S( N2+1, 1 ), LDST, WORK )
*
BQRB21 = DLANGE( 'F', N1, N2, T( N2+1, 1 ), LDST, WORK )
*
RAWACC = BQRA21.LE.THRESHA .AND. BQRB21.LE.THRESHB
*
CALL DLACPY( 'F', M, M, S, LDST, SCPY, LDST )
CALL DLACPY( 'F', M, M, T, LDST, TCPY, LDST )
CALL DLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
Expand All @@ -528,12 +540,7 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
*
* Compute F-norm(S21) in BRQA21. (T21 is 0.)
*
DSCALE = ZERO
DSUM = ONE
DO 30 I = 1, N2
CALL DLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
30 CONTINUE
BRQA21 = DSCALE*SQRT( DSUM )
BRQA21 = DLANGE( 'F', N1, N2, S( N2+1, 1 ), LDST, WORK )
*
* Triangularize the B-part by a QR factorization.
* Apply transformation (from right) to A-part, giving S.
Expand All @@ -552,23 +559,25 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
*
* Compute F-norm(S21) in BQRA21. (T21 is 0.)
*
DSCALE = ZERO
DSUM = ONE
DO 40 I = 1, N2
CALL DLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
40 CONTINUE
BQRA21 = DSCALE*SQRT( DSUM )
BQRA21 = DLANGE( 'F', N1, N2, SCPY( N2+1, 1 ), LDST, WORK )
*
* Decide which method to use.
* Weak stability test:
* F-norm(S21) <= O(EPS * F-norm((S)))
* and
* F-norm(T21) <= O(EPS * F-norm((T)))
*
* First check the raw Sylvester result (before RQ/QR cleanup).
* If both A and B off-diagonal blocks are below threshold the
* swap is accepted unconditionally -- the RQ/QR cleanup that
* follows cannot trigger a false rejection.
*
IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESHA ) THEN
CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST )
CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST )
CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
ELSE IF( BRQA21.GE.THRESHA ) THEN
ELSE IF( BRQA21.GE.THRESHA .AND. .NOT.RAWACC ) THEN
GO TO 70
END IF
*
Expand Down
48 changes: 28 additions & 20 deletions SRC/stgex2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,12 +246,11 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
PARAMETER ( WANDS = .TRUE. )
* ..
* .. Local Scalars ..
LOGICAL STRONG, WEAK
LOGICAL RAWACC, STRONG, WEAK
INTEGER I, IDUM, LINFO, M
REAL BQRA21, BRQA21, DDUM, DNORMA, DNORMB,
$ DSCALE,
$ DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM,
$ THRESHA, THRESHB
REAL BQRA21, BQRB21, BRQA21, DDUM, DNORMA, DNORMB,
$ DSCALE, DSUM, EPS, F, G, SA, SB, SCALE,
$ SMLNUM, THRESHA, THRESHB
* ..
* .. Local Arrays ..
INTEGER IWORK( LDST + 2 )
Expand All @@ -262,8 +261,8 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
$ TAUL( LDST ), TAUR( LDST ), TCPY( LDST, LDST )
* ..
* .. External Functions ..
REAL SLAMCH
EXTERNAL SLAMCH
REAL SLAMCH, SLANGE
EXTERNAL SLAMCH, SLANGE
* ..
* .. External Subroutines ..
EXTERNAL SGEMM, SGEQR2, SGERQ2, SLACPY, SLAGV2,
Expand Down Expand Up @@ -323,6 +322,8 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
* on 04/01/10.
* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
* Jim Demmel and Guillaume Revy. See forum post 1783.
* The weak test now checks both A and B off-diagonal blocks from the
* raw Sylvester result.
*
THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM )
THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM )
Expand Down Expand Up @@ -505,6 +506,16 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
CALL SGEMM( 'N', 'T', M, M, M, ONE, WORK, M, IR, LDST, ZERO,
$ T,
$ LDST )
*
* Compute F-norm(S21) and F-norm(T21) for the raw Sylvester
* result before the QR/RQ cleanup.
*
BQRA21 = SLANGE( 'F', N1, N2, S( N2+1, 1 ), LDST, WORK )
*
BQRB21 = SLANGE( 'F', N1, N2, T( N2+1, 1 ), LDST, WORK )
*
RAWACC = BQRA21.LE.THRESHA .AND. BQRB21.LE.THRESHB
*
CALL SLACPY( 'F', M, M, S, LDST, SCPY, LDST )
CALL SLACPY( 'F', M, M, T, LDST, TCPY, LDST )
CALL SLACPY( 'F', M, M, IR, LDST, IRCOP, LDST )
Expand All @@ -529,12 +540,7 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
*
* Compute F-norm(S21) in BRQA21. (T21 is 0.)
*
DSCALE = ZERO
DSUM = ONE
DO 30 I = 1, N2
CALL SLASSQ( N1, S( N2+1, I ), 1, DSCALE, DSUM )
30 CONTINUE
BRQA21 = DSCALE*SQRT( DSUM )
BRQA21 = SLANGE( 'F', N1, N2, S( N2+1, 1 ), LDST, WORK )
*
* Triangularize the B-part by a QR factorization.
* Apply transformation (from right) to A-part, giving S.
Expand All @@ -553,23 +559,25 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
*
* Compute F-norm(S21) in BQRA21. (T21 is 0.)
*
DSCALE = ZERO
DSUM = ONE
DO 40 I = 1, N2
CALL SLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM )
40 CONTINUE
BQRA21 = DSCALE*SQRT( DSUM )
BQRA21 = SLANGE( 'F', N1, N2, SCPY( N2+1, 1 ), LDST, WORK )
*
* Decide which method to use.
* Weak stability test:
* F-norm(S21) <= O(EPS * F-norm((S)))
* and
* F-norm(T21) <= O(EPS * F-norm((T)))
*
* First check the raw Sylvester result (before RQ/QR cleanup).
* If both A and B off-diagonal blocks are below threshold the
* swap is accepted unconditionally -- the RQ/QR cleanup that
* follows cannot trigger a false rejection.
*
IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESHA ) THEN
CALL SLACPY( 'F', M, M, SCPY, LDST, S, LDST )
CALL SLACPY( 'F', M, M, TCPY, LDST, T, LDST )
CALL SLACPY( 'F', M, M, IRCOP, LDST, IR, LDST )
CALL SLACPY( 'F', M, M, LICOP, LDST, LI, LDST )
ELSE IF( BRQA21.GE.THRESHA ) THEN
ELSE IF( BRQA21.GE.THRESHA .AND. .NOT.RAWACC ) THEN
GO TO 70
END IF
*
Expand Down
Loading