zgesdd.f

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

*
*              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 )
*
*              Copy A to U, generate Q
*              (Cworkspace: need 2*N, prefer N+N*NB)
*              (Rworkspace: 0)
*
               CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
               NB = ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNGBR', 'Q', M, N, N, 0, NB )
               CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
     $                      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 = NRWORK
               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 )
*
*              Multiply real matrix RWORK(IRVT) by P**H in VT,
*              storing the result in A, 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, A, LDA,
     $                      RWORK( NRWORK ) )
               CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
*
*              Multiply Q in U by real matrix RWORK(IRU), storing the
*              result in A, copying to U
*              (CWorkspace: need 0)
*              (Rworkspace: need N*N+2*M*N)
*
               NRWORK = IRVT
               CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
     $                      RWORK( NRWORK ) )
               CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
            ELSE
*
*              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 )
*
*              Copy A to U, generate Q
*              (Cworkspace: need 2*N, prefer N+N*NB)
*              (Rworkspace: 0)
*
               CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
               NB = ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNGBR', 'Q', M, M, N, 0, NB )
               CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
     $                      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 = NRWORK
               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 )
*
*              Multiply real matrix RWORK(IRVT) by P**H in VT,
*              storing the result in A, 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, A, LDA,
     $                      RWORK( NRWORK ) )
               CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
*
*              Multiply Q in U by real matrix RWORK(IRU), storing the
*              result in A, copying to U
*              (CWorkspace: 0)
*              (Rworkspace: need 3*N*N)
*
               NRWORK = IRVT
               CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
     $                      RWORK( NRWORK ) )
               CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
            END IF
*
         ELSE
*
*           M .LT. MNTHR2
*
*           Path 6 (M at least N, but not much larger)
*           Reduce to bidiagonal form without QR decomposition
*           Use ZUNMBR 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
               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 )
*
*              Copy real matrix RWORK(IRVT) to complex matrix VT
*              Overwrite VT by right singular vectors of A
*              (Cworkspace: need 2*N, prefer N+N*NB)
*              (Rworkspace: need 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, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
*
               IF( LWORK.GE.M*N+3*N ) THEN
*
*              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
*              Overwrite WORK(IU) by left singular vectors of A, copying
*              to A
*              (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)
*              (Rworkspace: need 0)
*
                  CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
     $                         LDWRKU )
                  CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
     $                         LDWRKU )
                  NB = ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 )
                  OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', M, N, N, 0, NB )
                  CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
     $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
                  CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
               ELSE
*
*                 Generate Q in A
*                 (Cworkspace: need 2*N, prefer N+N*NB)
*                 (Rworkspace: need 0)
*
                  NB = ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 )
                  OPS = OPS + DOPLA2( 'ZUNGBR', 'Q', M, N, N, 0, NB )
                  CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
     $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
*
*                 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 30 I = 1, M, LDWRKU
                     CHUNK = MIN( M-I+1, LDWRKU )
                     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 )
   30             CONTINUE
               END IF
*
            ELSE IF( WNTQS ) THEN
*
*              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 = NRWORK
               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 A
*              (CWorkspace: need 3*N, prefer 2*N+N*NB)
*              (RWorkspace: 0)
*
               CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )
               CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', M, N, N, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
     $                      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 A
*              (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', 'PRC', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRC', 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 )
            ELSE
*
*              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 = NRWORK
               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 )
*
*              Set the right corner of U to identity matrix
*
               CALL ZLASET( 'F', M, M, CZERO, CZERO, U, LDU )
               CALL ZLASET( 'F', M-N, M-N, CZERO, CONE, U( N+1, N+1 ),
     $                      LDU )
*
*              Copy real matrix RWORK(IRU) to complex matrix U
*              Overwrite U by left singular vectors of A
*              (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
*              (RWorkspace: 0)
*
               CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
               NB = ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'QLN', M, M, N, 0, NB )
               CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
     $                      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 A
*              (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', 'PRC', N, N, N, -1 )
               OPS = OPS + DOPLA2( 'ZUNMBR', 'PRC', 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 )
            END IF
*
         END IF
*
      ELSE
*
*        A has more columns than rows. If A has sufficiently more
*        columns than rows, first reduce using the LQ decomposition
*        (if sufficient workspace available)
*
         IF( N.GE.MNTHR1 ) THEN
*
            IF( WNTQN ) THEN
*
*              Path 1t (N much larger than M, JOBZ='N')
*              No singular vectors to be computed