phynexis-deps/lapack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

dbbcsd/sbbcsd/cbbcsd/zbbcsd: fix insufficiently accurate U2 in CS decomposition

Browse Source 
The Golub-Reinsch-SVD-style iteration in xBBCSD used an absolute
convergence threshold THRESH ~ 90*EPS to decide when off-diagonal
bulges are negligible and can be skipped.  For rows where the
diagonal entries are small (tiny singular values in the ratio),
this absolute threshold is too large relative to the diagonal,
causing the algorithm to stop chasing bulges prematurely and
producing inaccurate singular vectors (||X21 - U2 D2 V^*|| up to
50*EPS ||X21||).

Fix: change all bulge-convergence checks (RESTART flags in the
inner loop, initial bulge-chase decisions, and the IMAX-1 cleanup)
from absolute to relative by scaling THRESH by the adjacent diagonal
entries:

  RESTART11 = |B11E|^2+|BULGE|^2 <= (THRESH * MAX(|B11D(I-1)|,|B11D(I)|,UNFL))^2

This mirrors DBDSQR's relative convergence check |E| <= TOL*|D|.
Applied to all 4 bidiagonal blocks (B11/B21/B12/B22) at all 3
check-points in the iteration.

Fixes #965
This commit is contained in:
Julien Schueller
2026-06-11 19:15:11 +02:00
parent 077667a1be
commit a73f4a6524
4 changed files with 136 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SRC/cbbcsd.f
+34 -12
View File
@@ -616,7 +616,9 @@
*
* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
*
IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
IF( B11D(IMIN)**2+B11BULGE**2 .GT.
$ (THRESH*MAX( ABS(B11D(IMIN)),
$ ABS(B11D(IMIN+1)), UNFL ))**2 ) THEN
CALL SLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
$ RWORK(IU1CS+IMIN-1), R )
ELSE IF( MU .LE. NU ) THEN
@@ -626,7 +628,9 @@
CALL SLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
$ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
END IF
IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
IF( B21D(IMIN)**2+B21BULGE**2 .GT.
$ (THRESH*MAX( ABS(B21D(IMIN)),
$ ABS(B21D(IMIN+1)), UNFL ))**2 ) THEN
CALL SLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
$ RWORK(IU2CS+IMIN-1), R )
ELSE IF( NU .LT. MU ) THEN
@@ -690,10 +694,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11D(I-1)), ABS(B11D(I)),
$ UNFL ))**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21D(I-1)), ABS(B21D(I)),
$ UNFL ))**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(I-1)), ABS(B12D(I)),
$ UNFL ))**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(I-1)), ABS(B22D(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
@@ -775,10 +787,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11E(I)), ABS(B11D(I+1)),
$ UNFL ))**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12D(I)), ABS(B12E(I)),
$ UNFL ))**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21E(I)), ABS(B21D(I+1)),
$ UNFL ))**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22D(I)), ABS(B22E(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
@@ -866,8 +886,10 @@
*
* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
*
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(IMAX-1)), UNFL ))**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(IMAX-1)), UNFL ))**2
*
IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
CALL SLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),
SRC/dbbcsd.f
+34 -12
View File
@@ -616,7 +616,9 @@
*
* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
*
IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
IF( B11D(IMIN)**2+B11BULGE**2 .GT.
$ (THRESH*MAX( ABS(B11D(IMIN)),
$ ABS(B11D(IMIN+1)), UNFL ))**2 ) THEN
CALL DLARTGP( B11BULGE, B11D(IMIN), WORK(IU1SN+IMIN-1),
$ WORK(IU1CS+IMIN-1), R )
ELSE IF( MU .LE. NU ) THEN
@@ -626,7 +628,9 @@
CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
$ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
END IF
IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
IF( B21D(IMIN)**2+B21BULGE**2 .GT.
$ (THRESH*MAX( ABS(B21D(IMIN)),
$ ABS(B21D(IMIN+1)), UNFL ))**2 ) THEN
CALL DLARTGP( B21BULGE, B21D(IMIN), WORK(IU2SN+IMIN-1),
$ WORK(IU2CS+IMIN-1), R )
ELSE IF( NU .LT. MU ) THEN
@@ -690,10 +694,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11D(I-1)), ABS(B11D(I)),
$ UNFL ))**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21D(I-1)), ABS(B21D(I)),
$ UNFL ))**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(I-1)), ABS(B12D(I)),
$ UNFL ))**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(I-1)), ABS(B22D(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
@@ -776,10 +788,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11E(I)), ABS(B11D(I+1)),
$ UNFL ))**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12D(I)), ABS(B12E(I)),
$ UNFL ))**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21E(I)), ABS(B21D(I+1)),
$ UNFL ))**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22D(I)), ABS(B22E(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
@@ -863,8 +883,10 @@
*
* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
*
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(IMAX-1)), UNFL ))**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(IMAX-1)), UNFL ))**2
*
IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
CALL DLARTGP( Y2, Y1, WORK(IV2TSN+IMAX-1-1),
SRC/sbbcsd.f
+34 -12
View File
@@ -616,7 +616,9 @@
*
* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
*
IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
IF( B11D(IMIN)**2+B11BULGE**2 .GT.
$ (THRESH*MAX( ABS(B11D(IMIN)),
$ ABS(B11D(IMIN+1)), UNFL ))**2 ) THEN
CALL SLARTGP( B11BULGE, B11D(IMIN), WORK(IU1SN+IMIN-1),
$ WORK(IU1CS+IMIN-1), R )
ELSE IF( MU .LE. NU ) THEN
@@ -626,7 +628,9 @@
CALL SLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
$ WORK(IU1CS+IMIN-1), WORK(IU1SN+IMIN-1) )
END IF
IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
IF( B21D(IMIN)**2+B21BULGE**2 .GT.
$ (THRESH*MAX( ABS(B21D(IMIN)),
$ ABS(B21D(IMIN+1)), UNFL ))**2 ) THEN
CALL SLARTGP( B21BULGE, B21D(IMIN), WORK(IU2SN+IMIN-1),
$ WORK(IU2CS+IMIN-1), R )
ELSE IF( NU .LT. MU ) THEN
@@ -690,10 +694,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11D(I-1)), ABS(B11D(I)),
$ UNFL ))**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21D(I-1)), ABS(B21D(I)),
$ UNFL ))**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(I-1)), ABS(B12D(I)),
$ UNFL ))**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(I-1)), ABS(B22D(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
@@ -776,10 +788,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11E(I)), ABS(B11D(I+1)),
$ UNFL ))**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12D(I)), ABS(B12E(I)),
$ UNFL ))**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21E(I)), ABS(B21D(I+1)),
$ UNFL ))**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22D(I)), ABS(B22E(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
@@ -863,8 +883,10 @@
*
* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
*
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(IMAX-1)), UNFL ))**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(IMAX-1)), UNFL ))**2
*
IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
CALL SLARTGP( Y2, Y1, WORK(IV2TSN+IMAX-1-1),
SRC/zbbcsd.f
+34 -12
View File
@@ -615,7 +615,9 @@
*
* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
*
IF( B11D(IMIN)**2+B11BULGE**2 .GT. THRESH**2 ) THEN
IF( B11D(IMIN)**2+B11BULGE**2 .GT.
$ (THRESH*MAX( ABS(B11D(IMIN)),
$ ABS(B11D(IMIN+1)), UNFL ))**2 ) THEN
CALL DLARTGP( B11BULGE, B11D(IMIN), RWORK(IU1SN+IMIN-1),
$ RWORK(IU1CS+IMIN-1), R )
ELSE IF( MU .LE. NU ) THEN
@@ -625,7 +627,9 @@
CALL DLARTGS( B12D( IMIN ), B12E( IMIN ), NU,
$ RWORK(IU1CS+IMIN-1), RWORK(IU1SN+IMIN-1) )
END IF
IF( B21D(IMIN)**2+B21BULGE**2 .GT. THRESH**2 ) THEN
IF( B21D(IMIN)**2+B21BULGE**2 .GT.
$ (THRESH*MAX( ABS(B21D(IMIN)),
$ ABS(B21D(IMIN+1)), UNFL ))**2 ) THEN
CALL DLARTGP( B21BULGE, B21D(IMIN), RWORK(IU2SN+IMIN-1),
$ RWORK(IU2CS+IMIN-1), R )
ELSE IF( NU .LT. MU ) THEN
@@ -689,10 +693,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE. THRESH**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE. THRESH**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11E(I-1)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11D(I-1)), ABS(B11D(I)),
$ UNFL ))**2
RESTART21 = B21E(I-1)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21D(I-1)), ABS(B21D(I)),
$ UNFL ))**2
RESTART12 = B12D(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(I-1)), ABS(B12D(I)),
$ UNFL ))**2
RESTART22 = B22D(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(I-1)), ABS(B22D(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
@@ -774,10 +786,18 @@
* Determine if there are bulges to chase or if a new direct
* summand has been reached
*
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE. THRESH**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE. THRESH**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART11 = B11D(I)**2 + B11BULGE**2 .LE.
$ (THRESH*MAX( ABS(B11E(I)), ABS(B11D(I+1)),
$ UNFL ))**2
RESTART12 = B12E(I-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12D(I)), ABS(B12E(I)),
$ UNFL ))**2
RESTART21 = B21D(I)**2 + B21BULGE**2 .LE.
$ (THRESH*MAX( ABS(B21E(I)), ABS(B21D(I+1)),
$ UNFL ))**2
RESTART22 = B22E(I-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22D(I)), ABS(B22E(I)),
$ UNFL ))**2
*
* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
@@ -865,8 +885,10 @@
*
* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
*
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE. THRESH**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE. THRESH**2
RESTART12 = B12D(IMAX-1)**2 + B12BULGE**2 .LE.
$ (THRESH*MAX( ABS(B12E(IMAX-1)), UNFL ))**2
RESTART22 = B22D(IMAX-1)**2 + B22BULGE**2 .LE.
$ (THRESH*MAX( ABS(B22E(IMAX-1)), UNFL ))**2
*
IF( .NOT. RESTART12 .AND. .NOT. RESTART22 ) THEN
CALL DLARTGP( Y2, Y1, RWORK(IV2TSN+IMAX-1-1),