dgesdd.f.html

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

HTML
741 行
字号
12345
</span><span class="comment">*</span><span class="comment">              WORK(IU), storing result in WORK(IR) and copying to A
</span><span class="comment">*</span><span class="comment">              (Workspace: need 2*N*N, prefer N*N+M*N)
</span><span class="comment">*</span><span class="comment">
</span>               DO 10 I = 1, M, LDWRKR
                  CHUNK = MIN( M-I+1, LDWRKR )
                  CALL DGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, CHUNK, N, N, ONE, A( I, 1 ),
     $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
     $                        LDWRKR )
                  CALL <a name="DLACPY.586"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, CHUNK, N, WORK( IR ), LDWRKR,
     $                         A( I, 1 ), LDA )
   10          CONTINUE
<span class="comment">*</span><span class="comment">
</span>            ELSE IF( WNTQS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 3 (M much larger than N, JOBZ='S')
</span><span class="comment">*</span><span class="comment">              N left singular vectors to be computed in U and
</span><span class="comment">*</span><span class="comment">              N right singular vectors to be computed in VT
</span><span class="comment">*</span><span class="comment">
</span>               IR = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              WORK(IR) is N by N
</span><span class="comment">*</span><span class="comment">
</span>               LDWRKR = N
               ITAU = IR + LDWRKR*N
               NWORK = 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 N*N+2*N, prefer N*N+N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DGEQRF.607"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy R to WORK(IR), zeroing out below it
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLACPY.612"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="DLASET.613"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'L'</span>, N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
     $                      LDWRKR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Generate Q in A
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DORGQR.619"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = ITAU
               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 WORK(IR)
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DGEBRD.629"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
     $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Perform bidiagonal SVD, computing left singular vectors
</span><span class="comment">*</span><span class="comment">              of bidiagoal matrix in U and computing right singular
</span><span class="comment">*</span><span class="comment">              vectors of bidiagonal matrix in VT
</span><span class="comment">*</span><span class="comment">              (Workspace: need N+BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DBDSDC.638"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>( <span class="string">'U'</span>, <span class="string">'I'</span>, N, S, WORK( IE ), U, LDU, VT,
     $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Overwrite U by left singular vectors of R and VT
</span><span class="comment">*</span><span class="comment">              by right singular vectors of R
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DORMBR.646"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="DORMBR.650"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, <span class="string">'T'</span>, N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Multiply Q in A by left singular vectors of R in
</span><span class="comment">*</span><span class="comment">              WORK(IR), storing result in U
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLACPY.658"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'F'</span>, N, N, U, LDU, WORK( IR ), LDWRKR )
               CALL DGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, M, N, N, ONE, A, LDA, WORK( IR ),
     $                     LDWRKR, ZERO, U, LDU )
<span class="comment">*</span><span class="comment">
</span>            ELSE IF( WNTQA ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 4 (M much larger than N, JOBZ='A')
</span><span class="comment">*</span><span class="comment">              M left singular vectors to be computed in U and
</span><span class="comment">*</span><span class="comment">              N right singular vectors to be computed in VT
</span><span class="comment">*</span><span class="comment">
</span>               IU = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              WORK(IU) is N by N
</span><span class="comment">*</span><span class="comment">
</span>               LDWRKU = N
               ITAU = IU + LDWRKU*N
               NWORK = ITAU + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute A=Q*R, copying result to U
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DGEQRF.679"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL <a name="DLACPY.681"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'L'</span>, M, N, A, LDA, U, LDU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Generate Q in U
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
</span>               CALL <a name="DORGQR.685"></a><a href="dorgqr.f.html#DORGQR.1">DORGQR</a>( M, M, N, U, LDU, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Produce R in A, zeroing out other entries
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLASET.690"></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 )
               IE = ITAU
               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 N*N+4*N, prefer N*N+3*N+2*N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DGEBRD.699"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
     $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                      IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Perform bidiagonal SVD, computing left singular vectors
</span><span class="comment">*</span><span class="comment">              of bidiagonal matrix in WORK(IU) and computing right
</span><span class="comment">*</span><span class="comment">              singular vectors of bidiagonal matrix in VT
</span><span class="comment">*</span><span class="comment">              (Workspace: need N+N*N+BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DBDSDC.708"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>( <span class="string">'U'</span>, <span class="string">'I'</span>, N, S, WORK( IE ), WORK( IU ), N,
     $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
     $                      INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Overwrite WORK(IU) by left singular vectors of R and VT
</span><span class="comment">*</span><span class="comment">              by right singular vectors of R
</span><span class="comment">*</span><span class="comment">              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DORMBR.716"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, N, N, N, A, LDA,
     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               CALL <a name="DORMBR.719"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, <span class="string">'T'</span>, N, N, N, A, LDA,
     $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
<span class="comment">
12345

dgesdd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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