claqr5.f.html

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

</span>                     IF( CABS1( H( K+1, K ) )*
     $                   ( CABS1( VT( 2 ) )+CABS1( VT( 3 ) ) ).GT.ULP*
     $                   CABS1( VT( 1 ) )*( CABS1( H( K,
     $                   K ) )+CABS1( H( K+1, K+1 ) )+CABS1( 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="CLARFG.345"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
                        REFSUM = H( K+1, K ) +
     $                           H( K+2, K )*CONJG( VT( 2 ) ) +
     $                           H( K+3, K )*CONJG( VT( 3 ) )
                        H( K+1, K ) = H( K+1, K ) -
     $                                CONJG( 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
   10       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="CLAQR1.366"></a><a href="claqr1.f.html#CLAQR1.1">CLAQR1</a>( 2, H( K+1, K+1 ), LDH, S( 2*M22-1 ),
     $                         S( 2*M22 ), V( 1, M22 ) )
                  BETA = V( 1, M22 )
                  CALL <a name="CLARFG.369"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</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="CLARFG.373"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</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 30 J = MAX( KTOP, KRCOL ), JBOT
               MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
               DO 20 M = MTOP, MEND
                  K = KRCOL + 3*( M-1 )
                  REFSUM = CONJG( V( 1, M ) )*
     $                     ( H( K+1, J )+CONJG( V( 2, M ) )*H( K+2, J )+
     $                     CONJG( 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 )
   20          CONTINUE
   30       CONTINUE
            IF( BMP22 ) THEN
               K = KRCOL + 3*( M22-1 )
               DO 40 J = MAX( K+1, KTOP ), JBOT
                  REFSUM = CONJG( V( 1, M22 ) )*
     $                     ( H( K+1, J )+CONJG( 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 )
   40          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 80 M = MTOP, MBOT
               IF( V( 1, M ).NE.ZERO ) THEN
                  K = KRCOL + 3*( M-1 )
                  DO 50 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*CONJG( V( 2, M ) )
                     H( J, K+3 ) = H( J, K+3 ) -
     $                             REFSUM*CONJG( V( 3, M ) )
   50             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 60 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*CONJG( V( 2, M ) )
                        U( J, KMS+3 ) = U( J, KMS+3 ) -
     $                                  REFSUM*CONJG( V( 3, M ) )
   60                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 70 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*CONJG( V( 2, M ) )
                        Z( J, K+3 ) = Z( J, K+3 ) -
     $                                REFSUM*CONJG( V( 3, M ) )
   70                CONTINUE
                  END IF
               END IF
   80       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">
</span>            K = KRCOL + 3*( M22-1 )
            IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
               DO 90 J = JTOP, MIN( KBOT, K+3 )
                  REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
     $                     H( J, K+2 ) )