Skip to content

*bdsdc: use QR instead of D&C for bidiagonal SVD with vectors #1300

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 1 commit into
base: master
Choose a base branch
from
+32 −254
Open

*bdsdc: use QR instead of D&C for bidiagonal SVD with vectors #1300

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
143 changes: 16 additions & 127 deletions SRC/dbdsdc.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
143 changes: 16 additions & 127 deletions SRC/sbdsdc.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading