zgesvd.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 669 行 · 第 1/5 页

     $            MAXWRK = MAX( MAXWRK, 2*M+N*
     $                     <a name="ILAENV.519"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.519"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'P'</span>, N, N, M, -1 ) )
               IF( .NOT.WNTUN )
     $            MAXWRK = MAX( MAXWRK, 2*M+( M-1 )*
     $                     <a name="ILAENV.522"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.522"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'Q'</span>, M, M, M, -1 ) )
               MINWRK = 2*M + N
            END IF
         END IF
         MAXWRK = MAX( MAXWRK, MINWRK )
         WORK( 1 ) = MAXWRK
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -13
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.535"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGESVD.535"></a><a href="zgesvd.f.html#ZGESVD.1">ZGESVD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.549"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = SQRT( <a name="DLAMCH.550"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) ) / EPS
      BIGNUM = ONE / SMLNUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale A if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="ZLANGE.555"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, DUM )
      ISCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
         ISCL = 1
         CALL <a name="ZLASCL.559"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
      ELSE IF( ANRM.GT.BIGNUM ) THEN
         ISCL = 1
         CALL <a name="ZLASCL.562"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( M.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A has at least as many rows as columns. If A has sufficiently
</span><span class="comment">*</span><span class="comment">        more rows than columns, first reduce using the QR
</span><span class="comment">*</span><span class="comment">        decomposition (if sufficient workspace available)
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.GE.MNTHR ) THEN
<span class="comment">*</span><span class="comment">
</span>            IF( WNTUN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 1 (M much larger than N, JOBU='N')
</span><span class="comment">*</span><span class="comment">              No left singular vectors to be computed
</span><span class="comment">*</span><span class="comment">
</span>               ITAU = 1
               IWORK = ITAU + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute A=Q*R
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: need 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGEQRF.585"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
     $                      LWORK-IWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Zero out below R
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLASET.590"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'L'</span>, N-1, N-1, CZERO, CZERO, A( 2, 1 ),
     $                      LDA )
               IE = 1
               ITAUQ = 1
               ITAUP = ITAUQ + N
               IWORK = ITAUP + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Bidiagonalize R in A
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGEBRD.601"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
     $                      IERR )
               NCVT = 0
               IF( WNTVO .OR. WNTVAS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 If right singular vectors desired, generate P'.
</span><span class="comment">*</span><span class="comment">                 (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">                 (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="ZUNGBR.611"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>( <span class="string">'P'</span>, N, N, N, A, LDA, WORK( ITAUP ),
     $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
                  NCVT = N
               END IF
               IRWORK = IE + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Perform bidiagonal QR iteration, computing right
</span><span class="comment">*</span><span class="comment">              singular vectors of A in A if desired
</span><span class="comment">*</span><span class="comment">              (CWorkspace: 0)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZBDSQR.622"></a><a href="zbdsqr.f.html#ZBDSQR.1">ZBDSQR</a>( <span class="string">'U'</span>, N, NCVT, 0, 0, S, RWORK( IE ), A, LDA,
     $                      CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              If right singular vectors desired in VT, copy them there
</span><span class="comment">*</span><span class="comment">
</span>               IF( WNTVAS )
     $            CALL <a name="ZLACPY.628"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'F'</span>, N, N, A, LDA, VT, LDVT )
<span class="comment">*</span><span class="comment">
</span>            ELSE IF( WNTUO .AND. WNTVN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
</span><span class="comment">*</span><span class="comment">              N left singular vectors to be overwritten on A and
</span><span class="comment">*</span><span class="comment">              no right singular vectors to be computed
</span><span class="comment">*</span><span class="comment">
</span>               IF( LWORK.GE.N*N+3*N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Sufficient workspace for a fast algorithm
</span><span class="comment">*</span><span class="comment">
</span>                  IR = 1
                  IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
</span><span class="comment">*</span><span class="comment">
</span>                     LDWRKU = LDA
                     LDWRKR = LDA
                  ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    WORK(IU) is LDA by N, WORK(IR) is N by N
</span>