zgelsd.f.html

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

         MINWRK = MIN( MINWRK, MAXWRK )
         WORK( 1 ) = MAXWRK
         IWORK( 1 ) = LIWORK
         RWORK( 1 ) = LRWORK
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -12
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.289"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGELSD.289"></a><a href="zgelsd.f.html#ZGELSD.1">ZGELSD</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
         RANK = 0
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine parameters.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="DLAMCH.304"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SFMIN = <a name="DLAMCH.305"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SMLNUM = SFMIN / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="DLABAD.308"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale A if max entry outside range [SMLNUM,BIGNUM].
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="ZLANGE.312"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, RWORK )
      IASCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASCL.318"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
         IASCL = 1
      ELSE IF( ANRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASCL.324"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
         IASCL = 2
      ELSE IF( ANRM.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Matrix all zero. Return zero solution.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASET.330"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
         CALL <a name="DLASET.331"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'F'</span>, MINMN, 1, ZERO, ZERO, S, 1 )
         RANK = 0
         GO TO 10
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale B if max entry outside range [SMLNUM,BIGNUM].
</span><span class="comment">*</span><span class="comment">
</span>      BNRM = <a name="ZLANGE.338"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>( <span class="string">'M'</span>, M, NRHS, B, LDB, RWORK )
      IBSCL = 0
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASCL.344"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 1
      ELSE IF( BNRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASCL.350"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 2
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If M &lt; N make sure B(M+1:N,:) = 0
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LT.N )
     $   CALL <a name="ZLASET.357"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'F'</span>, N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Overdetermined case.
</span><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">        Path 1 - overdetermined or exactly determined.
</span><span class="comment">*</span><span class="comment">
</span>         MM = M
         IF( M.GE.MNTHR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Path 1a - overdetermined, with many more rows than columns
</span><span class="comment">*</span><span class="comment">
</span>            MM = N
            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">           (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">           (CWorkspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="ZGEQRF.378"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                   LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply B by transpose(Q).
</span><span class="comment">*</span><span class="comment">           (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">           (CWorkspace: need NRHS, prefer NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="ZUNMQR.385"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</a>( <span class="string">'L'</span>, <span class="string">'C'</span>, M, NRHS, N, A, LDA, WORK( ITAU ), B,
     $                   LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<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>            IF( N.GT.1 ) THEN
               CALL <a name="ZLASET.391"></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 )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         ITAUQ = 1
         ITAUP = ITAUQ + N
         NWORK = ITAUP + N
         IE = 1
         NRWORK = IE + 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">        (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGEBRD.406"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by transpose of left bidiagonalizing vectors of R.
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.413"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'C'</span>, MM, NRHS, N, A, LDA, WORK( ITAUQ ),
     $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve the bidiagonal least squares problem.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLALSD.418"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>( <span class="string">'U'</span>, SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
     $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
     $                IWORK, INFO )
         IF( INFO.NE.0 ) THEN
            GO TO 10
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by right bidiagonalizing vectors of R.