dsbgst.f.html

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

         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 170 K = KB, 1, -1
            IF( UPDATE ) THEN
               J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
            ELSE
               J2 = I - K - 1 + MAX( 1, K-I0+1 )*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 left
</span><span class="comment">*</span><span class="comment">
</span>            DO 140 L = KB - K, 1, -1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="DLARTV.408"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L, J2-L+1 ), INCA,
     $                         AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
     $                         WORK( J2-KA ), KA1 )
  140       CONTINUE
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            DO 150 J = J1, J2, -KA1
               WORK( J ) = WORK( J-KA )
               WORK( N+J ) = WORK( N+J-KA )
  150       CONTINUE
            DO 160 J = J2, 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)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J ) = WORK( J )*AB( 1, J+1 )
               AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
  160       CONTINUE
            IF( UPDATE ) THEN
               IF( I-K.LT.N-KA .AND. K.LE.KBT )
     $            WORK( I-K+KA ) = WORK( I-K )
            END IF
  170    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         DO 210 K = KB, 1, -1
            J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            IF( NR.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              generate rotations in 2nd set to annihilate elements
</span><span class="comment">*</span><span class="comment">              which have been created outside the band
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="DLARGV.441"></a><a href="dlargv.f.html#DLARGV.1">DLARGV</a>( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
     $                      WORK( N+J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 2nd set from the right
</span><span class="comment">*</span><span class="comment">
</span>               DO 180 L = 1, KA - 1
                  CALL <a name="DLARTV.447"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  180          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              apply rotations in 2nd 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="DLAR2V.455"></a><a href="dlar2v.f.html#DLAR2V.1">DLAR2V</a>( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, WORK( N+J2 ),
     $                      WORK( J2 ), KA1 )
<span class="comment">*</span><span class="comment">
</span>            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           start applying rotations in 2nd set from the left
</span><span class="comment">*</span><span class="comment">
</span>            DO 190 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="DLARTV.466"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  190       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 2nd set
</span><span class="comment">*</span><span class="comment">
</span>               DO 200 J = J2, J1, KA1
                  CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  200          CONTINUE
            END IF
  210    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         DO 230 K = 1, KB - 1
            J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           finish applying rotations in 1st set from the left
</span><span class="comment">*</span><span class="comment">
</span>            DO 220 L = KB - K, 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL <a name="DLARTV.490"></a><a href="dlartv.f.html#DLARTV.1">DLARTV</a>( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA,
     $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
  220       CONTINUE
  230    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( KB.GT.1 ) THEN
            DO 240 J = N - 1, I - KB + 2*KA + 1, -1
               WORK( N+J-M ) = WORK( N+J-KA-M )
               WORK( J-M ) = WORK( J-KA-M )
  240       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Transform A, working with the lower 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))**T * A * inv(S(i))
</span><span class="comment">*</span><span class="comment">
</span>            BII = BB( 1, I )
            DO 250 J = I, I1
               AB( J-I+1, I ) = AB( J-I+1, I ) / BII
  250       CONTINUE
            DO 260 J = MAX( 1, I-KA ), I
               AB( I-J+1, J ) = AB( I-J+1, J ) / BII
  260       CONTINUE
            DO 290 K = I - KBT, I - 1
               DO 270 J = I - KBT, K
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( I-J+1, J )*AB( I-K+1, K ) -
     $                             BB( I-K+1, K )*AB( I-J+1, J ) +
     $                             AB( 1, I )*BB( I-J+1, J )*
     $                             BB( I-K+1, K )
  270          CONTINUE
               DO 280 J = MAX( 1, I-KA ), I - KBT - 1
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( I-K+1, K )*AB( I-J+1, J )
  280          CONTINUE
  290       CONTINUE
            DO 310 J = I, I1
               DO 300 K = MAX( J-KA, I-KBT ), I - 1
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( I-K+1, K )*AB( J-I+1, I )
  300          CONTINUE
  310       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 DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
     $                       BB( KBT+1, I-KBT ), LDBB-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(i1,i) in RA1 for use in next loop over K
</span><span class="comment">*</span><span class="comment">
</span>            RA1 = AB( I1-I+1, 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 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 360 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-k+ka+1,i)
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="DLARTG.568"></a><a href="dlartg.f.html#DLARTG.1">DLARTG</a>( AB( KA1-K, I ), RA1, WORK( N+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+ka+1,i-k) 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( K+1, I-K )*RA1
                  WORK( I-K ) = WORK( N+I-K+KA-M )*T -
     $                          WORK( I-K+KA-M )*AB( KA1, I-K )
                  AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
     $                             WORK( N+I-K+KA-M )*AB( KA1, I-K )
                  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 320 J = J2T, J1, 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-m)
</span><span class="comment">*</span><span class="comment">
</span>               WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
               AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
  320       CONTINUE
<span class="comment">*</span><span class="comment">