sgelsd.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORMBR.406"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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, WLALSD ) ) 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, 4*M+M*LDA+WLALSD ) )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">        (Workspace: need 2*M, prefer M+M*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGELQF.424"></a><a href="sgelqf.f.html#SGELQF.1">SGELQF</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="SLACPY.430"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'L'</span>, M, M, A, LDA, WORK( IL ), LDWORK )
         CALL <a name="SLASET.431"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'U'</span>, M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
     $                LDWORK )
         IE = IL + LDWORK*M
         ITAUQ = IE + M
         ITAUP = ITAUQ + M
         NWORK = ITAUP + 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">        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBRD.441"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( M, M, WORK( IL ), LDWORK, S, WORK( 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">        (Workspace: 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="SORMBR.448"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'T'</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="SLALSD.454"></a><a href="slalsd.f.html#SLALSD.1">SLALSD</a>( <span class="string">'U'</span>, SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
     $                RCOND, RANK, WORK( NWORK ), 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="SORMBR.462"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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="SLASET.468"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, N-M, NRHS, ZERO, ZERO, 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">        (Workspace: need M+NRHS, prefer M+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORMLQ.474"></a><a href="sormlq.f.html#SORMLQ.1">SORMLQ</a>( <span class="string">'L'</span>, <span class="string">'T'</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>         IE = 1
         ITAUQ = IE + M
         ITAUP = ITAUQ + M
         NWORK = ITAUP + 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">        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBRD.489"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( M, N, A, LDA, S, WORK( 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">        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORMBR.496"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'T'</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="SLALSD.501"></a><a href="slalsd.f.html#SLALSD.1">SLALSD</a>( <span class="string">'L'</span>, SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
     $                RCOND, RANK, WORK( NWORK ), 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="SORMBR.509"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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="SLASCL.517"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="SLASCL.518"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
     $                INFO )
      ELSE IF( IASCL.EQ.2 ) THEN
         CALL <a name="SLASCL.521"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="SLASCL.522"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SLASCL.526"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         CALL <a name="SLASCL.528"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGELSD.536"></a><a href="sgelsd.f.html#SGELSD.1">SGELSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>