slaqr2.f.html

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

HTML
576 行
字号
123
</span>      REAL               <a name="SLAMCH.174"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.175"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SGEHRD.178"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>, SGEMM, <a name="SLABAD.178"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>, <a name="SLACPY.178"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>, <a name="SLAHQR.178"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>,
     $                   <a name="SLANV2.179"></a><a href="slanv2.f.html#SLANV2.1">SLANV2</a>, <a name="SLARF.179"></a><a href="slarf.f.html#SLARF.1">SLARF</a>, <a name="SLARFG.179"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>, <a name="SLASET.179"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, <a name="SORGHR.179"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>, <a name="STREXC.179"></a><a href="strexc.f.html#STREXC.1">STREXC</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, INT, MAX, MIN, REAL, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Estimate optimal workspace. ====
</span><span class="comment">*</span><span class="comment">
</span>      JW = MIN( NW, KBOT-KTOP+1 )
      IF( JW.LE.2 ) THEN
         LWKOPT = 1
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Workspace query call to <a name="SGEHRD.193"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEHRD.195"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
         LWK1 = INT( WORK( 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Workspace query call to <a name="SORGHR.198"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORGHR.200"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
         LWK2 = INT( WORK( 1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Optimal workspace ====
</span><span class="comment">*</span><span class="comment">
</span>         LWKOPT = JW + MAX( LWK1, LWK2 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Quick return in case of workspace query. ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( LWORK.EQ.-1 ) THEN
         WORK( 1 ) = REAL( LWKOPT )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Nothing to do ...
</span><span class="comment">*</span><span class="comment">     ... for an empty active block ... ====
</span>      NS = 0
      ND = 0
      IF( KTOP.GT.KBOT )
     $   RETURN
<span class="comment">*</span><span class="comment">     ... nor for an empty deflation window. ====
</span>      IF( NW.LT.1 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Machine constants ====
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="SLAMCH.227"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = ONE / SAFMIN
      CALL <a name="SLABAD.229"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="SLAMCH.230"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( REAL( N ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Setup deflation window ====
</span><span class="comment">*</span><span class="comment">
</span>      JW = MIN( NW, KBOT-KTOP+1 )
      KWTOP = KBOT - JW + 1
      IF( KWTOP.EQ.KTOP ) THEN
         S = ZERO
      ELSE
         S = H( KWTOP, KWTOP-1 )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( KBOT.EQ.KWTOP ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== 1-by-1 deflation window: not much to do ====
</span><span class="comment">*</span><span class="comment">
</span>         SR( KWTOP ) = H( KWTOP, KWTOP )
         SI( KWTOP ) = ZERO
         NS = 1
         ND = 0
         IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) )
     $        THEN
            NS = 0
            ND = 1
            IF( KWTOP.GT.KTOP )
     $         H( KWTOP, KWTOP-1 ) = ZERO
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Convert to spike-triangular form.  (In case of a
</span><span class="comment">*</span><span class="comment">     .    rare QR failure, this routine continues to do
</span><span class="comment">*</span><span class="comment">     .    aggressive early deflation using that part of
</span><span class="comment">*</span><span class="comment">     .    the deflation window that converged using INFQR
</span><span class="comment">*</span><span class="comment">     .    here and there to keep track.) ====
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLACPY.267"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'U'</span>, JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
      CALL SCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASET.270"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, JW, JW, ZERO, ONE, V, LDV )
      CALL <a name="SLAHQR.271"></a><a href="slahqr.f.html#SLAHQR.1">SLAHQR</a>( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
     $             SI( KWTOP ), 1, JW, V, LDV, INFQR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== <a name="STREXC.274"></a><a href="strexc.f.html#STREXC.1">STREXC</a> needs a clean margin near the diagonal ====
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 J = 1, JW - 3
         T( J+2, J ) = ZERO
         T( J+3, J ) = ZERO
   10 CONTINUE
      IF( JW.GT.2 )
     $   T( JW, JW-2 ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Deflation detection loop ====
</span><span class="comment">*</span><span class="comment">
</span>      NS = JW
      ILST = INFQR + 1
   20 CONTINUE
      IF( ILST.LE.NS ) THEN
         IF( NS.EQ.1 ) THEN
            BULGE = .FALSE.
         ELSE
            BULGE = T( NS, NS-1 ).NE.ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Small spike tip test for deflation ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( .NOT.BULGE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Real eigenvalue ====
</span><span class="comment">*</span><span class="comment">
</span>            FOO = ABS( T( NS, NS ) )
            IF( FOO.EQ.ZERO )
     $         FOO = ABS( S )
            IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Deflatable ====
</span><span class="comment">*</span><span class="comment">
</span>               NS = NS - 1
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Undeflatable.   Move it up out of the way.
</span><span class="comment">*</span><span class="comment">              .    (<a name="STREXC.312"></a><a href="strexc.f.html#STREXC.1">STREXC</a> can not fail in this case.) ====
</span><span class="comment">*</span><span class="comment">
</span>               IFST = NS
               CALL <a name="STREXC.315"></a><a href="strexc.f.html#STREXC.1">STREXC</a>( <span class="string">'V'</span>, JW, T, LDT, V, LDV, IFST, ILST, WORK,
     $                      INFO )
               ILST = ILST + 1
            END IF
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Complex conjugate pair ====
</span><span class="comment">*</span><span class="comment">
</span>            FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
     $            SQRT( ABS( T( NS-1, NS ) ) )
            IF( FOO.EQ.ZERO )
     $         FOO = ABS( S )
            IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE.
     $          MAX( SMLNUM, ULP*FOO ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Deflatable ====
</span><span class="comment">*</span><span class="comment">
</span>               NS = NS - 2
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Undflatable. Move them up out of the way.
</span><span class="comment">*</span><span class="comment">              .    Fortunately, <a name="STREXC.336"></a><a href="strexc.f.html#STREXC.1">STREXC</a> does the right thing with
</span><span class="comment">*</span><span class="comment">              .    ILST in case of a rare exchange failure. ====
</span><span class="comment">*</span><span class="comment">
</span>               IFST = NS
               CALL <a name="STREXC.340"></a><a href="strexc.f.html#STREXC.1">STREXC</a>( <span class="string">'V'</span>, JW, T, LDT, V, LDV, IFST, ILST, WORK,
     $                      INFO )
               ILST = ILST + 2
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== End deflation detection loop ====
</span><span class="comment">*</span><span class="comment">
</span>         GO TO 20
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Return to Hessenberg form ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( NS.EQ.0 )
     $   S = ZERO
<span class="comment">*</span><span class="comment">
</span>      IF( NS.LT.JW ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== sorting diagonal blocks of T improves accuracy for
</span><span class="comment">*</span><span class="comment">        .    graded matrices.  Bubble sort deals well with
</span><span class="comment">*</span><span class="comment">        .    exchange failures. ====
</span><span class="comment">*</span><span class="comment">
</span>         SORTED = .false.
         I = NS + 1
   30    CONTINUE
         IF( SORTED )
123

slaqr2.f.html - 源码说明

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

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

⌨️ 快捷键说明

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