zgesdd.f.html

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

HTML
713 行
字号
12345
               IR = IU + LDWRKU*N
               IF( LWORK.GE.M*N+N*N+3*N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 WORK(IR) is M by N
</span><span class="comment">*</span><span class="comment">
</span>                  LDWRKR = M
               ELSE
                  LDWRKR = ( LWORK-N*N-3*N ) / N
               END IF
               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">              (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGEQRF.570"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</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="ZLACPY.575"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="ZLASET.576"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'L'</span>, N-1, N-1, CZERO, CZERO, 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">              (CWorkspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZUNGQR.583"></a><a href="zungqr.f.html#ZUNGQR.1">ZUNGQR</a>( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = 1
               ITAUQ = ITAU
               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">              (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*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="ZGEBRD.594"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( N, N, WORK( IR ), LDWRKR, S, RWORK( 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 R in WORK(IRU) and computing right singular vectors
</span><span class="comment">*</span><span class="comment">              of R in WORK(IRVT)
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need 0)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: need BDSPAC)
</span><span class="comment">*</span><span class="comment">
</span>               IRU = IE + N
               IRVT = IRU + N*N
               NRWORK = IRVT + N*N
               CALL <a name="DBDSDC.607"></a><a href="dbdsdc.f.html#DBDSDC.1">DBDSDC</a>( <span class="string">'U'</span>, <span class="string">'I'</span>, N, S, RWORK( IE ), RWORK( IRU ),
     $                      N, RWORK( IRVT ), N, DUM, IDUM,
     $                      RWORK( NRWORK ), IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
</span><span class="comment">*</span><span class="comment">              Overwrite WORK(IU) by the left singular vectors of R
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACP2.616"></a><a href="zlacp2.f.html#ZLACP2.1">ZLACP2</a>( <span class="string">'F'</span>, N, N, RWORK( IRU ), N, WORK( IU ),
     $                      LDWRKU )
               CALL <a name="ZUNMBR.618"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, N, N, N, WORK( IR ), LDWRKR,
     $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy real matrix RWORK(IRVT) to complex matrix VT
</span><span class="comment">*</span><span class="comment">              Overwrite VT by the right singular vectors of R
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACP2.627"></a><a href="zlacp2.f.html#ZLACP2.1">ZLACP2</a>( <span class="string">'F'</span>, N, N, RWORK( IRVT ), N, VT, LDVT )
               CALL <a name="ZUNMBR.628"></a><a href="zunmbr.f.html#ZUNMBR.1">ZUNMBR</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, <span class="string">'C'</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(IU), storing result in WORK(IR) and copying to A
</span><span class="comment">*</span><span class="comment">              (CWorkspace: need 2*N*N, prefer N*N+M*N)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               DO 10 I = 1, M, LDWRKR
                  CHUNK = MIN( M-I+1, LDWRKR )
                  CALL ZGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, CHUNK, N, N, CONE, A( I, 1 ),
     $                        LDA, WORK( IU ), LDWRKU, CZERO,
     $                        WORK( IR ), LDWRKR )
                  CALL <a name="ZLACPY.642"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</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">              (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZGEQRF.664"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</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="ZLACPY.669"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="ZLASET.670"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'L'</span>, N-1, N-1, CZERO, CZERO, 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">              (CWorkspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">              (RWorkspace: 0)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZUNGQR.677"></a><a href="zungqr.f.html#ZUNGQR.1">ZUNGQR</a>( M, N, N, A, LDA, WORK( ITAU ),
     $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
               IE = 1
               ITAUQ = ITAU
               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">              (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*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="ZGEBRD.688"></a><a href="zgebrd.f.html#ZGEBRD.1">ZGEBRD</a>( N, N, WORK( IR ), LDWRKR, S, RWORK( 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 RWORK(IRU) and computing right
</span><span class="comment">*</span><span class="
12345

zgesdd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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