sgesdd.f.html

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

HTML
742 行
字号
12345
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.431"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGESDD.431"></a><a href="sgesdd.f.html#SGESDD.1">SGESDD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine constants
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.445"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = SQRT( <a name="SLAMCH.446"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) ) / EPS
      BIGNUM = ONE / SMLNUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale A if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="SLANGE.451"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, DUM )
      ISCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
         ISCL = 1
         CALL <a name="SLASCL.455"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
      ELSE IF( ANRM.GT.BIGNUM ) THEN
         ISCL = 1
         CALL <a name="SLASCL.458"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( M.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        A has at least as many rows as columns. If A has sufficiently
</span><span class="comment">*</span><span class="comment">        more rows than columns, first reduce using the QR
</span><span class="comment">*</span><span class="comment">        decomposition (if sufficient workspace available)
</span><span class="comment">*</span><span class="comment">
</span>         IF( M.GE.MNTHR ) THEN
<span class="comment">*</span><span class="comment">
</span>            IF( WNTQN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 1 (M much larger than N, JOBZ='N')
</span><span class="comment">*</span><span class="comment">              No singular vectors to be computed
</span><span class="comment">*</span><span class="comment">
</span>               ITAU = 1
               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 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SGEQRF.480"></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">              Zero out below R
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASET.485"></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 = 1
               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 4*N, prefer 3*N+2*N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SGEBRD.494"></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 )
               NWORK = IE + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Perform bidiagonal SVD, computing singular values only
</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.502"></a><a href="sbdsdc.f.html#SBDSDC.1">SBDSDC</a>( <span class="string">'U'</span>, <span class="string">'N'</span>, N, S, WORK( IE ), DUM, 1, DUM, 1,
     $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span>            ELSE IF( WNTQO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 2 (M much larger than N, JOBZ = 'O')
</span><span class="comment">*</span><span class="comment">              N left singular vectors to be overwritten on A 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 LDWRKR by N
</span><span class="comment">*</span><span class="comment">
</span>               IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
                  LDWRKR = LDA
               ELSE
                  LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / 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">              (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.526"></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.531"></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.532"></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.538"></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 VT, copying result to 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.548"></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">              WORK(IU) is N by N
</span><span class="comment">*</span><span class="comment">
</span>               IU = NWORK
               NWORK = IU + N*N
<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.562"></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
</span><span class="comment">*</span><span class="comment">              and VT by right singular vectors of R
</span><span class="comment">*</span><span class="comment">              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SORMBR.570"></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 ), WORK( IU ), N, WORK( NWORK ),
     $                      LWORK-NWORK+1, IERR )
               CALL <a name="SORMBR.573"></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
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 + -
显示快捷键?