sgelss.f.html

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

HTML
642 行
字号
1234
         IBSCL = 2
      END IF
<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
            IWORK = 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">           (Workspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SGEQRF.325"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
     $                   LWORK-IWORK+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="SORMQR.331"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>( <span class="string">'L'</span>, <span class="string">'T'</span>, M, NRHS, N, A, LDA, WORK( ITAU ), B,
     $                   LDB, WORK( IWORK ), LWORK-IWORK+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 )
     $         CALL <a name="SLASET.337"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'L'</span>, N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IE = 1
         ITAUQ = IE + N
         ITAUP = ITAUQ + N
         IWORK = 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="SGEBRD.348"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+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="SORMBR.355"></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>, MM, NRHS, N, A, LDA, WORK( ITAUQ ),
     $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate right bidiagonalizing vectors of R in A
</span><span class="comment">*</span><span class="comment">        (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORGBR.361"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>( <span class="string">'P'</span>, N, N, N, A, LDA, WORK( ITAUP ),
     $                WORK( IWORK ), LWORK-IWORK+1, INFO )
         IWORK = IE + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Perform bidiagonal QR iteration
</span><span class="comment">*</span><span class="comment">          multiply B by transpose of left singular vectors
</span><span class="comment">*</span><span class="comment">          compute right singular vectors in A
</span><span class="comment">*</span><span class="comment">        (Workspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SBDSQR.370"></a><a href="sbdsqr.f.html#SBDSQR.1">SBDSQR</a>( <span class="string">'U'</span>, N, N, 0, NRHS, S, WORK( IE ), A, LDA, VDUM,
     $                1, B, LDB, WORK( IWORK ), INFO )
         IF( INFO.NE.0 )
     $      GO TO 70
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by reciprocals of singular values
</span><span class="comment">*</span><span class="comment">
</span>         THR = MAX( RCOND*S( 1 ), SFMIN )
         IF( RCOND.LT.ZERO )
     $      THR = MAX( EPS*S( 1 ), SFMIN )
         RANK = 0
         DO 10 I = 1, N
            IF( S( I ).GT.THR ) THEN
               CALL <a name="SRSCL.383"></a><a href="srscl.f.html#SRSCL.1">SRSCL</a>( NRHS, S( I ), B( I, 1 ), LDB )
               RANK = RANK + 1
            ELSE
               CALL <a name="SLASET.386"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
            END IF
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by right singular vectors
</span><span class="comment">*</span><span class="comment">        (Workspace: need N, prefer N*NRHS)
</span><span class="comment">*</span><span class="comment">
</span>         IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, A, LDA, B, LDB, ZERO,
     $                  WORK, LDB )
            CALL <a name="SLACPY.396"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'G'</span>, N, NRHS, WORK, LDB, B, LDB )
         ELSE IF( NRHS.GT.1 ) THEN
            CHUNK = LWORK / N
            DO 20 I = 1, NRHS, CHUNK
               BL = MIN( NRHS-I+1, CHUNK )
               CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, BL, N, ONE, A, LDA, B( 1, I ),
     $                     LDB, ZERO, WORK, N )
               CALL <a name="SLACPY.403"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'G'</span>, N, BL, WORK, N, B( 1, I ), LDB )
   20       CONTINUE
         ELSE
            CALL SGEMV( <span class="string">'T'</span>, N, N, ONE, A, LDA, B, 1, ZERO, WORK, 1 )
            CALL SCOPY( N, WORK, 1, B, 1 )
         END IF
<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
         IWORK = 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.425"></a><a href="sgelqf.f.html#SGELQF.1">SGELQF</a>( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
     $                LWORK-IWORK+1, INFO )
         IL = IWORK
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy L to WORK(IL), zeroing out above it
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLACPY.431"></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.432"></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
         IWORK = 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.442"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
     $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
     $                LWORK-IWORK+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.449"></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( IWORK ),
     $                LWORK-IWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate right bidiagonalizing vectors of R in WORK(IL)
</span><span class="comment">*</span><span class="comment">        (Workspace: need M*M+5*M-1, prefer M*M+4*M+(M-1)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORGBR.456"></a><a href="sorgbr.f.html#SORGBR.1">SORGBR</a>( <span class="string">'P'</span>, M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
     $                WORK( IWORK ), LWORK-IWORK+1, INFO )
         IWORK = IE + M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Perform bidiagonal QR iteration,
</span><span class="comment">*</span><span class="comment">           computing right singular vectors of L in WORK(IL) and
</span><span class="comment">*</span><span class="comment">           multiplying B by transpose of left singular vectors
</span><span class="comment">*</span><span class="comment">        (Workspace: need M*M+M+BDSPAC)
</span><span class="comment">*</span><span class="comment">
1234

sgelss.f.html - 源码说明

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

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

⌨️ 快捷键说明

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