chbgst.f.html

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

</span><span class="comment">*</span><span class="comment">     A in phase 1, and up toward the top of A in phase 2, by applying
</span><span class="comment">*</span><span class="comment">     plane rotations.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
</span><span class="comment">*</span><span class="comment">     of them are linearly independent, so annihilating a bulge requires
</span><span class="comment">*</span><span class="comment">     only 2*kb-1 plane rotations. The rotations are divided into a 1st
</span><span class="comment">*</span><span class="comment">     set of kb-1 rotations, and a 2nd set of kb rotations.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Wherever possible, rotations are generated and applied in vector
</span><span class="comment">*</span><span class="comment">     operations of length NR between the indices J1 and J2 (sometimes
</span><span class="comment">*</span><span class="comment">     replaced by modified values NRT, J1T or J2T).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The real cosines and complex sines of the rotations are stored in
</span><span class="comment">*</span><span class="comment">     the arrays RWORK and WORK, those of the 1st set in elements
</span><span class="comment">*</span><span class="comment">     2:m-kb-1, and those of the 2nd set in elements m-kb+1:n.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The bulges are not formed explicitly; nonzero elements outside the
</span><span class="comment">*</span><span class="comment">     band are created only when they are required for generating new
</span><span class="comment">*</span><span class="comment">     rotations; they are stored in the array WORK, in positions where
</span><span class="comment">*</span><span class="comment">     they are later overwritten by the sines of the rotations which
</span><span class="comment">*</span><span class="comment">     annihilate them.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     **************************** Phase 1 *****************************
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The logical structure of this phase is:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     UPDATE = .TRUE.
</span><span class="comment">*</span><span class="comment">     DO I = N, M + 1, -1
</span><span class="comment">*</span><span class="comment">        use S(i) to update A and create a new bulge
</span><span class="comment">*</span><span class="comment">        apply rotations to push all bulges KA positions downward
</span><span class="comment">*</span><span class="comment">     END DO
</span><span class="comment">*</span><span class="comment">     UPDATE = .FALSE.
</span><span class="comment">*</span><span class="comment">     DO I = M + KA + 1, N - 1
</span><span class="comment">*</span><span class="comment">        apply rotations to push all bulges KA positions downward
</span><span class="comment">*</span><span class="comment">     END DO
</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 = N + 1
   10 CONTINUE
      IF( UPDATE ) THEN
         I = I - 1
         KBT = MIN( KB, I-1 )
         I0 = I - 1
         I1 = MIN( N, I+KA )
         I2 = I - KBT + KA1
         IF( I.LT.M+1 ) THEN
            UPDATE = .FALSE.
            I = I + 1
            I0 = M
            IF( KA.EQ.0 )
     $         GO TO 480
            GO TO 10
         END IF
      ELSE
         I = I + KA
         IF( I.GT.N-1 )
     $      GO TO 480
      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 = REAL( BB( KB1, I ) )
            AB( KA1, I ) = ( REAL( AB( KA1, I ) ) / BII ) / BII
            DO 20 J = I + 1, I1
               AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
   20       CONTINUE
            DO 30 J = MAX( 1, I-KA ), I - 1
               AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
   30       CONTINUE
            DO 60 K = I - KBT, I - 1
               DO 40 J = I - KBT, K
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               BB( J-I+KB1, I )*
     $                               CONJG( AB( K-I+KA1, I ) ) -
     $                               CONJG( BB( K-I+KB1, I ) )*
     $                               AB( J-I+KA1, I ) +
     $                               REAL( AB( KA1, I ) )*
     $                               BB( J-I+KB1, I )*
     $                               CONJG( BB( K-I+KB1, I ) )
   40          CONTINUE
               DO 50 J = MAX( 1, I-KA ), I - KBT - 1
                  AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
     $                               CONJG( BB( K-I+KB1, I ) )*
     $                               AB( J-I+KA1, I )
   50          CONTINUE
   60       CONTINUE
            DO 80 J = I, I1
               DO 70 K = MAX( J-KA, I-KBT ), I - 1
                  AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
     $                               BB( K-I+KB1, I )*AB( I-J+KA1, J )
   70          CONTINUE
   80       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 CSSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL CGERC( N-M, KBT, -CONE, X( M+1, I ), 1,
     $                        BB( KB1-KBT, I ), 1, X( M+1, I-KBT ),
     $                        LDX )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           store a(i,i1) in RA1 for use in next loop over K
</span><span class="comment">*</span><span class="comment">
</span>            RA1 = AB( I-I1+KA1, I1 )
         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 down toward the bottom of the
</span><span class="comment">*</span><span class="comment">        band
</span><span class="comment">*</span><span class="comment">
</span>         DO 130 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+KA.LT.N .AND. I-K.GT.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 generate rotation to annihilate a(i,i-k+ka+1)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLARTG.317"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>( AB( K+1, I-K+KA ), RA1,
     $                         RWORK( I-K+KA-M ), WORK( I-K+KA-M ), RA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 create nonzero element a(i-k,i-k+ka+1) outside the
</span><span class="comment">*</span><span class="comment">                 band and store it in WORK(i-k)
</span><span class="comment">*</span><span class="comment">
</span>                  T = -BB( KB1-K, I )*RA1
                  WORK( I-K ) = RWORK( I-K+KA-M )*T -
     $                          CONJG( WORK( I-K+KA-M ) )*
     $                          AB( 1, I-K+KA )
                  AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
     $                              RWORK( I-K+KA-M )*AB( 1, I-K+KA )
                  RA1 = RA
               END IF
            END IF
            J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( UPDATE ) THEN
               J2T = MAX( J2, I+2*KA-K+1 )
            ELSE
               J2T = J2
            END IF
            NRT = ( N-J2T+KA ) / KA1
            DO 90 J = J2T, J1, KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              create nonzero element a(j-ka,j+1) outside the band
</span><span class="comment">*</span><span class="comment">              and store it in WORK(j-m)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
               AB( 1, J+1 ) = RWORK( J-M )*AB( 1, J+1 )
   90       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="CLARGV.354"></a><a href="clargv.f.html#CLARGV.1">CLARGV</a>( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
     $                      RWORK( J2T-M ), 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 right
</span><span class="comment">*</span><span class="comment">
</span>               DO 100 L = 1, KA - 1
                  CALL <a name="CLARTV.361"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, RWORK( J2-M ),
     $                         WORK( J2-M ), KA1 )
  100          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="CLAR2V.369"></a><a href="clar2v.f.html#CLAR2V.1">CLAR2V</a>( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, RWORK( J2-M ),
     $                      WORK( J2-M ), KA1 )
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACGV.373"></a><a href="clacgv.f.html#CLACGV.1">CLACGV</a>( NR, WORK( J2-M ), 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 left
</span><span class="comment">*</span><span class="comment">
</span>            DO 110 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="CLARTV.381"></a><a href="clartv.f.html#CLARTV.1">CLARTV</a>( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA, RWORK( J2-M ),
     $                         WORK( J2-M ), KA1 )
  110       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">