sgesdd.f.html

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

HTML
742 行
字号
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 SGEMM( <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="SLACPY.586"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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="SGEQRF.607"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</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="SLACPY.612"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="SLASET.613"></a><a href="slaset.f.html#SLASET.1">SLASET</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="SORGQR.619"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</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="SGEBRD.629"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</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="SBDSDC.638"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</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="SORMBR.646"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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="SORMBR.650"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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="SLACPY.658"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'F'</span>, N, N, U, LDU, WORK( IR ), LDWRKR )
               CALL SGEMM( <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="SGEQRF.679"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL <a name="SLACPY.681"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</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="SORGQR.685"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</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="SLASET.690"></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 )
               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="SGEBRD.699"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</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="SBDSDC.708"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</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="SORMBR.716"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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="SORMBR.719"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</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">*</span><span class="comment">
</span><span class="comment
12345

sgesdd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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