cgelss.f.html

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

HTML
659 行
字号
1234
</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">           (CWorkspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">           (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CGEQRF.322"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</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">           (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">           (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CUNMQR.329"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>( <span class="string">'L'</span>, <span class="string">'C'</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="CLASET.335"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'L'</span>, N-1, N-1, CZERO, CZERO, A( 2, 1 ),
     $                      LDA )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IE = 1
         ITAUQ = 1
         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">        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CGEBRD.348"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</a>( MM, N, A, LDA, S, RWORK( 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">        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNMBR.356"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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( 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">        (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNGBR.363"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</a>( <span class="string">'P'</span>, N, N, N, A, LDA, WORK( ITAUP ),
     $                WORK( IWORK ), LWORK-IWORK+1, INFO )
         IRWORK = 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">        (CWorkspace: none)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CBDSQR.373"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a>( <span class="string">'U'</span>, N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
     $                1, B, LDB, RWORK( IRWORK ), 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="CSRSCL.386"></a><a href="csrscl.f.html#CSRSCL.1">CSRSCL</a>( NRHS, S( I ), B( I, 1 ), LDB )
               RANK = RANK + 1
            ELSE
               CALL <a name="CLASET.389"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, 1, NRHS, CZERO, CZERO, 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">        (CWorkspace: need N, prefer N*NRHS)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
            CALL CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, N, NRHS, N, CONE, A, LDA, B, LDB,
     $                  CZERO, WORK, LDB )
            CALL <a name="CLACPY.400"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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 CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, N, BL, N, CONE, A, LDA, B( 1, I ),
     $                     LDB, CZERO, WORK, N )
               CALL <a name="CLACPY.407"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'G'</span>, N, BL, WORK, N, B( 1, I ), LDB )
   20       CONTINUE
         ELSE
            CALL CGEMV( <span class="string">'C'</span>, N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
            CALL CCOPY( N, WORK, 1, B, 1 )
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
     $          THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Underdetermined case, M much less than N
</span><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.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
     $      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">        (CWorkspace: need 2*M, prefer M+M*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CGELQF.432"></a><a href="cgelqf.f.html#CGELQF.1">CGELQF</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="CLACPY.438"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'L'</span>, M, M, A, LDA, WORK( IL ), LDWORK )
         CALL <a name="CLASET.439"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'U'</span>, M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
     $                LDWORK )
         IE = 1
         ITAUQ = IL + LDWORK*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">        (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need M)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CGEBRD.450"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</a>( M, M, WORK( IL ), LDWORK, S, RWORK( 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">        (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNMBR.458"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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( 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">        (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: none)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNGBR.466"></a><a href="cungbr.f.html#CUNGBR.1">CUNGBR</a>( <span class="string">'P'</span>, M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
     $                WORK( IWORK ), LWORK-IWORK+1, INFO )
         IRWORK = IE + M
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Perform bidiagonal QR iteration, computing right singular
</span><span class="comment">*</span><span class="comment">        vectors of L in WORK(IL) and multiplying B by transpose of
</span><span class="comment">*</span><span class="comment">        left singular vectors
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need M*M)
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CBDSQR.476"></a><a href="cbdsqr.f.html#CBDSQR.1">CBDSQR</a>( <span class="string">'U'</span>, M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
1234

cgelss.f.html - 源码说明

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

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

⌨️ 快捷键说明

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