slaqr5.f.html

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

</span>            K = KRCOL + 3*( M22-1 )
            IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
               DO 100 J = JTOP, MIN( KBOT, K+3 )
                  REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
     $                     H( J, K+2 ) )
                  H( J, K+1 ) = H( J, K+1 ) - REFSUM
                  H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M22 )
  100          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ACCUM ) THEN
                  KMS = K - INCOL
                  DO 110 J = MAX( 1, KTOP-INCOL ), KDU
                     REFSUM = V( 1, M22 )*( U( J, KMS+1 )+V( 2, M22 )*
     $                        U( J, KMS+2 ) )
                     U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
                     U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M22 )
  110             CONTINUE
               ELSE IF( WANTZ ) THEN
                  DO 120 J = ILOZ, IHIZ
                     REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
     $                        Z( J, K+2 ) )
                     Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
                     Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M22 )
  120             CONTINUE
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Vigilant deflation check ====
</span><span class="comment">*</span><span class="comment">
</span>            MSTART = MTOP
            IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
     $         MSTART = MSTART + 1
            MEND = MBOT
            IF( BMP22 )
     $         MEND = MEND + 1
            IF( KRCOL.EQ.KBOT-2 )
     $         MEND = MEND + 1
            DO 130 M = MSTART, MEND
               K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== The following convergence test requires that
</span><span class="comment">*</span><span class="comment">              .    the tradition small-compared-to-nearby-diagonals
</span><span class="comment">*</span><span class="comment">              .    criterion and the Ahues &amp; Tisseur (LAWN 122, 1997)
</span><span class="comment">*</span><span class="comment">              .    criteria both be satisfied.  The latter improves
</span><span class="comment">*</span><span class="comment">              .    accuracy in some examples. Falling back on an
</span><span class="comment">*</span><span class="comment">              .    alternate convergence criterion when TST1 or TST2
</span><span class="comment">*</span><span class="comment">              .    is zero (as done here) is traditional but probably 
</span><span class="comment">*</span><span class="comment">              .    unnecessary. ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( H( K+1, K ).NE.ZERO ) THEN
                  TST1 = ABS( H( K, K ) ) + ABS( H( K+1, K+1 ) )
                  IF( TST1.EQ.ZERO ) THEN
                     IF( K.GE.KTOP+1 )
     $                  TST1 = TST1 + ABS( H( K, K-1 ) )
                     IF( K.GE.KTOP+2 )
     $                  TST1 = TST1 + ABS( H( K, K-2 ) )
                     IF( K.GE.KTOP+3 )
     $                  TST1 = TST1 + ABS( H( K, K-3 ) )
                     IF( K.LE.KBOT-2 )
     $                  TST1 = TST1 + ABS( H( K+2, K+1 ) )
                     IF( K.LE.KBOT-3 )
     $                  TST1 = TST1 + ABS( H( K+3, K+1 ) )
                     IF( K.LE.KBOT-4 )
     $                  TST1 = TST1 + ABS( H( K+4, K+1 ) )
                  END IF
                  IF( ABS( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
     $                 THEN
                     H12 = MAX( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )
                     H21 = MIN( ABS( H( K+1, K ) ), ABS( H( K, K+1 ) ) )
                     H11 = MAX( ABS( H( K+1, K+1 ) ),
     $                     ABS( H( K, K )-H( K+1, K+1 ) ) )
                     H22 = MIN( ABS( H( K+1, K+1 ) ),
     $                     ABS( H( K, K )-H( K+1, K+1 ) ) )
                     SCL = H11 + H12
                     TST2 = H22*( H11 / SCL )
<span class="comment">*</span><span class="comment">
</span>                     IF( TST2.EQ.ZERO .OR. H21*( H12 / SCL ).LE.
     $                   MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
                  END IF
               END IF
  130       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Fill in the last row of each bulge. ====
</span><span class="comment">*</span><span class="comment">
</span>            MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
            DO 140 M = MTOP, MEND
               K = KRCOL + 3*( M-1 )
               REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
               H( K+4, K+1 ) = -REFSUM
               H( K+4, K+2 ) = -REFSUM*V( 2, M )
               H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*V( 3, M )
  140       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== End of near-the-diagonal bulge chase. ====
</span><span class="comment">*</span><span class="comment">
</span>  150    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Use U (if accumulated) to update far-from-diagonal
</span><span class="comment">*</span><span class="comment">        .    entries in H.  If required, use U to update Z as
</span><span class="comment">*</span><span class="comment">        .    well. ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( ACCUM ) THEN
            IF( WANTT ) THEN
               JTOP = 1
               JBOT = N
            ELSE
               JTOP = KTOP
               JBOT = KBOT
            END IF
            IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
     $          ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Updates not exploiting the 2-by-2 block
</span><span class="comment">*</span><span class="comment">              .    structure of U.  K1 and NU keep track of
</span><span class="comment">*</span><span class="comment">              .    the location and size of U in the special
</span><span class="comment">*</span><span class="comment">              .    cases of introducing bulges and chasing
</span><span class="comment">*</span><span class="comment">              .    bulges off the bottom.  In these special
</span><span class="comment">*</span><span class="comment">              .    cases and in case the number of shifts
</span><span class="comment">*</span><span class="comment">              .    is NS = 2, there is no 2-by-2 block
</span><span class="comment">*</span><span class="comment">              .    structure to exploit.  ====
</span><span class="comment">*</span><span class="comment">
</span>               K1 = MAX( 1, KTOP-INCOL )
               NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Horizontal Multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 160 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
                  JLEN = MIN( NH, JBOT-JCOL+1 )
                  CALL SGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, NU, JLEN, NU, ONE, U( K1, K1 ),
     $                        LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
     $                        LDWH )
                  CALL <a name="SLACPY.617"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'ALL'</span>, NU, JLEN, WH, LDWH,
     $                         H( INCOL+K1, JCOL ), LDH )
  160          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Vertical multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 170 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
                  JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
                  CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, NU, NU, ONE,
     $                        H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
     $                        LDU, ZERO, WV, LDWV )
                  CALL <a name="SLACPY.628"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'ALL'</span>, JLEN, NU, WV, LDWV,
     $                         H( JROW, INCOL+K1 ), LDH )
  170          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Z multiply (also vertical) ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( WANTZ ) THEN
                  DO 180 JROW = ILOZ, IHIZ, NV
                     JLEN = MIN( NV, IHIZ-JROW+1 )
                     CALL SGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, NU, NU, ONE,
     $                           Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
     $                           LDU, ZERO, WV, LDWV )
                     CALL <a name="SLACPY.640"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'ALL'</span>, JLEN, NU, WV, LDWV,
     $                            Z( JROW, INCOL+K1 ), LDZ )
  180             CONTINUE
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Updates exploiting U's 2-by-2 block structure.
</span><span class="comment">*</span><span class="comment">              .    (I2, I4, J2, J4 are the last rows and columns
</span><span class="comment">*</span><span class="comment">              .    of the blocks.) ====
</span><span class="comment">*</span><span class="comment">
</span>               I2 = ( KDU+1 ) / 2
               I4 = KDU
               J2 = I4 - I2
               J4 = KDU