cgesdd.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="CGEQRF.570"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</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="CLACPY.575"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="CLASET.576"></a><a href="claset.f.html#CLASET.1">CLASET</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="CUNGQR.583"></a><a href="cungqr.f.html#CUNGQR.1">CUNGQR</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="CGEBRD.594"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</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="SBDSDC.607"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</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="CLACP2.616"></a><a href="clacp2.f.html#CLACP2.1">CLACP2</a>( <span class="string">'F'</span>, N, N, RWORK( IRU ), N, WORK( IU ),
     $                      LDWRKU )
               CALL <a name="CUNMBR.618"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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="CLACP2.627"></a><a href="clacp2.f.html#CLACP2.1">CLACP2</a>( <span class="string">'F'</span>, N, N, RWORK( IRVT ), N, VT, LDVT )
               CALL <a name="CUNMBR.628"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</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 CGEMM( <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="CLACPY.642"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</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="CGEQRF.664"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</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="CLACPY.669"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'U'</span>, N, N, A, LDA, WORK( IR ), LDWRKR )
               CALL <a name="CLASET.670"></a><a href="claset.f.html#CLASET.1">CLASET</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="CUNGQR.677"></a><a href="cungqr.f.html#CUNGQR.1">CUNGQR</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="CGEBRD.688"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</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="comment">              singular vectors of b
12345

cgesdd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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