From 0c217b68710970643c6f0fc648564d30e917a164 Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Sat, 13 Jun 2026 13:09:33 +0200 Subject: [PATCH] *tgex2: loosen swap acceptance threshold to fix false rejections The weak stability test in *TGEX2 rejects a block swap when the off-diagonal block norm BQRA21 exceeds THRESHA = TWENTY*EPS*SA. For certain matrices (reported by MathWorks from MATLAB's ORDQZ), this threshold is too tight -- BQRA21 can be ~5x larger than TWENTY*EPS*SA due to round-off in the preceding QZ steps, even though the swap is mathematically valid (verified via Sylvester equation-based proof). Increase the safety factor from TWENTY (20) to HUNDRED (100), which fixes all known failing cases (3x3 and 4x4 pencils) while remaining well below the BRQA21 bound that guards against truly invalid swaps. Closes #1201 --- SRC/ctgex2.f | 10 ++++++---- SRC/dtgex2.f | 10 ++++++---- SRC/stgex2.f | 10 ++++++---- SRC/ztgex2.f | 10 ++++++---- 4 files changed, 24 insertions(+), 16 deletions(-) diff --git a/SRC/ctgex2.f b/SRC/ctgex2.f index 99b18109b..c289dd674 100644 --- a/SRC/ctgex2.f +++ b/SRC/ctgex2.f @@ -206,8 +206,8 @@ SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, COMPLEX CZERO, CONE PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), $ CONE = ( 1.0E+0, 0.0E+0 ) ) - REAL TWENTY - PARAMETER ( TWENTY = 2.0E+1 ) + REAL HUNDRED + PARAMETER ( HUNDRED = 1.0E+2 ) INTEGER LDST PARAMETER ( LDST = 2 ) LOGICAL WANDS @@ -273,9 +273,11 @@ SUBROUTINE CTGEX2( 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. +* Then from TWENTY to HUNDRED on 03/20/26. +* "Bug" reported by Heiko Weichelt and Bobby Cheng from MathWorks. * - THRESHA = MAX( TWENTY*EPS*SA, SMLNUM ) - THRESHB = MAX( TWENTY*EPS*SB, SMLNUM ) + THRESHA = MAX( HUNDRED*EPS*SA, SMLNUM ) + THRESHB = MAX( HUNDRED*EPS*SB, SMLNUM ) * * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks * using Givens rotations and perform the swap tentatively. diff --git a/SRC/dtgex2.f b/SRC/dtgex2.f index 8c01a4f5e..1fbbb118f 100644 --- a/SRC/dtgex2.f +++ b/SRC/dtgex2.f @@ -238,8 +238,8 @@ SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) - DOUBLE PRECISION TWENTY - PARAMETER ( TWENTY = 2.0D+01 ) + DOUBLE PRECISION HUNDRED + PARAMETER ( HUNDRED = 1.0D+02 ) INTEGER LDST PARAMETER ( LDST = 4 ) LOGICAL WANDS @@ -322,9 +322,11 @@ 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. +* Then from TWENTY to HUNDRED on 03/20/26. +* "Bug" reported by Heiko Weichelt and Bobby Cheng from MathWorks. * - THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM ) - THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM ) + THRESHA = MAX( HUNDRED*EPS*DNORMA, SMLNUM ) + THRESHB = MAX( HUNDRED*EPS*DNORMB, SMLNUM ) * IF( M.EQ.2 ) THEN * diff --git a/SRC/stgex2.f b/SRC/stgex2.f index d6a11841d..ca5e1c577 100644 --- a/SRC/stgex2.f +++ b/SRC/stgex2.f @@ -238,8 +238,8 @@ SUBROUTINE STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) - REAL TWENTY - PARAMETER ( TWENTY = 2.0E+01 ) + REAL HUNDRED + PARAMETER ( HUNDRED = 1.0E+02 ) INTEGER LDST PARAMETER ( LDST = 4 ) LOGICAL WANDS @@ -323,9 +323,11 @@ 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. +* Then from TWENTY to HUNDRED on 03/20/26. +* "Bug" reported by Heiko Weichelt and Bobby Cheng from MathWorks. * - THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM ) - THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM ) + THRESHA = MAX( HUNDRED*EPS*DNORMA, SMLNUM ) + THRESHB = MAX( HUNDRED*EPS*DNORMB, SMLNUM ) * IF( M.EQ.2 ) THEN * diff --git a/SRC/ztgex2.f b/SRC/ztgex2.f index 5f97848d4..6fa209791 100644 --- a/SRC/ztgex2.f +++ b/SRC/ztgex2.f @@ -206,8 +206,8 @@ SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), $ CONE = ( 1.0D+0, 0.0D+0 ) ) - DOUBLE PRECISION TWENTY - PARAMETER ( TWENTY = 2.0D+1 ) + DOUBLE PRECISION HUNDRED + PARAMETER ( HUNDRED = 1.0D+2 ) INTEGER LDST PARAMETER ( LDST = 2 ) LOGICAL WANDS @@ -273,9 +273,11 @@ SUBROUTINE ZTGEX2( 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. +* Then from TWENTY to HUNDRED on 03/20/26. +* "Bug" reported by Heiko Weichelt and Bobby Cheng from MathWorks. * - THRESHA = MAX( TWENTY*EPS*SA, SMLNUM ) - THRESHB = MAX( TWENTY*EPS*SB, SMLNUM ) + THRESHA = MAX( HUNDRED*EPS*SA, SMLNUM ) + THRESHB = MAX( HUNDRED*EPS*SB, SMLNUM ) * * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks * using Givens rotations and perform the swap tentatively.