slaqr5.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="SLAQR1.319"></a><a href="slaqr1.f.html#SLAQR1.1">SLAQR1</a>( 3, H( K+1, K+1 ), LDH, SR( 2*M-1 ),
     $                            SI( 2*M-1 ), SR( 2*M ), SI( 2*M ),
     $                            VT )
                     SCL = ABS( VT( 1 ) ) + ABS( VT( 2 ) ) +
     $                     ABS( VT( 3 ) )
                     IF( SCL.NE.ZERO ) THEN
                        VT( 1 ) = VT( 1 ) / SCL
                        VT( 2 ) = VT( 2 ) / SCL
                        VT( 3 ) = VT( 3 ) / SCL
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== The following is the traditional and
</span><span class="comment">*</span><span class="comment">                    .    conservative two-small-subdiagonals
</span><span class="comment">*</span><span class="comment">                    .    test.  ====
</span><span class="comment">*</span><span class="comment">                    .
</span>                     IF( ABS( H( K+1, K ) )*( ABS( VT( 2 ) )+
     $                   ABS( VT( 3 ) ) ).GT.ULP*ABS( VT( 1 ) )*
     $                   ( ABS( H( K, K ) )+ABS( H( K+1,
     $                   K+1 ) )+ABS( H( K+2, K+2 ) ) ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       ==== Starting a new bulge here would
</span><span class="comment">*</span><span class="comment">                       .    create non-negligible fill.   If
</span><span class="comment">*</span><span class="comment">                       .    the old reflector is diagonal (only
</span><span class="comment">*</span><span class="comment">                       .    possible with underflows), then
</span><span class="comment">*</span><span class="comment">                       .    change it to I.  Otherwise, use
</span><span class="comment">*</span><span class="comment">                       .    it with trepidation. ====
</span><span class="comment">*</span><span class="comment">
</span>                        IF( V( 2, M ).EQ.ZERO .AND. V( 3, M ).EQ.ZERO )
     $                       THEN
                           V( 1, M ) = ZERO
                        ELSE
                           H( K+1, K ) = BETA
                           H( K+2, K ) = ZERO
                           H( K+3, K ) = ZERO
                        END IF
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       ==== Stating a new bulge here would
</span><span class="comment">*</span><span class="comment">                       .    create only negligible fill.
</span><span class="comment">*</span><span class="comment">                       .    Replace the old reflector with
</span><span class="comment">*</span><span class="comment">                       .    the new one. ====
</span><span class="comment">*</span><span class="comment">
</span>                        ALPHA = VT( 1 )
                        CALL <a name="SLARFG.362"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
                        REFSUM = H( K+1, K ) + H( K+2, K )*VT( 2 ) +
     $                           H( K+3, K )*VT( 3 )
                        H( K+1, K ) = H( K+1, K ) - VT( 1 )*REFSUM
                        H( K+2, K ) = ZERO
                        H( K+3, K ) = ZERO
                        V( 1, M ) = VT( 1 )
                        V( 2, M ) = VT( 2 )
                        V( 3, M ) = VT( 3 )
                     END IF
                  END IF
               END IF
   20       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Generate a 2-by-2 reflection, if needed. ====
</span><span class="comment">*</span><span class="comment">
</span>            K = KRCOL + 3*( M22-1 )
            IF( BMP22 ) THEN
               IF( K.EQ.KTOP-1 ) THEN
                  CALL <a name="SLAQR1.381"></a><a href="slaqr1.f.html#SLAQR1.1">SLAQR1</a>( 2, H( K+1, K+1 ), LDH, SR( 2*M22-1 ),
     $                         SI( 2*M22-1 ), SR( 2*M22 ), SI( 2*M22 ),
     $                         V( 1, M22 ) )
                  BETA = V( 1, M22 )
                  CALL <a name="SLARFG.385"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
               ELSE
                  BETA = H( K+1, K )
                  V( 2, M22 ) = H( K+2, K )
                  CALL <a name="SLARFG.389"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
                  H( K+1, K ) = BETA
                  H( K+2, K ) = ZERO
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Initialize V(1,M22) here to avoid possible undefined
</span><span class="comment">*</span><span class="comment">              .    variable problems later. ====
</span><span class="comment">*</span><span class="comment">
</span>               V( 1, M22 ) = ZERO
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Multiply H by reflections from the left ====
</span><span class="comment">*</span><span class="comment">
</span>            IF( ACCUM ) THEN
               JBOT = MIN( NDCOL, KBOT )
            ELSE IF( WANTT ) THEN
               JBOT = N
            ELSE
               JBOT = KBOT
            END IF
            DO 40 J = MAX( KTOP, KRCOL ), JBOT
               MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
               DO 30 M = MTOP, MEND
                  K = KRCOL + 3*( M-1 )
                  REFSUM = V( 1, M )*( H( K+1, J )+V( 2, M )*
     $                     H( K+2, J )+V( 3, M )*H( K+3, J ) )
                  H( K+1, J ) = H( K+1, J ) - REFSUM
                  H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
                  H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
   30          CONTINUE
   40       CONTINUE
            IF( BMP22 ) THEN
               K = KRCOL + 3*( M22-1 )
               DO 50 J = MAX( K+1, KTOP ), JBOT
                  REFSUM = V( 1, M22 )*( H( K+1, J )+V( 2, M22 )*
     $                     H( K+2, J ) )
                  H( K+1, J ) = H( K+1, J ) - REFSUM
                  H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
   50          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Multiply H by reflections from the right.
</span><span class="comment">*</span><span class="comment">           .    Delay filling in the last row until the
</span><span class="comment">*</span><span class="comment">           .    vigilant deflation check is complete. ====
</span><span class="comment">*</span><span class="comment">
</span>            IF( ACCUM ) THEN
               JTOP = MAX( KTOP, INCOL )
            ELSE IF( WANTT ) THEN
               JTOP = 1
            ELSE
               JTOP = KTOP
            END IF
            DO 90 M = MTOP, MBOT
               IF( V( 1, M ).NE.ZERO ) THEN
                  K = KRCOL + 3*( M-1 )
                  DO 60 J = JTOP, MIN( KBOT, K+3 )
                     REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
     $                        H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
                     H( J, K+1 ) = H( J, K+1 ) - REFSUM
                     H( J, K+2 ) = H( J, K+2 ) - REFSUM*V( 2, M )
                     H( J, K+3 ) = H( J, K+3 ) - REFSUM*V( 3, M )
   60             CONTINUE
<span class="comment">*</span><span class="comment">
</span>                  IF( ACCUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Accumulate U. (If necessary, update Z later
</span><span class="comment">*</span><span class="comment">                    .    with with an efficient matrix-matrix
</span><span class="comment">*</span><span class="comment">                    .    multiply.) ====
</span><span class="comment">*</span><span class="comment">
</span>                     KMS = K - INCOL
                     DO 70 J = MAX( 1, KTOP-INCOL ), KDU
                        REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
     $                           U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
                        U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
                        U( J, KMS+2 ) = U( J, KMS+2 ) - REFSUM*V( 2, M )
                        U( J, KMS+3 ) = U( J, KMS+3 ) - REFSUM*V( 3, M )
   70                CONTINUE
                  ELSE IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== U is not accumulated, so update Z
</span><span class="comment">*</span><span class="comment">                    .    now by multiplying by reflections
</span><span class="comment">*</span><span class="comment">                    .    from the right. ====
</span><span class="comment">*</span><span class="comment">
</span>                     DO 80 J = ILOZ, IHIZ
                        REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
     $                           Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
                        Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
                        Z( J, K+2 ) = Z( J, K+2 ) - REFSUM*V( 2, M )
                        Z( J, K+3 ) = Z( J, K+3 ) - REFSUM*V( 3, M )
   80                CONTINUE
                  END IF
               END IF
   90       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Special case: 2-by-2 reflection (if needed) ====
</span><span class="comment">*</span><span class="comment">