slaqr2.f.html

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

     $      GO TO 50
         SORTED = .true.
<span class="comment">*</span><span class="comment">
</span>         KEND = I - 1
         I = INFQR + 1
         IF( I.EQ.NS ) THEN
            K = I + 1
         ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
            K = I + 1
         ELSE
            K = I + 2
         END IF
   40    CONTINUE
         IF( K.LE.KEND ) THEN
            IF( K.EQ.I+1 ) THEN
               EVI = ABS( T( I, I ) )
            ELSE
               EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )*
     $               SQRT( ABS( T( I, I+1 ) ) )
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( K.EQ.KEND ) THEN
               EVK = ABS( T( K, K ) )
            ELSE IF( T( K+1, K ).EQ.ZERO ) THEN
               EVK = ABS( T( K, K ) )
            ELSE
               EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )*
     $               SQRT( ABS( T( K, K+1 ) ) )
            END IF
<span class="comment">*</span><span class="comment">
</span>            IF( EVI.GE.EVK ) THEN
               I = K
            ELSE
               SORTED = .false.
               IFST = I
               ILST = K
               CALL <a name="STREXC.402"></a><a href="strexc.f.html#STREXC.1">STREXC</a>( <span class="string">'V'</span>, JW, T, LDT, V, LDV, IFST, ILST, WORK,
     $                      INFO )
               IF( INFO.EQ.0 ) THEN
                  I = ILST
               ELSE
                  I = K
               END IF
            END IF
            IF( I.EQ.KEND ) THEN
               K = I + 1
            ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
               K = I + 1
            ELSE
               K = I + 2
            END IF
            GO TO 40
         END IF
         GO TO 30
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Restore shift/eigenvalue array from T ====
</span><span class="comment">*</span><span class="comment">
</span>      I = JW
   60 CONTINUE
      IF( I.GE.INFQR+1 ) THEN
         IF( I.EQ.INFQR+1 ) THEN
            SR( KWTOP+I-1 ) = T( I, I )
            SI( KWTOP+I-1 ) = ZERO
            I = I - 1
         ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN
            SR( KWTOP+I-1 ) = T( I, I )
            SI( KWTOP+I-1 ) = ZERO
            I = I - 1
         ELSE
            AA = T( I-1, I-1 )
            CC = T( I, I-1 )
            BB = T( I-1, I )
            DD = T( I, I )
            CALL <a name="SLANV2.441"></a><a href="slanv2.f.html#SLANV2.1">SLANV2</a>( AA, BB, CC, DD, SR( KWTOP+I-2 ),
     $                   SI( KWTOP+I-2 ), SR( KWTOP+I-1 ),
     $                   SI( KWTOP+I-1 ), CS, SN )
            I = I - 2
         END IF
         GO TO 60
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
         IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Reflect spike back into lower triangle ====
</span><span class="comment">*</span><span class="comment">
</span>            CALL SCOPY( NS, V, LDV, WORK, 1 )
            BETA = WORK( 1 )
            CALL <a name="SLARFG.456"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( NS, BETA, WORK( 2 ), 1, TAU )
            WORK( 1 ) = ONE
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLASET.459"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'L'</span>, JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLARF.461"></a><a href="slarf.f.html#SLARF.1">SLARF</a>( <span class="string">'L'</span>, NS, JW, WORK, 1, TAU, T, LDT,
     $                  WORK( JW+1 ) )
            CALL <a name="SLARF.463"></a><a href="slarf.f.html#SLARF.1">SLARF</a>( <span class="string">'R'</span>, NS, NS, WORK, 1, TAU, T, LDT,
     $                  WORK( JW+1 ) )
            CALL <a name="SLARF.465"></a><a href="slarf.f.html#SLARF.1">SLARF</a>( <span class="string">'R'</span>, JW, NS, WORK, 1, TAU, V, LDV,
     $                  WORK( JW+1 ) )
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="SGEHRD.468"></a><a href="sgehrd.f.html#SGEHRD.1">SGEHRD</a>( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
     $                   LWORK-JW, INFO )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Copy updated reduced window into place ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( KWTOP.GT.1 )
     $      H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
         CALL <a name="SLACPY.476"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'U'</span>, JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
         CALL SCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
     $               LDH+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Accumulate orthogonal matrix in order update
</span><span class="comment">*</span><span class="comment">        .    H and Z, if requested.  (A modified version
</span><span class="comment">*</span><span class="comment">        .    of  <a name="SORGHR.482"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a> that accumulates block Householder
</span><span class="comment">*</span><span class="comment">        .    transformations into V directly might be
</span><span class="comment">*</span><span class="comment">        .    marginally more efficient than the following.) ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
            CALL <a name="SORGHR.487"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
     $                   LWORK-JW, INFO )
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JW, NS, NS, ONE, V, LDV, T, LDT, ZERO,
     $                  WV, LDWV )
            CALL <a name="SLACPY.491"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, JW, NS, WV, LDWV, V, LDV )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update vertical slab in H ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTT ) THEN
            LTOP = 1
         ELSE
            LTOP = KTOP
         END IF
         DO 70 KROW = LTOP, KWTOP - 1, NV
            KLN = MIN( NV, KWTOP-KROW )
            CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, KLN, JW, JW, ONE, H( KROW, KWTOP ),
     $                  LDH, V, LDV, ZERO, WV, LDWV )
            CALL <a name="SLACPY.505"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
   70    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update horizontal slab in H ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTT ) THEN
            DO 80 KCOL = KBOT + 1, N, NH
               KLN = MIN( NH, N-KCOL+1 )
               CALL SGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, JW, KLN, JW, ONE, V, LDV,
     $                     H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
               CALL <a name="SLACPY.515"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, JW, KLN, T, LDT, H( KWTOP, KCOL ),
     $                      LDH )
   80       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Update vertical slab in Z ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( WANTZ ) THEN
            DO 90 KROW = ILOZ, IHIZ, NV
               KLN = MIN( NV, IHIZ-KROW+1 )
               CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, KLN, JW, JW, ONE, Z( KROW, KWTOP ),
     $                     LDZ, V, LDV, ZERO, WV, LDWV )
               CALL <a name="SLACPY.527"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'A'</span>, KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
     $                      LDZ )
   90       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Return the number of deflations ... ====
</span><span class="comment">*</span><span class="comment">
</span>      ND = JW - NS
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== ... and the number of shifts. (Subtracting
</span><span class="comment">*</span><span class="comment">     .    INFQR from the spike length takes care
</span><span class="comment">*</span><span class="comment">     .    of the case of a rare QR failure while
</span><span class="comment">*</span><span class="comment">     .    calculating eigenvalues of the deflation
</span><span class="comment">*</span><span class="comment">     .    window.)  ====
</span><span class="comment">*</span><span class="comment">
</span>      NS = NS - INFQR
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      ==== Return optimal workspace. ====
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = REAL( LWKOPT )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="SLAQR2.549"></a><a href="slaqr2.f.html#SLAQR2.1">SLAQR2</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>