cgesdd.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="CUNGBR.414"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.416"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.421"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.423"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.426"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.428"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CGEBRD.434"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</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="CUNMBR.439"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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="CUNMBR.441"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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="CUNGBR.446"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.448"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.451"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CUNGBR.453"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</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="CGESDD.468"></a><a href="cgesdd.f.html#CGESDD.1">CGESDD</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="SLAMCH.479"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = SQRT( <a name="SLAMCH.480"></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="CLANGE.485"></a><a href="clange.f.html#CLANGE.1">CLANGE</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="CLASCL.489"></a><a href="clascl.f.html#CLASCL.1">CLASCL</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="CLASCL.492"></a><a href="clascl.f.html#CLASCL.1">CLASCL</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="CGEQRF.515"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</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="CLASET.520"></a><a href="claset.f.html#CLASET.1">CLASET</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="CGEBRD.531"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</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="SBDSDC.540"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</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