dlaqr3.f.html

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

HTML
586 行
字号
123
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, INT, MAX, MIN, 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="DGEHRD.191"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DGEHRD.193"></a><a href="dgehrd.f.html#DGEHRD.1">DGEHRD</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="DORGHR.196"></a><a href="dorghr.f.html#DORGHR.1">DORGHR</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DORGHR.198"></a><a href="dorghr.f.html#DORGHR.1">DORGHR</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">        ==== Workspace query call to <a name="DLAQR4.201"></a><a href="dlaqr4.f.html#DLAQR4.1">DLAQR4</a> ====
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLAQR4.203"></a><a href="dlaqr4.f.html#DLAQR4.1">DLAQR4</a>( .true., .true., JW, 1, JW, T, LDT, SR, SI, 1, JW,
     $                V, LDV, WORK, -1, INFQR )
         LWK3 = 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 = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
      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 ) = DBLE( 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="DLAMCH.231"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = ONE / SAFMIN
      CALL <a name="DLABAD.233"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="DLAMCH.234"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( DBLE( 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="DLACPY.271"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'U'</span>, JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
      CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLASET.274"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'A'</span>, JW, JW, ZERO, ONE, V, LDV )
      NMIN = <a name="ILAENV.275"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 12, <span class="string">'<a name="DLAQR3.275"></a><a href="dlaqr3.f.html#DLAQR3.1">DLAQR3</a>'</span>, <span class="string">'SV'</span>, JW, 1, JW, LWORK )
      IF( JW.GT.NMIN ) THEN
         CALL <a name="DLAQR4.277"></a><a href="dlaqr4.f.html#DLAQR4.1">DLAQR4</a>( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
     $                SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR )
      ELSE
         CALL <a name="DLAHQR.280"></a><a href="dlahqr.f.html#DLAHQR.1">DLAHQR</a>( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
     $                SI( KWTOP ), 1, JW, V, LDV, INFQR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== <a name="DTREXC.284"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</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="DTREXC.322"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</a> can not fail in this case.) ====
</span><span class="comment">*</span><span class="comment">
</span>               IFST = NS
               CALL <a name="DTREXC.325"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</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="DTREXC.346"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</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="DTREXC.350"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</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
123

dlaqr3.f.html - 源码说明

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

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

⌨️ 快捷键说明

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