sgesvd.f.html

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

            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 10t(N greater than M, but not much larger)
</span><span class="comment">*</span><span class="comment">
</span>               MAXWRK = 3*M + ( M+N )*<a name="ILAENV.526"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SGEBRD.526"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>'</span>, <span class="string">' '</span>, M, N,
     $                  -1, -1 )
               IF( WNTVS .OR. WNTVO )
     $            MAXWRK = MAX( MAXWRK, 3*M+M*
     $                     <a name="ILAENV.530"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SORGBR.530"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>'</span>, <span class="string">'P'</span>, M, N, M, -1 ) )
               IF( WNTVA )
     $            MAXWRK = MAX( MAXWRK, 3*M+N*
     $                     <a name="ILAENV.533"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SORGBR.533"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>'</span>, <span class="string">'P'</span>, N, N, M, -1 ) )
               IF( .NOT.WNTUN )
     $            MAXWRK = MAX( MAXWRK, 3*M+( M-1 )*
     $                     <a name="ILAENV.536"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="SORGBR.536"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>'</span>, <span class="string">'Q'</span>, M, M, M, -1 ) )
               MAXWRK = MAX( MAXWRK, BDSPAC )
               MINWRK = MAX( 3*M+N, BDSPAC )
            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.550"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGESVD.550"></a><a href="sgesvd.f.html#SGESVD.1">SGESVD</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="SLAMCH.564"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = SQRT( <a name="SLAMCH.565"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="SLANGE.570"></a><a href="slange.f.html#SLANGE.1">SLANGE</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="SLASCL.574"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SLASCL.577"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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">              (Workspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SGEQRF.599"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</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="SLASET.604"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'L'</span>, N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
               IE = 1
               ITAUQ = IE + N
               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">              (Workspace: need 4*N, prefer 3*N+2*N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SGEBRD.613"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( N, N, A, LDA, S, WORK( 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">                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SORGBR.622"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>( <span class="string">'P'</span>, N, N, N, A, LDA, WORK( ITAUP ),
     $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
                  NCVT = N
               END IF
               IWORK = 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">              (Workspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SBDSQR.632"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a>( <span class="string">'U'</span>, N, NCVT, 0, 0, S, WORK( IE ), A, LDA,
     $                      DUM, 1, DUM, 1, WORK( IWORK ), 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="SLACPY.638"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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+MAX( 4*N, BDSPAC ) ) 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+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
                     LD