ssbgst.f

famous linear algebra library (LAPACK) ports to windows

   90       CONTINUE
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
     $                      WORK( N+J2T-M ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the right
*
               DO 100 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  100          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
     $                      WORK( J2-M ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the left
*
            DO 110 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA,
     $                         WORK( N+J2-M ), WORK( J2-M ), KA1 )
  110       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 120 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J-M ), WORK( J-M ) )
  120          CONTINUE
            END IF
  130    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.LE.N .AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i-kbt,i-kbt+ka+1) outside the
*              band and store it in WORK(i-kbt)
*
               WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
            END IF
         END IF
*
         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
*
*           finish applying rotations in 2nd set from the left
*
            DO 140 L = KB - K, 1, -1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( 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
*
*              create nonzero element a(j-ka,j+1) outside the band
*              and store it in WORK(j)
*
               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
*
         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
*
*              generate rotations in 2nd set to annihilate elements
*              which have been created outside the band
*
               CALL SLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
     $                      WORK( N+J2 ), KA1 )
*
*              apply rotations in 2nd set from the right
*
               DO 180 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J2 ), INCA,
     $                         AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  180          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
     $                      AB( KA, J2+1 ), INCA, WORK( N+J2 ),
     $                      WORK( J2 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the left
*
            DO 190 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J2+KA1-L ), INCA,
     $                         AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
     $                         WORK( J2 ), KA1 )
  190       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 200 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  200          CONTINUE
            END IF
  210    CONTINUE
*
         DO 230 K = 1, KB - 1
            J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
*
*           finish applying rotations in 1st set from the left
*
            DO 220 L = KB - K, 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( 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
*
         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
*
      ELSE
*
*        Transform A, working with the lower triangle
*
         IF( UPDATE ) THEN
*
*           Form  inv(S(i))**T * A * inv(S(i))
*
            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
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( N-M, KBT, -ONE, X( M+1, I ), 1,
     $                       BB( KBT+1, I-KBT ), LDBB-1,
     $                       X( M+1, I-KBT ), LDX )
            END IF
*
*           store a(i1,i) in RA1 for use in next loop over K
*
            RA1 = AB( I1-I+1, I )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions down toward the bottom of the
*        band
*
         DO 360 K = 1, KB - 1
            IF( UPDATE ) THEN
*
*              Determine the rotations which would annihilate the bulge
*              which has in theory just been created
*
               IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
*
*                 generate rotation to annihilate a(i-k+ka+1,i)
*
                  CALL SLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
     $                         WORK( I-K+KA-M ), RA )
*
*                 create nonzero element a(i-k+ka+1,i-k) outside the
*                 band and store it in WORK(i-k)
*
                  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
*
*              create nonzero element a(j+1,j-ka) outside the band
*              and store it in WORK(j-m)
*
               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
*
*           generate rotations in 1st set to annihilate elements which
*           have been created outside the band
*
            IF( NRT.GT.0 )
     $         CALL SLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
     $                      KA1, WORK( N+J2T-M ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the left
*
               DO 330 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J2-L ), INCA,
     $                         AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  330          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
     $                      INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the right
*
            DO 340 L = KA - 1, KB - K + 1, -1
               NRT = ( N-J2+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
     $                         AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
     $                         WORK( J2-M ), KA1 )
  340       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 350 J = J2, J1, KA1
                  CALL SROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
     $                       WORK( N+J-M ), WORK( J-M ) )
  350          CONTINUE
            END IF
  360    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.LE.N .AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i-kbt+ka+1,i-kbt) outside the
*              band and store it in WORK(i-kbt)
*
               WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
            END IF
         END IF
*
         DO 400 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
*
*           finish applying rotations in 2nd set from the right
*
            DO 370 L = KB - K, 1, -1
               NRT = ( N-J2+KA+L ) / KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
     $                         AB( KA1-L, J2-KA+1 ), INCA,
     $                         WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
  370       CONTINUE
            NR = ( N-J2+KA ) / KA1
            J1 = J2 + ( NR-1 )*KA1
            DO 380 J = J1, J2, -KA1
               WORK( J ) = WORK( J-KA )