zgesdd.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

*
               ITAU = 1
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (CWorkspace: need 2*M, prefer M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGELQF', M, N, 0, 0, NB )
               CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Zero out above L
*
               CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
     $                      LDA )
               IE = 1
               ITAUQ = 1
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in A
*              (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
*              (RWorkspace: need M)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', M, M, 0, 0, NB )
               CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
               NRWORK = IE + M
*
*              Perform bidiagonal SVD, compute singular values only
*              (CWorkspace: 0)
*              (RWorkspace: need BDSPAC)
*
               CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
*
            ELSE IF( WNTQO ) THEN
*
*              Path 2t (N much larger than M, JOBZ='O')
*              M right singular vectors to be overwritten on A and
*              M left singular vectors to be computed in U
*
               IVT = 1
               LDWKVT = M
*
*              WORK(IVT) is M by M
*
               IL = IVT + LDWKVT*M
               IF( LWORK.GE.M*N+M*M+3*M ) THEN
*
*                 WORK(IL) M by N
*
                  LDWRKL = M
                  CHUNK = N
               ELSE
*
*                 WORK(IL) is M by CHUNK
*
                  LDWRKL = M
                  CHUNK = ( LWORK-M*M-3*M ) / M
               END IF
               ITAU = IL + LDWRKL*CHUNK
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (CWorkspace: need 2*M, prefer M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGELQF', M, N, 0, 0, NB )
               CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy L to WORK(IL), zeroing about above it
*
               CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
               CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
     $                      WORK( IL+LDWRKL ), LDWRKL )
*
*              Generate Q in A
*              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGLQ', ' ', M, N, M, -1 )
               OPS = OPS + DOPLA( 'ZUNGLQ', M, N, M, 0, NB )
               CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = 1
               ITAUQ = ITAU
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in WORK(IL)
*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
*              (RWorkspace: need M)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', M, M, 0, 0, NB )
               CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in RWORK(IRU) and computing right
*              singular vectors of bidiagonal matrix in RWORK(IRVT)
*              (CWorkspace: need 0)
*              (RWorkspace: need BDSPAC)
*
               IRU = IE + M
               IRVT = IRU + M*M
               NRWORK = IRVT + M*M
               CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
     $                      M, RWORK( IRVT ), M, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
*
*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
*              Overwrite WORK(IU) by the left singular vectors of L
*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', M, M, M, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
*              Overwrite WORK(IVT) by the right singular vectors of L
*              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
     $                      LDWKVT )
               NB = ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRC', M, M, M, 0, NB )
               CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Multiply right singular vectors of L in WORK(IL) by Q
*              in A, storing result in WORK(IL) and copying to A
*              (CWorkspace: need 2*M*M, prefer M*M+M*N))
*              (RWorkspace: 0)
*
               DO 40 I = 1, N, CHUNK
                  BLK = MIN( N-I+1, CHUNK )
                  OPS = OPS + DOPBL3( 'ZGEMM ', M, BLK, M )
                  CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
     $                        A( 1, I ), LDA, CZERO, WORK( IL ),
     $                        LDWRKL )
                  CALL ZLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
     $                         A( 1, I ), LDA )
   40          CONTINUE
*
            ELSE IF( WNTQS ) THEN
*
*             Path 3t (N much larger than M, JOBZ='S')
*             M right singular vectors to be computed in VT and
*             M left singular vectors to be computed in U
*
               IL = 1
*
*              WORK(IL) is M by M
*
               LDWRKL = M
               ITAU = IL + LDWRKL*M
               NWORK = ITAU + M
*
*              Compute A=L*Q
*              (CWorkspace: need 2*M, prefer M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGELQF', M, N, 0, 0, NB )
               CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy L to WORK(IL), zeroing out above it
*
               CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
               CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
     $                      WORK( IL+LDWRKL ), LDWRKL )
*
*              Generate Q in A
*              (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGLQ', ' ', M, N, M, -1 )
               OPS = OPS + DOPLA( 'ZUNGLQ', M, N, M, 0, NB )
               CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = 1
               ITAUQ = ITAU
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in WORK(IL)
*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
*              (RWorkspace: need M)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', M, M, 0, 0, NB )
               CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in RWORK(IRU) and computing right
*              singular vectors of bidiagonal matrix in RWORK(IRVT)
*              (CWorkspace: need 0)
*              (RWorkspace: need BDSPAC)
*
               IRU = IE + M
               IRVT = IRU + M*M
               NRWORK = IRVT + M*M
               CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
     $                      M, RWORK( IRVT ), M, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
*
*              Copy real matrix RWORK(IRU) to complex matrix U
*              Overwrite U by left singular vectors of L
*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', M, M, M, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy real matrix RWORK(IRVT) to complex matrix VT
*              Overwrite VT by left singular vectors of L
*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
               NB = ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRC', M, M, M, 0, NB )
               CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy VT to WORK(IL), multiply right singular vectors of L
*              in WORK(IL) by Q in A, storing result in VT
*              (CWorkspace: need M*M)
*              (RWorkspace: 0)
*
               CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
               OPS = OPS + DOPBL3( 'ZGEMM ', M, N, M )
               CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
     $                     A, LDA, CZERO, VT, LDVT )
*
            ELSE IF( WNTQA ) THEN
*
*              Path 9t (N much larger than M, JOBZ='A')
*              N right singular vectors to be computed in VT and
*              M left singular vectors to be computed in U
*
               IVT = 1
*
*              WORK(IVT) is M by M
*
               LDWKVT = M
               ITAU = IVT + LDWKVT*M
               NWORK = ITAU + M
*
*              Compute A=L*Q, copying result to VT
*              (CWorkspace: need 2*M, prefer M+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGELQF', M, N, 0, 0, NB )
               CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
*
*              Generate Q in VT
*              (CWorkspace: need M+N, prefer M+N*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGLQ', ' ', N, N, M, -1 )
               OPS = OPS + DOPLA( 'ZUNGLQ', N, N, M, 0, NB )
               CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Produce L in A, zeroing out above it
*
               CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
     $                      LDA )
               IE = 1
               ITAUQ = ITAU
               ITAUP = ITAUQ + M
               NWORK = ITAUP + M
*
*              Bidiagonalize L in A
*              (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
*              (RWorkspace: need M)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', M, M, 0, 0, NB )
               CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
*
*              Perform bidiagonal SVD, computing left singular vectors
*              of bidiagonal matrix in RWORK(IRU) and computing right
*              singular vectors of bidiagonal matrix in RWORK(IRVT)
*              (CWorkspace: need 0)
*              (RWorkspace: need BDSPAC)
*
               IRU = IE + M
               IRVT = IRU + M*M
               NRWORK = IRVT + M*M
               CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
     $                      M, RWO