zgesdd.f

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

*              Multiply Q in A by left singular vectors of R in
*              WORK(IU), storing result in WORK(IR) and copying to A
*              (CWorkspace: need 2*N*N, prefer N*N+M*N)
*              (RWorkspace: 0)
*
               DO 10 I = 1, M, LDWRKR
                  CHUNK = MIN( M-I+1, LDWRKR )
                  OPS = OPS + DOPBL3( 'ZGEMM ', CHUNK, N, N )
                  CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
     $                        LDA, WORK( IU ), LDWRKU, CZERO,
     $                        WORK( IR ), LDWRKR )
                  CALL ZLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
     $                         A( I, 1 ), LDA )
   10          CONTINUE
*
            ELSE IF( WNTQS ) THEN
*
*              Path 3 (M much larger than N, JOBZ='S')
*              N left singular vectors to be computed in U and
*              N right singular vectors to be computed in VT
*
               IR = 1
*
*              WORK(IR) is N by N
*
               LDWRKR = N
               ITAU = IR + LDWRKR*N
               NWORK = ITAU + N
*
*              Compute A=Q*R
*              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEQRF', M, N, 0, 0, NB )
               CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy R to WORK(IR), zeroing out below it
*
               CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
     $                      LDWRKR )
*
*              Generate Q in A
*              (CWorkspace: need 2*N, prefer N+N*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGQR', ' ', M, N, N, -1 )
               OPS = OPS + DOPLA( 'ZUNGQR', M, N, N, 0, NB )
               CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = 1
               ITAUQ = ITAU
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in WORK(IR)
*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
*              (RWorkspace: need N)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', N, N, 0, 0, NB )
               CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, 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 + N
               IRVT = IRU + N*N
               NRWORK = IRVT + N*N
               CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
     $                      N, RWORK( IRVT ), N, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
*
*              Copy real matrix RWORK(IRU) to complex matrix U
*              Overwrite U by left singular vectors of R
*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', N, N, N, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Copy real matrix RWORK(IRVT) to complex matrix VT
*              Overwrite VT by right singular vectors of R
*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
               NB = ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRC', N, N, N, 0, NB )
               CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply Q in A by left singular vectors of R in
*              WORK(IR), storing result in U
*              (CWorkspace: need N*N)
*              (RWorkspace: 0)
*
               CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
               OPS = OPS + DOPBL3( 'ZGEMM ', M, N, N )
               CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
     $                     LDWRKR, CZERO, U, LDU )
*
            ELSE IF( WNTQA ) THEN
*
*              Path 4 (M much larger than N, JOBZ='A')
*              M left singular vectors to be computed in U and
*              N right singular vectors to be computed in VT
*
               IU = 1
*
*              WORK(IU) is N by N
*
               LDWRKU = N
               ITAU = IU + LDWRKU*N
               NWORK = ITAU + N
*
*              Compute A=Q*R, copying result to U
*              (CWorkspace: need 2*N, prefer N+N*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEQRF', M, N, 0, 0, NB )
               CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
*
*              Generate Q in U
*              (CWorkspace: need N+M, prefer N+M*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGQR', ' ', M, M, N, -1 )
               OPS = OPS + DOPLA( 'ZUNGQR', M, M, N, 0, NB )
               CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Produce R in A, zeroing out below it
*
               CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
     $                      LDA )
               IE = 1
               ITAUQ = ITAU
               ITAUP = ITAUQ + N
               NWORK = ITAUP + N
*
*              Bidiagonalize R in A
*              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
*              (RWorkspace: need N)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', N, N, 0, 0, NB )
               CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
               IRU = IE + N
               IRVT = IRU + N*N
               NRWORK = IRVT + N*N
*
*              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)
*
               CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
     $                      N, RWORK( IRVT ), N, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
*
*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
*              Overwrite WORK(IU) by left singular vectors of R
*              (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
     $                      LDWRKU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', N, N, N, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Copy real matrix RWORK(IRVT) to complex matrix VT
*              Overwrite VT by right singular vectors of R
*              (CWorkspace: need 3*N, prefer 2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
               NB = ILAENV( 1, 'ZUNMBR', 'PRT', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRT', N, N, N, 0, NB )
               CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
*              Multiply Q in U by left singular vectors of R in
*              WORK(IU), storing result in A
*              (CWorkspace: need N*N)
*              (RWorkspace: 0)
*
               OPS = OPS + DOPBL3( 'ZGEMM ', M, N, N )
               CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
     $                     LDWRKU, CZERO, A, LDA )
*
*              Copy left singular vectors of A from A to U
*
               CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
*
            END IF
*
         ELSE IF( M.GE.MNTHR2 ) THEN
*
*           MNTHR2 <= M < MNTHR1
*
*           Path 5 (M much larger than N, but not as much as MNTHR1)
*           Reduce to bidiagonal form without QR decomposition, use
*           ZUNGBR and matrix multiplication to compute singular vectors
*
            IE = 1
            NRWORK = IE + N
            ITAUQ = 1
            ITAUP = ITAUQ + N
            NWORK = ITAUP + N
*
*           Bidiagonalize A
*           (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
*           (RWorkspace: need N)
*
               NB = ILAENV( 1, 'ZGEBRD', ' ', M, N, -1, -1 )
               OPS = OPS + DOPLA( 'ZGEBRD', M, N, 0, 0, NB )
            CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                   IERR )
            IF( WNTQN ) THEN
*
*              Compute singular values only
*              (Cworkspace: 0)
*              (Rworkspace: need BDSPAC)
*
               CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
            ELSE IF( WNTQO ) THEN
               IU = NWORK
               IRU = NRWORK
               IRVT = IRU + N*N
               NRWORK = IRVT + N*N
*
*              Copy A to VT, generate P**H
*              (Cworkspace: need 2*N, prefer N+N*NB)
*              (Rworkspace: 0)
*
               CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
               NB = ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNGBR', 'P', N, N, N, 0, NB )
               CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*              Generate Q in A
*              (CWorkspace: need 2*N, prefer N+N*NB)
*              (RWorkspace: 0)
*
               NB = ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNGBR', 'Q', N, N, N, 0, NB )
               CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
*
               IF( LWORK.GE.M*N+3*N ) THEN
*
*                 WORK( IU ) is M by N
*
                  LDWRKU = M
               ELSE
*
*                 WORK(IU) is LDWRKU by N
*
                  LDWRKU = ( LWORK-3*N ) / N
               END IF
               NWORK = IU + LDWRKU*N
*
*              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)
*
               CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
     $                      N, RWORK( IRVT ), N, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
*
*              Multiply real matrix RWORK(IRVT) by P**H in VT,
*              storing the result in WORK(IU), copying to VT
*              (Cworkspace: need 0)
*              (Rworkspace: need 3*N*N)
*
               OPS = OPS + DBLE( 4*N*N*N )
               CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
     $                      WORK( IU ), LDWRKU, RWORK( NRWORK ) )
               CALL ZLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
*
*              Multiply Q in A by real matrix RWORK(IRU), storing the
*              result in WORK(IU), copying to A
*              (CWorkspace: need N*N, prefer M*N)
*              (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
*
               NRWORK = IRVT
               DO 20 I = 1, M, LDWRKU
                  CHUNK = MIN( M-I+1, LDWRKU )
                  OPS = OPS + DBLE( 4*CHUNK*N*N )
                  CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
     $                         N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
                  CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
     $                         A( I, 1 ), LDA )
   20          CONTINUE
*
            ELSE IF( WNTQS ) THEN