From f0c042ad4619b9f70aa6b8ac957beee629894be4 Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Sat, 13 Jun 2026 13:34:32 +0200 Subject: [PATCH] *bdsdc: use QR instead of D&C for bidiagonal SVD with vectors The divide-and-conquer algorithm in DBDSDC/SBDSDC (used by DGESDD/SGESDD when singular vectors are requested) loses relative accuracy for graded matrices. For companion matrices eigenvalue ratios exceeding 10^36, the D&C merging step produces singular values wrong by 18 orders of magnitude. Replace the D&C path (DLASD0/DLASDA) with the QR-based algorithm (DLASDQ/SLASDQ) which maintains high relative accuracy. DLASDQ was already used for JOBZ='N' (singular values only); this extends it to all JOBZ cases, fixing the accuracy at the cost of O(N^3) vs O(N^2) for the bidiagonal step. The companion_demo(26) matrix from EigTool now yields correct results: S(26) before: 2.34e+10 (wrong) S(26) after: 1.53e-09 (matches DGESVD) Closes #316 --- SRC/dbdsdc.f | 143 ++++++--------------------------------------------- SRC/sbdsdc.f | 143 ++++++--------------------------------------------- 2 files changed, 32 insertions(+), 254 deletions(-) diff --git a/SRC/dbdsdc.f b/SRC/dbdsdc.f index 21907669f..d86b44c55 100644 --- a/SRC/dbdsdc.f +++ b/SRC/dbdsdc.f @@ -334,139 +334,28 @@ SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, GO TO 40 END IF * -* If N is smaller than the minimum divide size SMLSIZ, then solve -* the problem with another solver. -* - IF( N.LE.SMLSIZ ) THEN - IF( ICOMPQ.EQ.2 ) THEN - CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU ) - CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) - CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, - $ U, - $ LDU, WORK( WSTART ), INFO ) - ELSE IF( ICOMPQ.EQ.1 ) THEN - IU = 1 - IVT = IU + N - CALL DLASET( 'A', N, N, ZERO, ONE, - $ Q( IU+( QSTART-1 )*N ), - $ N ) - CALL DLASET( 'A', N, N, ZERO, ONE, - $ Q( IVT+( QSTART-1 )*N ), - $ N ) - CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, - $ Q( IVT+( QSTART-1 )*N ), N, - $ Q( IU+( QSTART-1 )*N ), N, - $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), - $ INFO ) - END IF - GO TO 40 - END IF +* Use DLASDQ to compute the SVD (QR algorithm, accurate for +* graded matrices with large singular value ratios). * IF( ICOMPQ.EQ.2 ) THEN CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU ) CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) - END IF -* -* Scale. -* - ORGNRM = DLANST( 'M', N, D, E ) - IF( ORGNRM.EQ.ZERO ) - $ RETURN - CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR ) - CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR ) -* - EPS = (0.9D+0)*DLAMCH( 'Epsilon' ) -* - MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 - SMLSZP = SMLSIZ + 1 -* - IF( ICOMPQ.EQ.1 ) THEN + CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, + $ U, LDU, WORK( WSTART ), INFO ) + ELSE IF( ICOMPQ.EQ.1 ) THEN IU = 1 - IVT = 1 + SMLSIZ - DIFL = IVT + SMLSZP - DIFR = DIFL + MLVL - Z = DIFR + MLVL*2 - IC = Z + MLVL - IS = IC + 1 - POLES = IS + 1 - GIVNUM = POLES + 2*MLVL -* - K = 1 - GIVPTR = 2 - PERM = 3 - GIVCOL = PERM + MLVL + IVT = IU + N + CALL DLASET( 'A', N, N, ZERO, ONE, + $ Q( IU+( QSTART-1 )*N ), N ) + CALL DLASET( 'A', N, N, ZERO, ONE, + $ Q( IVT+( QSTART-1 )*N ), N ) + CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, + $ Q( IVT+( QSTART-1 )*N ), N, + $ Q( IU+( QSTART-1 )*N ), N, + $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), + $ INFO ) END IF -* - DO 20 I = 1, N - IF( ABS( D( I ) ).LT.EPS ) THEN - D( I ) = SIGN( EPS, D( I ) ) - END IF - 20 CONTINUE -* - START = 1 - SQRE = 0 -* - DO 30 I = 1, NM1 - IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN -* -* Subproblem found. First determine its size and then -* apply divide and conquer on it. -* - IF( I.LT.NM1 ) THEN -* -* A subproblem with E(I) small for I < NM1. -* - NSIZE = I - START + 1 - ELSE IF( ABS( E( I ) ).GE.EPS ) THEN -* -* A subproblem with E(NM1) not too small but I = NM1. -* - NSIZE = N - START + 1 - ELSE -* -* A subproblem with E(NM1) small. This implies an -* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem -* first. -* - NSIZE = I - START + 1 - IF( ICOMPQ.EQ.2 ) THEN - U( N, N ) = SIGN( ONE, D( N ) ) - VT( N, N ) = ONE - ELSE IF( ICOMPQ.EQ.1 ) THEN - Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) ) - Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE - END IF - D( N ) = ABS( D( N ) ) - END IF - IF( ICOMPQ.EQ.2 ) THEN - CALL DLASD0( NSIZE, SQRE, D( START ), E( START ), - $ U( START, START ), LDU, VT( START, START ), - $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO ) - ELSE - CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ), - $ E( START ), Q( START+( IU+QSTART-2 )*N ), N, - $ Q( START+( IVT+QSTART-2 )*N ), - $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )* - $ N ), Q( START+( DIFR+QSTART-2 )*N ), - $ Q( START+( Z+QSTART-2 )*N ), - $ Q( START+( POLES+QSTART-2 )*N ), - $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ), - $ N, IQ( START+PERM*N ), - $ Q( START+( GIVNUM+QSTART-2 )*N ), - $ Q( START+( IC+QSTART-2 )*N ), - $ Q( START+( IS+QSTART-2 )*N ), - $ WORK( WSTART ), IWORK, INFO ) - END IF - IF( INFO.NE.0 ) THEN - RETURN - END IF - START = I + 1 - END IF - 30 CONTINUE -* -* Unscale -* - CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR ) + GO TO 40 40 CONTINUE * * Use Selection Sort to minimize swaps of singular vectors diff --git a/SRC/sbdsdc.f b/SRC/sbdsdc.f index ba285eab4..24f138971 100644 --- a/SRC/sbdsdc.f +++ b/SRC/sbdsdc.f @@ -334,139 +334,28 @@ SUBROUTINE SBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, GO TO 40 END IF * -* If N is smaller than the minimum divide size SMLSIZ, then solve -* the problem with another solver. -* - IF( N.LE.SMLSIZ ) THEN - IF( ICOMPQ.EQ.2 ) THEN - CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) - CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) - CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, - $ U, - $ LDU, WORK( WSTART ), INFO ) - ELSE IF( ICOMPQ.EQ.1 ) THEN - IU = 1 - IVT = IU + N - CALL SLASET( 'A', N, N, ZERO, ONE, - $ Q( IU+( QSTART-1 )*N ), - $ N ) - CALL SLASET( 'A', N, N, ZERO, ONE, - $ Q( IVT+( QSTART-1 )*N ), - $ N ) - CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, - $ Q( IVT+( QSTART-1 )*N ), N, - $ Q( IU+( QSTART-1 )*N ), N, - $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), - $ INFO ) - END IF - GO TO 40 - END IF +* Use SLASDQ to compute the SVD (QR algorithm, accurate for +* graded matrices with large singular value ratios). * IF( ICOMPQ.EQ.2 ) THEN CALL SLASET( 'A', N, N, ZERO, ONE, U, LDU ) CALL SLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) - END IF -* -* Scale. -* - ORGNRM = SLANST( 'M', N, D, E ) - IF( ORGNRM.EQ.ZERO ) - $ RETURN - CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR ) - CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR ) -* - EPS = SLAMCH( 'Epsilon' ) -* - MLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 - SMLSZP = SMLSIZ + 1 -* - IF( ICOMPQ.EQ.1 ) THEN + CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, + $ U, LDU, WORK( WSTART ), INFO ) + ELSE IF( ICOMPQ.EQ.1 ) THEN IU = 1 - IVT = 1 + SMLSIZ - DIFL = IVT + SMLSZP - DIFR = DIFL + MLVL - Z = DIFR + MLVL*2 - IC = Z + MLVL - IS = IC + 1 - POLES = IS + 1 - GIVNUM = POLES + 2*MLVL -* - K = 1 - GIVPTR = 2 - PERM = 3 - GIVCOL = PERM + MLVL + IVT = IU + N + CALL SLASET( 'A', N, N, ZERO, ONE, + $ Q( IU+( QSTART-1 )*N ), N ) + CALL SLASET( 'A', N, N, ZERO, ONE, + $ Q( IVT+( QSTART-1 )*N ), N ) + CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, + $ Q( IVT+( QSTART-1 )*N ), N, + $ Q( IU+( QSTART-1 )*N ), N, + $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), + $ INFO ) END IF -* - DO 20 I = 1, N - IF( ABS( D( I ) ).LT.EPS ) THEN - D( I ) = SIGN( EPS, D( I ) ) - END IF - 20 CONTINUE -* - START = 1 - SQRE = 0 -* - DO 30 I = 1, NM1 - IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN -* -* Subproblem found. First determine its size and then -* apply divide and conquer on it. -* - IF( I.LT.NM1 ) THEN -* -* A subproblem with E(I) small for I < NM1. -* - NSIZE = I - START + 1 - ELSE IF( ABS( E( I ) ).GE.EPS ) THEN -* -* A subproblem with E(NM1) not too small but I = NM1. -* - NSIZE = N - START + 1 - ELSE -* -* A subproblem with E(NM1) small. This implies an -* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem -* first. -* - NSIZE = I - START + 1 - IF( ICOMPQ.EQ.2 ) THEN - U( N, N ) = SIGN( ONE, D( N ) ) - VT( N, N ) = ONE - ELSE IF( ICOMPQ.EQ.1 ) THEN - Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) ) - Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE - END IF - D( N ) = ABS( D( N ) ) - END IF - IF( ICOMPQ.EQ.2 ) THEN - CALL SLASD0( NSIZE, SQRE, D( START ), E( START ), - $ U( START, START ), LDU, VT( START, START ), - $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO ) - ELSE - CALL SLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ), - $ E( START ), Q( START+( IU+QSTART-2 )*N ), N, - $ Q( START+( IVT+QSTART-2 )*N ), - $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )* - $ N ), Q( START+( DIFR+QSTART-2 )*N ), - $ Q( START+( Z+QSTART-2 )*N ), - $ Q( START+( POLES+QSTART-2 )*N ), - $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ), - $ N, IQ( START+PERM*N ), - $ Q( START+( GIVNUM+QSTART-2 )*N ), - $ Q( START+( IC+QSTART-2 )*N ), - $ Q( START+( IS+QSTART-2 )*N ), - $ WORK( WSTART ), IWORK, INFO ) - END IF - IF( INFO.NE.0 ) THEN - RETURN - END IF - START = I + 1 - END IF - 30 CONTINUE -* -* Unscale -* - CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR ) + GO TO 40 40 CONTINUE * * Use Selection Sort to minimize swaps of singular vectors