zhbgst.f.html

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

</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     To avoid duplicating code, the two loops are merged.
</span><span class="comment">*</span><span class="comment">
</span>      UPDATE = .TRUE.
      I = 0
  490 CONTINUE
      IF( UPDATE ) THEN
         I = I + 1
         KBT = MIN( KB, M-I )
         I0 = I + 1
         I1 = MAX( 1, I-KA )
         I2 = I + KBT - KA1
         IF( I.GT.M ) THEN
            UPDATE = .FALSE.
            I = I - 1
            I0 = M + 1
            IF( KA.EQ.0 )
     $         RETURN
            GO TO 490
         END IF
      ELSE
         I = I - KA
         IF( I.LT.2 )
     $      RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( I.LT.M-KBT ) THEN
         NX = M
      ELSE
         NX = N
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Transform A, working with the upper triangle
</span><span class="comment">*</span><span class="comment">
</span>         IF( UPDATE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Form  inv(S(i))**H * A * inv(S(i))
</span><span class="comment">*</span><span class="comment">
</span>            BII = DBLE( BB( KB1, I ) )
            AB( KA1, I ) = ( DBLE( AB( KA1, I ) ) / BII ) / BII
            DO 500 J = I1, I - 1
               AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
  500       CONTINUE
            DO 510 J = I + 1, MIN( N, I+KA )
               AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
  510       CONTINUE
            DO 540 K = I + 1, I + KBT
               DO 520 J = K, I + KBT
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( I-J+KB1, J )*
     $                               DCONJG( AB( I-K+KA1, K ) ) -
     $                               DCONJG( BB( I-K+KB1, K ) )*
     $                               AB( I-J+KA1, J ) +
     $                               DBLE( AB( KA1, I ) )*
     $                               BB( I-J+KB1, J )*
     $                               DCONJG( BB( I-K+KB1, K ) )
  520          CONTINUE
               DO 530 J = I + KBT + 1, MIN( N, I+KA )
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               DCONJG( BB( I-K+KB1, K ) )*
     $                               AB( I-J+KA1, J )
  530          CONTINUE
  540       CONTINUE
            DO 560 J = I1, I
               DO 550 K = I + 1, MIN( J+KA, I+KBT )
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( I-K+KB1, K )*AB( J-I+KA1, I )
  550          CONTINUE
  560       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              post-multiply X by inv(S(i))
</span><span class="comment">*</span><span class="comment">
</span>               CALL ZDSCAL( NX, ONE / BII, X( 1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL ZGERU( NX, KBT, -CONE, X( 1, I ), 1,
     $                        BB( KB, I+1 ), LDBB-1, X( 1, I+1 ), LDX )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           store a(i1,i) in RA1 for use in next loop over K
</span><span class="comment">*</span><span class="comment">
</span>            RA1 = AB( I1-I+KA1, I )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate and apply vectors of rotations to chase all the
</span><span class="comment">*</span><span class="comment">        existing bulges KA positions up toward the top of the band
</span><span class="comment">*</span><span class="comment">
</span>         DO 610 K = 1, KB - 1
            IF( UPDATE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Determine the rotations which would annihilate the bulge
</span><span class="comment">*</span><span class="comment">              which has in theory just been created
</span><span class="comment">*</span><span class="comment">
</span>               IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 generate rotation to annihilate a(i+k-ka-1,i)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="ZLARTG.898"></a><a href="zlartg.f.html#ZLARTG.1">ZLARTG</a>( AB( K+1, I ), RA1, RWORK( I+K-KA ),
     $                         WORK( I+K-KA ), RA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 create nonzero element a(i+k-ka-1,i+k) outside the
</span><span class="comment">*</span><span class="comment">                 band and store it in WORK(m-kb+i+k)
</span><span class="comment">*</span><span class="comment">
</span>                  T = -BB( KB1-K, I+K )*RA1
                  WORK( M-KB+I+K ) = RWORK( I+K-KA )*T -
     $                               DCONJG( WORK( I+K-KA ) )*
     $                               AB( 1, I+K )
                  AB( 1, I+K ) = WORK( I+K-KA )*T +
     $                           RWORK( I+K-KA )*AB( 1, I+K )
                  RA1 = RA
               END IF
            END IF
            J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MIN( J2, I-2*KA+K-1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( J2T+KA-1 ) / KA1
            DO 570 J = J1, J2T, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(j-1,j+ka) outside the band
</span><span class="comment">*</span><span class="comment">              and store it in WORK(j)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
               AB( 1, J+KA-1 ) = RWORK( J )*AB( 1, J+KA-1 )
  570       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           generate rotations in 1st set to annihilate elements which
</span><span class="comment">*</span><span class="comment">           have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span>            IF( NRT.GT.0 )
     $         CALL <a name="ZLARGV.935"></a><a href="zlargv.f.html#ZLARGV.1">ZLARGV</a>( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
     $                      RWORK( J1 ), KA1 )
            IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 1st set from the left
</span><span class="comment">*</span><span class="comment">
</span>               DO 580 L = 1, KA - 1
                  CALL <a name="ZLARTV.942"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NR, AB( KA1-L, J1+L ), INCA,
     $                         AB( KA-L, J1+L ), INCA, RWORK( J1 ),
     $                         WORK( J1 ), KA1 )
  580          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 1st set from both sides to diagonal
</span><span class="comment">*</span><span class="comment">              blocks
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLAR2V.950"></a><a href="zlar2v.f.html#ZLAR2V.1">ZLAR2V</a>( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
     $                      AB( KA, J1 ), INCA, RWORK( J1 ), WORK( J1 ),
     $                      KA1 )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="ZLACGV.954"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( NR, WORK( J1 ), KA1 )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           start applying rotations in 1st set from the right
</span><span class="comment">*</span><span class="comment">
</span>            DO 590 L = KA - 1, KB - K + 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL <a name="ZLARTV.963"></a><a href="zlartv.f.html#ZLARTV.1">ZLARTV</a>( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA, RWORK( J1T ),
     $                         WORK( J1T ), KA1 )
  590       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( WANTX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              post-multiply X by product of rotations in 1st set
</span><span class="comment">*</span><span class="comment">
</span>               DO 600 J = J1, J2, KA1
                  CALL <a name="ZROT.973"></a><a href="zrot.f.html#ZROT.1">ZROT</a>( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       RWORK( J ), WORK( J ) )
  600          CONTINUE
            END IF
  610    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( UPDATE ) THEN
            IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(i+kbt-ka-1,i+kbt) outside the
</span><span class="comment">*</span><span class="comment">              band and store it in WORK(m-kb+i+kbt)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 650 K = KB, 1, -1
            IF( UPDATE ) THEN
               J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
            ELSE
               J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           finish applying rotations in 2nd set from the right
</span><span class="comment">*</span><span class="com