ssbgst.f

famous linear algebra library (LAPACK) ports to windows

     $                      KA1, WORK( N+M-KB+J1 ), KA1 )
*
*              apply rotations in 2nd set from the left
*
               DO 660 L = 1, KA - 1
                  CALL SLARTV( NR, AB( KA1-L, J1+L ), INCA,
     $                         AB( KA-L, J1+L ), INCA,
     $                         WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
  660          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
     $                      AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
     $                      WORK( M-KB+J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the right
*
            DO 670 L = KA - 1, KB - K + 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA,
     $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
     $                         KA1 )
  670       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 680 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
  680          CONTINUE
            END IF
  690    CONTINUE
*
         DO 710 K = 1, KB - 1
            J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
*
*           finish applying rotations in 1st set from the right
*
            DO 700 L = KB - K, 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( L, J1T ), INCA,
     $                         AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
     $                         WORK( J1T ), KA1 )
  700       CONTINUE
  710    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
               WORK( N+J ) = WORK( N+J+KA )
               WORK( J ) = WORK( J+KA )
  720       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 730 J = I1, I
               AB( I-J+1, J ) = AB( I-J+1, J ) / BII
  730       CONTINUE
            DO 740 J = I, MIN( N, I+KA )
               AB( J-I+1, I ) = AB( J-I+1, I ) / BII
  740       CONTINUE
            DO 770 K = I + 1, I + KBT
               DO 750 J = K, I + KBT
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( J-I+1, I )*AB( K-I+1, I ) -
     $                             BB( K-I+1, I )*AB( J-I+1, I ) +
     $                             AB( 1, I )*BB( J-I+1, I )*
     $                             BB( K-I+1, I )
  750          CONTINUE
               DO 760 J = I + KBT + 1, MIN( N, I+KA )
                  AB( J-K+1, K ) = AB( J-K+1, K ) -
     $                             BB( K-I+1, I )*AB( J-I+1, I )
  760          CONTINUE
  770       CONTINUE
            DO 790 J = I1, I
               DO 780 K = I + 1, MIN( J+KA, I+KBT )
                  AB( K-J+1, J ) = AB( K-J+1, J ) -
     $                             BB( K-I+1, I )*AB( I-J+1, J )
  780          CONTINUE
  790       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by inv(S(i))
*
               CALL SSCAL( NX, ONE / BII, X( 1, I ), 1 )
               IF( KBT.GT.0 )
     $            CALL SGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
     $                       X( 1, I+1 ), LDX )
            END IF
*
*           store a(i,i1) in RA1 for use in next loop over K
*
            RA1 = AB( I-I1+1, I1 )
         END IF
*
*        Generate and apply vectors of rotations to chase all the
*        existing bulges KA positions up toward the top of the band
*
         DO 840 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-KA1.GT.0 .AND. I+K.LT.M ) THEN
*
*                 generate rotation to annihilate a(i,i+k-ka-1)
*
                  CALL SLARTG( AB( KA1-K, I+K-KA ), RA1,
     $                         WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
*
*                 create nonzero element a(i+k,i+k-ka-1) outside the
*                 band and store it in WORK(m-kb+i+k)
*
                  T = -BB( K+1, I )*RA1
                  WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
     $                               WORK( I+K-KA )*AB( KA1, I+K-KA )
                  AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
     $                                WORK( N+I+K-KA )*AB( KA1, I+K-KA )
                  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 800 J = J1, J2T, KA1
*
*              create nonzero element a(j+ka,j-1) outside the band
*              and store it in WORK(j)
*
               WORK( J ) = WORK( J )*AB( KA1, J-1 )
               AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
  800       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, J1 ), INCA, WORK( J1 ), KA1,
     $                      WORK( N+J1 ), KA1 )
            IF( NR.GT.0 ) THEN
*
*              apply rotations in 1st set from the right
*
               DO 810 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
     $                         INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
  810          CONTINUE
*
*              apply rotations in 1st set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
     $                      AB( 2, J1-1 ), INCA, WORK( N+J1 ),
     $                      WORK( J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 1st set from the left
*
            DO 820 L = KA - 1, KB - K + 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+J1T ), WORK( J1T ), KA1 )
  820       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 1st set
*
               DO 830 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+J ), WORK( J ) )
  830          CONTINUE
            END IF
  840    CONTINUE
*
         IF( UPDATE ) THEN
            IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
*
*              create nonzero element a(i+kbt,i+kbt-ka-1) outside the
*              band and store it in WORK(m-kb+i+kbt)
*
               WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
            END IF
         END IF
*
         DO 880 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
*
*           finish applying rotations in 2nd set from the left
*
            DO 850 L = KB - K, 1, -1
               NRT = ( J2+KA+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
     $                         AB( KA1-L, J1T+L-1 ), INCA,
     $                         WORK( N+M-KB+J1T+KA ),
     $                         WORK( M-KB+J1T+KA ), KA1 )
  850       CONTINUE
            NR = ( J2+KA-1 ) / KA1
            J1 = J2 - ( NR-1 )*KA1
            DO 860 J = J1, J2, KA1
               WORK( M-KB+J ) = WORK( M-KB+J+KA )
               WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
  860       CONTINUE
            DO 870 J = J1, J2, KA1
*
*              create nonzero element a(j+ka,j-1) outside the band
*              and store it in WORK(m-kb+j)
*
               WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
               AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
  870       CONTINUE
            IF( UPDATE ) THEN
               IF( I+K.GT.KA1 .AND. K.LE.KBT )
     $            WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
            END IF
  880    CONTINUE
*
         DO 920 K = KB, 1, -1
            J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
            NR = ( J2+KA-1 ) / 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( KA1, J1 ), INCA, WORK( M-KB+J1 ),
     $                      KA1, WORK( N+M-KB+J1 ), KA1 )
*
*              apply rotations in 2nd set from the right
*
               DO 890 L = 1, KA - 1
                  CALL SLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
     $                         INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
     $                         KA1 )
  890          CONTINUE
*
*              apply rotations in 2nd set from both sides to diagonal
*              blocks
*
               CALL SLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
     $                      AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
     $                      WORK( M-KB+J1 ), KA1 )
*
            END IF
*
*           start applying rotations in 2nd set from the left
*
            DO 900 L = KA - 1, KB - K + 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
     $                         KA1 )
  900       CONTINUE
*
            IF( WANTX ) THEN
*
*              post-multiply X by product of rotations in 2nd set
*
               DO 910 J = J1, J2, KA1
                  CALL SROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
     $                       WORK( N+M-KB+J ), WORK( M-KB+J ) )
  910          CONTINUE
            END IF
  920    CONTINUE
*
         DO 940 K = 1, KB - 1
            J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
*
*           finish applying rotations in 1st set from the left
*
            DO 930 L = KB - K, 1, -1
               NRT = ( J2+L-1 ) / KA1
               J1T = J2 - ( NRT-1 )*KA1
               IF( NRT.GT.0 )
     $            CALL SLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
     $                         AB( KA1-L, J1T-KA1+L ), INCA,
     $                         WORK( N+J1T ), WORK( J1T ), KA1 )
  930       CONTINUE
  940    CONTINUE
*
         IF( KB.GT.1 ) THEN
            DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
               WORK( N+J ) = WORK( N+J+KA )
               WORK( J ) = WORK( J+KA )
  950       CONTINUE
         END IF
*
      END IF
*
      GO TO 490
*
*     End of SSBGST
*
      END