zlalsd.f.html

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

HTML
625 行
字号
123
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
               NSUB = NSUB + 1
               IWORK( NSUB ) = N
               IWORK( SIZEI+NSUB-1 ) = 1
               CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
            END IF
            ST1 = ST - 1
            IF( NSIZE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a 1-by-1 subproblem and is not solved
</span><span class="comment">*</span><span class="comment">              explicitly.
</span><span class="comment">*</span><span class="comment">
</span>               CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
            ELSE IF( NSIZE.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a small subproblem and is solved by <a name="DLASDQ.417"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLASET.419"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, NSIZE, NSIZE, ZERO, ONE,
     $                      RWORK( VT+ST1 ), N )
               CALL <a name="DLASET.421"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, NSIZE, NSIZE, ZERO, ONE,
     $                      RWORK( U+ST1 ), N )
               CALL <a name="DLASDQ.423"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a>( <span class="string">'U'</span>, 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
     $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
     $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
     $                      INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              In the real version, B is passed to <a name="DLASDQ.431"></a><a href="dlasdq.f.html#DLASDQ.1">DLASDQ</a> and multiplied
</span><span class="comment">*</span><span class="comment">              internally by Q'. Here B is complex and that product is
</span><span class="comment">*</span><span class="comment">              computed below in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span>               J = IRWB - 1
               DO 190 JCOL = 1, NRHS
                  DO 180 JROW = ST, ST + NSIZE - 1
                     J = J + 1
                     RWORK( J ) = DBLE( B( JROW, JCOL ) )
  180             CONTINUE
  190          CONTINUE
               CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
     $                     ZERO, RWORK( IRWRB ), NSIZE )
               J = IRWB - 1
               DO 210 JCOL = 1, NRHS
                  DO 200 JROW = ST, ST + NSIZE - 1
                     J = J + 1
                     RWORK( J ) = DIMAG( B( JROW, JCOL ) )
  200             CONTINUE
  210          CONTINUE
               CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
     $                     ZERO, RWORK( IRWIB ), NSIZE )
               JREAL = IRWRB - 1
               JIMAG = IRWIB - 1
               DO 230 JCOL = 1, NRHS
                  DO 220 JROW = ST, ST + NSIZE - 1
                     JREAL = JREAL + 1
                     JIMAG = JIMAG + 1
                     B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                                 RWORK( JIMAG ) )
  220             CONTINUE
  230          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACPY.466"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'A'</span>, NSIZE, NRHS, B( ST, 1 ), LDB,
     $                      WORK( BX+ST1 ), N )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A large problem. Solve it using divide and conquer.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLASDA.472"></a><a href="dlasda.f.html#DLASDA.1">DLASDA</a>( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
     $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
     $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
     $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
     $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
     $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
     $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
     $                      RWORK( S+ST1 ), RWORK( NRWORK ),
     $                      IWORK( IWK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
               BXST = BX + ST1
               CALL <a name="ZLALSA.485"></a><a href="zlalsa.f.html#ZLALSA.1">ZLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
     $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
     $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
     $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
     $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
     $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
     $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
     $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            ST = I + 1
         END IF
  240 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply the singular values and treat the tiny ones as zero.
</span><span class="comment">*</span><span class="comment">
</span>      TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      DO 250 I = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Some of the elements in D can be negative because 1-by-1
</span><span class="comment">*</span><span class="comment">        subproblems were not solved explicitly.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( D( I ) ).LE.TOL ) THEN
            CALL <a name="ZLASET.512"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
         ELSE
            RANK = RANK + 1
            CALL <a name="ZLASCL.515"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS,
     $                   WORK( BX+I-1 ), N, INFO )
         END IF
         D( I ) = ABS( D( I ) )
  250 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now apply back the right singular vectors.
</span><span class="comment">*</span><span class="comment">
</span>      ICMPQ2 = 1
      DO 320 I = 1, NSUB
         ST = IWORK( I )
         ST1 = ST - 1
         NSIZE = IWORK( SIZEI+I-1 )
         BXST = BX + ST1
         IF( NSIZE.EQ.1 ) THEN
            CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
         ELSE IF( NSIZE.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Since B and BX are complex, the following call to DGEMM
</span><span class="comment">*</span><span class="comment">           is performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
</span><span class="comment">*</span><span class="comment">    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
</span><span class="comment">*</span><span class="comment">    $                  B( ST, 1 ), LDB )
</span><span class="comment">*</span><span class="comment">
</span>            J = BXST - N - 1
            JREAL = IRWB - 1
            DO 270 JCOL = 1, NRHS
               J = J + N
               DO 260 JROW = 1, NSIZE
                  JREAL = JREAL + 1
                  RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
  260          CONTINUE
  270       CONTINUE
            CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
     $                  RWORK( IRWRB ), NSIZE )
            J = BXST - N - 1
            JIMAG = IRWB - 1
            DO 290 JCOL = 1, NRHS
               J = J + N
               DO 280 JROW = 1, NSIZE
                  JIMAG = JIMAG + 1
                  RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
  280          CONTINUE
  290       CONTINUE
            CALL DGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
     $                  RWORK( IRWIB ), NSIZE )
            JREAL = IRWRB - 1
            JIMAG = IRWIB - 1
            DO 310 JCOL = 1, NRHS
               DO 300 JROW = ST, ST + NSIZE - 1
                  JREAL = JREAL + 1
                  JIMAG = JIMAG + 1
                  B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
     $                              RWORK( JIMAG ) )
  300          CONTINUE
  310       CONTINUE
         ELSE
            CALL <a name="ZLALSA.575"></a><a href="zlalsa.f.html#ZLALSA.1">ZLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
     $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
     $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
     $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
     $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
     $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
     $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
     $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
            IF( INFO.NE.0 ) THEN
               RETURN
            END IF
         END IF
  320 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Unscale and sort the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASCL.592"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
      CALL <a name="DLASRT.593"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
      CALL <a name="ZLASCL.594"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZLALSD.598"></a><a href="zlalsd.f.html#ZLALSD.1">ZLALSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
123

zlalsd.f.html - 源码说明

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

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

⌨️ 快捷键说明

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