dgelsd.f.html

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

HTML
553 行
字号
123
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DGEQRF.353"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</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">           (Workspace: need N+NRHS, prefer N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DORMQR.359"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>( <span class="string">'L'</span>, <span class="string">'T'</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="DLASET.365"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'L'</span>, N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IE = 1
         ITAUQ = IE + N
         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">        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGEBRD.377"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>( MM, 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 of R.
</span><span class="comment">*</span><span class="comment">        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DORMBR.384"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'T'</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="DLALSD.389"></a><a href="dlalsd.f.html#DLALSD.1">DLALSD</a>( <span class="string">'U'</span>, SMLSIZ, N, 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 R.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DORMBR.397"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</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="DGELQF.415"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</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="DLACPY.421"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'L'</span>, M, M, A, LDA, WORK( IL ), LDWORK )
         CALL <a name="DLASET.422"></a><a href="dlaset.f.html#DLASET.1">DLASET</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="DGEBRD.432"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</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="DORMBR.439"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</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="DLALSD.445"></a><a href="dlalsd.f.html#DLALSD.1">DLALSD</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="DORMBR.453"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</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="DLASET.459"></a><a href="dlaset.f.html#DLASET.1">DLASET</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="DORMLQ.465"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</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="DGEBRD.480"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</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="DORMBR.487"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</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="DLALSD.492"></a><a href="dlalsd.f.html#DLALSD.1">DLALSD</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="DORMBR.500"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</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="DLASCL.508"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="DLASCL.509"></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="DLASCL.512"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
         CALL <a name="DLASCL.513"></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="DLASCL.517"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         CALL <a name="DLASCL.519"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DGELSD.526"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
123

dgelsd.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 dgelsd.f.html 源码文件，采用 HTML 编程语言编写，共 553 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?