zgesdd.f.html

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

               MINWRK = 2*M + N
               IF( WNTQO ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.414"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.414"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'P'</span>, M, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.416"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.416"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'Q'</span>, M, M, N, -1 ) )
                  MAXWRK = MAXWRK + M*N
                  MINWRK = MINWRK + M*M
               ELSE IF( WNTQS ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.421"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.421"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'P'</span>, M, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.423"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.423"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'Q'</span>, M, M, N, -1 ) )
               ELSE IF( WNTQA ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+N*
     $                     <a name="ILAENV.426"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.426"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'P'</span>, N, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.428"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.428"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'Q'</span>, M, M, N, -1 ) )
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 6t (N greater than M, but not much larger)
</span><span class="comment">*</span><span class="comment">
</span>               MAXWRK = 2*M + ( M+N )*<a name="ILAENV.434"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZGEBRD.434"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>'</span>, <span class="string">' '</span>, M, N,
     $                  -1, -1 )
               MINWRK = 2*M + N
               IF( WNTQO ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.439"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMBR.439"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>'</span>, <span class="string">'PRC'</span>, M, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.441"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMBR.441"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>'</span>, <span class="string">'QLN'</span>, M, M, N, -1 ) )
                  MAXWRK = MAXWRK + M*N
                  MINWRK = MINWRK + M*M
               ELSE IF( WNTQS ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.446"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.446"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'PRC'</span>, M, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.448"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.448"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'QLN'</span>, M, M, N, -1 ) )
               ELSE IF( WNTQA ) THEN
                  MAXWRK = MAX( MAXWRK, 2*M+N*
     $                     <a name="ILAENV.451"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.451"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'PRC'</span>, N, N, M, -1 ) )
                  MAXWRK = MAX( MAXWRK, 2*M+M*
     $                     <a name="ILAENV.453"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNGBR.453"></a><a href="zungbr.f.html#ZUNGBR.1">ZUNGBR</a>'</span>, <span class="string">'QLN'</span>, M, M, N, -1 ) )
               END IF
            END IF
         END IF
         MAXWRK = MAX( MAXWRK, MINWRK )
      END IF
      IF( INFO.EQ.0 ) THEN
         WORK( 1 ) = MAXWRK
         IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
     $      INFO = -13
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick returns
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.468"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGESDD.468"></a><a href="zgesdd.f.html#ZGESDD.1">ZGESDD</a>'</span>, -INFO )
         RETURN
      END IF
      IF( LWORK.EQ.LQUERV )
     $   RETURN
      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.479"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = SQRT( <a name="DLAMCH.480"></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.485"></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.489"></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.492"></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.MNTHR1 ) THEN
<span class="comment">*</span><span class="comment">
</span>            IF( WNTQN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 1 (M much larger than N, JOBZ='N')
</span><span class="comment">*</span><span class="comment">              No singular vectors to be computed
</span><span class="comment">*</span><span class="comment">
</span>               ITAU = 1
               NWORK = 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.515"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+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.520"></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
               NWORK = 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.531"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
               NRWORK = IE + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Perform bidiagonal SVD, compute singular values only
</span><span class="comment">*</span><span class="comment">              (CWorkspace: 0)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: need BDSPAN)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DBDSDC.540"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>( <span class="string">'U'</span>, <span class="string">'N'</span>, N, S, RWORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span>            ELSE IF( WNTQO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 2 (M much larger than N, JOBZ='O')
</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">              N right singular vectors to be computed in VT
</span><span class="comment">*</span><span class="comment">
</span>               IU = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              WORK(IU) is N by N
</span><span class="comment">*</span><span class="comment">
</span>               LDWRKU = N