📄 sgesvd.f
字号:
*
* Copy R to WORK(IR), zeroing out below it
*
CALL SLACPY( 'U', N, N, A, LDA, WORK( IR ),
$ LDWRKR )
CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
$ WORK( IR+1 ), LDWRKR )
*
* Generate Q in A
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IE = ITAU
ITAUQ = IE + N
ITAUP = ITAUQ + N
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IR)
* (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
*
CALL SGEBRD( N, N, WORK( IR ), LDWRKR, S,
$ WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate left vectors bidiagonalizing R in WORK(IR)
* (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
*
CALL SORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
IWORK = IE + N
*
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IR)
* (Workspace: need N*N+BDSPAC)
*
CALL SBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
$ 1, WORK( IR ), LDWRKR, DUM, 1,
$ WORK( IWORK ), INFO )
*
* Multiply Q in A by left singular vectors of R in
* WORK(IR), storing result in U
* (Workspace: need N*N)
*
CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
$ WORK( IR ), LDWRKR, ZERO, U, LDU )
*
ELSE
*
* Insufficient workspace for a fast algorithm
*
ITAU = 1
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
* (Workspace: need 2*N, prefer N+N*NB)
*
CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
* (Workspace: need 2*N, prefer N+N*NB)
*
CALL SORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IE = ITAU
ITAUQ = IE + N
ITAUP = ITAUQ + N
IWORK = ITAUP + N
*
* Zero out below R in A
*
CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
$ LDA )
*
* Bidiagonalize R in A
* (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply Q in U by left vectors bidiagonalizing R
* (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
IWORK = IE + N
*
* Perform bidiagonal QR iteration, computing left
* singular vectors of A in U
* (Workspace: need BDSPAC)
*
CALL SBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM,
$ 1, U, LDU, DUM, 1, WORK( IWORK ),
$ INFO )
*
END IF
*
ELSE IF( WNTVO ) THEN
*
* Path 5 (M much larger than N, JOBU='S', JOBVT='O')
* N left singular vectors to be computed in U and
* N right singular vectors to be overwritten on A
*
IF( LWORK.GE.2*N*N+MAX( 4*N, BDSPAC ) ) THEN
*
* Sufficient workspace for a fast algorithm
*
IU = 1
IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is LDA by N
*
LDWRKU = LDA
IR = IU + LDWRKU*N
LDWRKR = LDA
ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
*
* WORK(IU) is LDA by N and WORK(IR) is N by N
*
LDWRKU = LDA
IR = IU + LDWRKU*N
LDWRKR = N
ELSE
*
* WORK(IU) is N by N and WORK(IR) is N by N
*
LDWRKU = N
IR = IU + LDWRKU*N
LDWRKR = N
END IF
ITAU = IR + LDWRKR*N
IWORK = ITAU + N
*
* Compute A=Q*R
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
*
CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Copy R to WORK(IU), zeroing out below it
*
CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
$ LDWRKU )
CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
$ WORK( IU+1 ), LDWRKU )
*
* Generate Q in A
* (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
*
CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IE = ITAU
ITAUQ = IE + N
ITAUP = ITAUQ + N
IWORK = ITAUP + N
*
* Bidiagonalize R in WORK(IU), copying result to
* WORK(IR)
* (Workspace: need 2*N*N+4*N,
* prefer 2*N*N+3*N+2*N*NB)
*
CALL SGEBRD( N, N, WORK( IU ), LDWRKU, S,
$ WORK( IE ), WORK( ITAUQ ),
$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
CALL SLACPY( 'U', N, N, WORK( IU ), LDWRKU,
$ WORK( IR ), LDWRKR )
*
* Generate left bidiagonalizing vectors in WORK(IU)
* (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
*
CALL SORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
$ WORK( ITAUQ ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right bidiagonalizing vectors in WORK(IR)
* (Workspace: need 2*N*N+4*N-1,
* prefer 2*N*N+3*N+(N-1)*NB)
*
CALL SORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
$ WORK( ITAUP ), WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
IWORK = IE + N
*
* Perform bidiagonal QR iteration, computing left
* singular vectors of R in WORK(IU) and computing
* right singular vectors of R in WORK(IR)
* (Workspace: need 2*N*N+BDSPAC)
*
CALL SBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
$ WORK( IR ), LDWRKR, WORK( IU ),
$ LDWRKU, DUM, 1, WORK( IWORK ), INFO )
*
* Multiply Q in A by left singular vectors of R in
* WORK(IU), storing result in U
* (Workspace: need N*N)
*
CALL SGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
$ WORK( IU ), LDWRKU, ZERO, U, LDU )
*
* Copy right singular vectors of R to A
* (Workspace: need N*N)
*
CALL SLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
$ LDA )
*
ELSE
*
* Insufficient workspace for a fast algorithm
*
ITAU = 1
IWORK = ITAU + N
*
* Compute A=Q*R, copying result to U
* (Workspace: need 2*N, prefer N+N*NB)
*
CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
CALL SLACPY( 'L', M, N, A, LDA, U, LDU )
*
* Generate Q in U
* (Workspace: need 2*N, prefer N+N*NB)
*
CALL SORGQR( M, N, N, U, LDU, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IE = ITAU
ITAUQ = IE + N
ITAUP = ITAUQ + N
IWORK = ITAUP + N
*
* Zero out below R in A
*
CALL SLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
$ LDA )
*
* Bidiagonalize R in A
* (Workspace: need 4*N, prefer 3*N+2*N*NB)
*
CALL SGEBRD( N, N, A, LDA, S, WORK( IE ),
$ WORK( ITAUQ ), WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Multiply Q in U by left vectors bidiagonalizing R
* (Workspace: need 3*N+M, prefer 3*N+M*NB)
*
CALL SORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
$ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
$ LWORK-IWORK+1, IERR )
*
* Generate right vectors bidiagonalizing R in A
* (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
*
CALL SORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IWORK = IE + N
*
* Perform bidiagonal QR iteration, computing left
* singular vectors of A in U and computing right
* singular vectors of A in A
* (Workspace: need BDSPAC)
*
CALL SBDSQR( 'U', N, N, M, 0, S, WORK( IE ), A,
$ LDA, U, LDU, DUM, 1, WORK( IWORK ),
$ INFO )
*
END IF
*
ELSE IF( WNTVAS ) THEN
*
* Path 6 (M much larger than N, JOBU='S', JOBVT='S'
* or 'A')
* N left singular vectors to be computed in U and
* N right singular vectors to be computed in VT
*
IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
*
* Sufficient workspace for a fast algorithm
*
IU = 1
IF( LWORK.GE.WRKBL+LDA*N ) THEN
*
* WORK(IU) is LDA by N
*
LDWRKU = LDA
ELSE
*
* WORK(IU) is N by N
*
LDWRKU = N
END IF
ITAU = IU + LDWRKU*N
IWORK = ITAU + N
*
* Compute A=Q*R
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
CALL SGEQRF( M, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
*
* Copy R to WORK(IU), zeroing out below it
*
CALL SLACPY( 'U', N, N, A, LDA, WORK( IU ),
$ LDWRKU )
CALL SLASET( 'L', N-1, N-1, ZERO, ZERO,
$ WORK( IU+1 ), LDWRKU )
*
* Generate Q in A
* (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
*
CALL SORGQR( M, N, N, A, LDA, WORK( ITAU ),
$ WORK( IWORK ), LWORK-IWORK+1, IERR )
IE = ITAU
ITAUQ = IE + N
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -