zgelsd.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.427"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, N, NRHS, N, A, LDA, WORK( ITAUP ),
     $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
     $         MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Path 2a - underdetermined, with many more columns than rows
</span><span class="comment">*</span><span class="comment">        and sufficient workspace for an efficient algorithm.
</span><span class="comment">*</span><span class="comment">
</span>         LDWORK = M
         IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
     $       M*LDA+M+M*NRHS ) )LDWORK = LDA
         ITAU = 1
         NWORK = M + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute A=L*Q.
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*M, prefer M+M*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGELQF.445"></a><a href="zgelqf.f.html#ZGELQF.1">ZGELQF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                LWORK-NWORK+1, INFO )
         IL = NWORK
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy L to WORK(IL), zeroing out above its diagonal.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLACPY.451"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'L'</span>, M, M, A, LDA, WORK( IL ), LDWORK )
         CALL <a name="ZLASET.452"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'U'</span>, M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
     $                LDWORK )
         ITAUQ = IL + LDWORK*M
         ITAUP = ITAUQ + M
         NWORK = ITAUP + M
         IE = 1
         NRWORK = IE + M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Bidiagonalize L in WORK(IL).
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need M)
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGEBRD.464"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( M, M, WORK( IL ), LDWORK, 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 L.
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.471"></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>, M, NRHS, M, WORK( IL ), LDWORK,
     $                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.477"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>( <span class="string">'U'</span>, SMLSIZ, M, 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 L.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.486"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, M, NRHS, M, WORK( IL ), LDWORK,
     $                WORK( ITAUP ), 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 first M rows of B.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZLASET.492"></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 )
         NWORK = ITAU + M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply transpose(Q) by B.
</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="ZUNMLQ.498"></a><a href="zunmlq.f.html#ZUNMLQ.1">ZUNMLQ</a>( <span class="string">'L'</span>, <span class="string">'C'</span>, N, NRHS, M, A, LDA, WORK( ITAU ), B,
     $                LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Path 2 - remaining underdetermined cases.
</span><span class="comment">*</span><span class="comment">
</span>         ITAUQ = 1
         ITAUP = ITAUQ + M
         NWORK = ITAUP + M
         IE = 1
         NRWORK = IE + M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Bidiagonalize A.
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need M)
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGEBRD.515"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( M, 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.
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.522"></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>, M, 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.527"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>( <span class="string">'L'</span>, SMLSIZ, M, 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 A.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZUNMBR.536"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, N, NRHS, M, A, LDA, WORK( ITAUP ),
     $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling.
</span><span class="comment">*</span><span class="comment">
</span>      IF( IASCL.EQ.1 ) THEN
         CALL <a name="ZLASCL.544"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="DLASCL.545"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
     $                INFO )
      ELSE IF( IASCL.EQ.2 ) THEN
         CALL <a name="ZLASCL.548"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="DLASCL.549"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
     $                INFO )
      END IF
      IF( IBSCL.EQ.1 ) THEN
         CALL <a name="ZLASCL.553"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         CALL <a name="ZLASCL.555"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      WORK( 1 ) = MAXWRK
      IWORK( 1 ) = LIWORK
      RWORK( 1 ) = LRWORK
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZGELSD.564"></a><a href="zgelsd.f.html#ZGELSD.1">ZGELSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>