#include "blaswrap.h"
#include "f2c.h"

/* Subroutine */ int dsbgst_(char *vect, char *uplo, integer *n, integer *ka, 
	integer *kb, doublereal *ab, integer *ldab, doublereal *bb, integer *
	ldbb, doublereal *x, integer *ldx, doublereal *work, integer *info)
{
/*  -- LAPACK routine (version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       June 30, 1999   


    Purpose   
    =======   

    DSBGST reduces a real symmetric-definite banded generalized   
    eigenproblem  A*x = lambda*B*x  to standard form  C*y = lambda*y,   
    such that C has the same bandwidth as A.   

    B must have been previously factorized as S**T*S by DPBSTF, using a   
    split Cholesky factorization. A is overwritten by C = X**T*A*X, where   
    X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the   
    bandwidth of A.   

    Arguments   
    =========   

    VECT    (input) CHARACTER*1   
            = 'N':  do not form the transformation matrix X;   
            = 'V':  form X.   

    UPLO    (input) CHARACTER*1   
            = 'U':  Upper triangle of A is stored;   
            = 'L':  Lower triangle of A is stored.   

    N       (input) INTEGER   
            The order of the matrices A and B.  N >= 0.   

    KA      (input) INTEGER   
            The number of superdiagonals of the matrix A if UPLO = 'U',   
            or the number of subdiagonals if UPLO = 'L'.  KA >= 0.   

    KB      (input) INTEGER   
            The number of superdiagonals of the matrix B if UPLO = 'U',   
            or the number of subdiagonals if UPLO = 'L'.  KA >= KB >= 0.   

    AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)   
            On entry, the upper or lower triangle of the symmetric band   
            matrix A, stored in the first ka+1 rows of the array.  The   
            j-th column of A is stored in the j-th column of the array AB   
            as follows:   
            if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;   
            if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).   

            On exit, the transformed matrix X**T*A*X, stored in the same   
            format as A.   

    LDAB    (input) INTEGER   
            The leading dimension of the array AB.  LDAB >= KA+1.   

    BB      (input) DOUBLE PRECISION array, dimension (LDBB,N)   
            The banded factor S from the split Cholesky factorization of   
            B, as returned by DPBSTF, stored in the first KB+1 rows of   
            the array.   

    LDBB    (input) INTEGER   
            The leading dimension of the array BB.  LDBB >= KB+1.   

    X       (output) DOUBLE PRECISION array, dimension (LDX,N)   
            If VECT = 'V', the n-by-n matrix X.   
            If VECT = 'N', the array X is not referenced.   

    LDX     (input) INTEGER   
            The leading dimension of the array X.   
            LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.   

    WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)   

    INFO    (output) INTEGER   
            = 0:  successful exit   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   

    =====================================================================   


       Test the input parameters   

       Parameter adjustments */
    /* Table of constant values */
    static doublereal c_b8 = 0.;
    static doublereal c_b9 = 1.;
    static integer c__1 = 1;
    static doublereal c_b20 = -1.;
    
    /* System generated locals */
    integer ab_dim1, ab_offset, bb_dim1, bb_offset, x_dim1, x_offset, i__1, 
	    i__2, i__3, i__4;
    doublereal d__1;
    /* Local variables */
    static integer inca;
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *), drot_(integer *, doublereal *, integer *, doublereal *
	    , integer *, doublereal *, doublereal *);
    static integer i__, j, k, l, m;
    static doublereal t;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
	    integer *);
    extern logical lsame_(char *, char *);
    static integer i0, i1;
    static logical upper;
    static integer i2, j1, j2;
    static logical wantx;
    extern /* Subroutine */ int dlar2v_(integer *, doublereal *, doublereal *,
	     doublereal *, integer *, doublereal *, doublereal *, integer *);
    static doublereal ra;
    static integer nr, nx;
    extern /* Subroutine */ int dlaset_(char *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *), 
	    dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *), xerbla_(char *, integer *), dlargv_(
	    integer *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static logical update;
    extern /* Subroutine */ int dlartv_(integer *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, doublereal *, integer *);
    static integer ka1, kb1;
    static doublereal ra1;
    static integer j1t, j2t;
    static doublereal bii;
    static integer kbt, nrt;
#define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1]
#define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1]
#define bb_ref(a_1,a_2) bb[(a_2)*bb_dim1 + a_1]


    ab_dim1 = *ldab;
    ab_offset = 1 + ab_dim1 * 1;
    ab -= ab_offset;
    bb_dim1 = *ldbb;
    bb_offset = 1 + bb_dim1 * 1;
    bb -= bb_offset;
    x_dim1 = *ldx;
    x_offset = 1 + x_dim1 * 1;
    x -= x_offset;
    --work;

    /* Function Body */
    wantx = lsame_(vect, "V");
    upper = lsame_(uplo, "U");
    ka1 = *ka + 1;
    kb1 = *kb + 1;
    *info = 0;
    if (! wantx && ! lsame_(vect, "N")) {
	*info = -1;
    } else if (! upper && ! lsame_(uplo, "L")) {
	*info = -2;
    } else if (*n < 0) {
	*info = -3;
    } else if (*ka < 0) {
	*info = -4;
    } else if (*kb < 0) {
	*info = -5;
    } else if (*ldab < *ka + 1) {
	*info = -7;
    } else if (*ldbb < *kb + 1) {
	*info = -9;
    } else if (*ldx < 1 || wantx && *ldx < max(1,*n)) {
	*info = -11;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("DSBGST", &i__1);
	return 0;
    }

/*     Quick return if possible */

    if (*n == 0) {
	return 0;
    }

    inca = *ldab * ka1;

/*     Initialize X to the unit matrix, if needed */

    if (wantx) {
	dlaset_("Full", n, n, &c_b8, &c_b9, &x[x_offset], ldx);
    }

/*     Set M to the splitting point m. It must be the same value as is   
       used in DPBSTF. The chosen value allows the arrays WORK and RWORK   
       to be of dimension (N). */

    m = (*n + *kb) / 2;

/*     The routine works in two phases, corresponding to the two halves   
       of the split Cholesky factorization of B as S**T*S where   

       S = ( U    )   
           ( M  L )   

       with U upper triangular of order m, and L lower triangular of   
       order n-m. S has the same bandwidth as B.   

       S is treated as a product of elementary matrices:   

       S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)   

       where S(i) is determined by the i-th row of S.   

       In phase 1, the index i takes the values n, n-1, ... , m+1;   
       in phase 2, it takes the values 1, 2, ... , m.   

       For each value of i, the current matrix A is updated by forming   
       inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside   
       the band of A. The bulge is then pushed down toward the bottom of   
       A in phase 1, and up toward the top of A in phase 2, by applying   
       plane rotations.   

       There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1   
       of them are linearly independent, so annihilating a bulge requires   
       only 2*kb-1 plane rotations. The rotations are divided into a 1st   
       set of kb-1 rotations, and a 2nd set of kb rotations.   

       Wherever possible, rotations are generated and applied in vector   
       operations of length NR between the indices J1 and J2 (sometimes   
       replaced by modified values NRT, J1T or J2T).   

       The cosines and sines of the rotations are stored in the array   
       WORK. The cosines of the 1st set of rotations are stored in   
       elements n+2:n+m-kb-1 and the sines of the 1st set in elements   
       2:m-kb-1; the cosines of the 2nd set are stored in elements   
       n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.   

       The bulges are not formed explicitly; nonzero elements outside the   
       band are created only when they are required for generating new   
       rotations; they are stored in the array WORK, in positions where   
       they are later overwritten by the sines of the rotations which   
       annihilate them.   

       **************************** Phase 1 *****************************   

       The logical structure of this phase is:   

       UPDATE = .TRUE.   
       DO I = N, M + 1, -1   
          use S(i) to update A and create a new bulge   
          apply rotations to push all bulges KA positions downward   
       END DO   
       UPDATE = .FALSE.   
       DO I = M + KA + 1, N - 1   
          apply rotations to push all bulges KA positions downward   
       END DO   

       To avoid duplicating code, the two loops are merged. */

    update = TRUE_;
    i__ = *n + 1;
L10:
    if (update) {
	--i__;
/* Computing MIN */
	i__1 = *kb, i__2 = i__ - 1;
	kbt = min(i__1,i__2);
	i0 = i__ - 1;
/* Computing MIN */
	i__1 = *n, i__2 = i__ + *ka;
	i1 = min(i__1,i__2);
	i2 = i__ - kbt + ka1;
	if (i__ < m + 1) {
	    update = FALSE_;
	    ++i__;
	    i0 = m;
	    if (*ka == 0) {
		goto L480;
	    }
	    goto L10;
	}
    } else {
	i__ += *ka;
	if (i__ > *n - 1) {
	    goto L480;
	}
    }

    if (upper) {

/*        Transform A, working with the upper triangle */

	if (update) {

/*           Form  inv(S(i))**T * A * inv(S(i)) */

	    bii = bb_ref(kb1, i__);
	    i__1 = i1;
	    for (j = i__; j <= i__1; ++j) {
		ab_ref(i__ - j + ka1, j) = ab_ref(i__ - j + ka1, j) / bii;
/* L20: */
	    }
/* Computing MAX */
	    i__1 = 1, i__2 = i__ - *ka;
	    i__3 = i__;
	    for (j = max(i__1,i__2); j <= i__3; ++j) {
		ab_ref(j - i__ + ka1, i__) = ab_ref(j - i__ + ka1, i__) / bii;
/* L30: */
	    }
	    i__3 = i__ - 1;
	    for (k = i__ - kbt; k <= i__3; ++k) {
		i__1 = k;
		for (j = i__ - kbt; j <= i__1; ++j) {
		    ab_ref(j - k + ka1, k) = ab_ref(j - k + ka1, k) - bb_ref(
			    j - i__ + kb1, i__) * ab_ref(k - i__ + ka1, i__) 
			    - bb_ref(k - i__ + kb1, i__) * ab_ref(j - i__ + 
			    ka1, i__) + ab_ref(ka1, i__) * bb_ref(j - i__ + 
			    kb1, i__) * bb_ref(k - i__ + kb1, i__);
/* L40: */
		}
/* Computing MAX */
		i__1 = 1, i__2 = i__ - *ka;
		i__4 = i__ - kbt - 1;
		for (j = max(i__1,i__2); j <= i__4; ++j) {
		    ab_ref(j - k + ka1, k) = ab_ref(j - k + ka1, k) - bb_ref(
			    k - i__ + kb1, i__) * ab_ref(j - i__ + ka1, i__);
/* L50: */
		}
/* L60: */
	    }
	    i__3 = i1;
	    for (j = i__; j <= i__3; ++j) {
/* Computing MAX */
		i__4 = j - *ka, i__1 = i__ - kbt;
		i__2 = i__ - 1;
		for (k = max(i__4,i__1); k <= i__2; ++k) {
		    ab_ref(k - j + ka1, j) = ab_ref(k - j + ka1, j) - bb_ref(
			    k - i__ + kb1, i__) * ab_ref(i__ - j + ka1, j);
/* L70: */
		}
/* L80: */
	    }

	    if (wantx) {

/*              post-multiply X by inv(S(i)) */

		i__3 = *n - m;
		d__1 = 1. / bii;
		dscal_(&i__3, &d__1, &x_ref(m + 1, i__), &c__1);
		if (kbt > 0) {
		    i__3 = *n - m;
		    dger_(&i__3, &kbt, &c_b20, &x_ref(m + 1, i__), &c__1, &
			    bb_ref(kb1 - kbt, i__), &c__1, &x_ref(m + 1, i__ 
			    - kbt), ldx);
		}
	    }

/*           store a(i,i1) in RA1 for use in next loop over K */

	    ra1 = ab_ref(i__ - i1 + ka1, i1);
	}

/*        Generate and apply vectors of rotations to chase all the   
          existing bulges KA positions down toward the bottom of the   
          band */

	i__3 = *kb - 1;
	for (k = 1; k <= i__3; ++k) {
	    if (update) {

/*              Determine the rotations which would annihilate the bulge   
                which has in theory just been created */

		if (i__ - k + *ka < *n && i__ - k > 1) {

/*                 generate rotation to annihilate a(i,i-k+ka+1) */

		    dlartg_(&ab_ref(k + 1, i__ - k + *ka), &ra1, &work[*n + 
			    i__ - k + *ka - m], &work[i__ - k + *ka - m], &ra)
			    ;

/*                 create nonzero element a(i-k,i-k+ka+1) outside the   
                   band and store it in WORK(i-k) */

		    t = -bb_ref(kb1 - k, i__) * ra1;
		    work[i__ - k] = work[*n + i__ - k + *ka - m] * t - work[
			    i__ - k + *ka - m] * ab_ref(1, i__ - k + *ka);
		    ab_ref(1, i__ - k + *ka) = work[i__ - k + *ka - m] * t + 
			    work[*n + i__ - k + *ka - m] * ab_ref(1, i__ - k 
			    + *ka);
		    ra1 = ra;
		}
	    }
/* Computing MAX */
	    i__2 = 1, i__4 = k - i0 + 2;
	    j2 = i__ - k - 1 + max(i__2,i__4) * ka1;
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    if (update) {
/* Computing MAX */
		i__2 = j2, i__4 = i__ + (*ka << 1) - k + 1;
		j2t = max(i__2,i__4);
	    } else {
		j2t = j2;
	    }
	    nrt = (*n - j2t + *ka) / ka1;
	    i__2 = j1;
	    i__4 = ka1;
	    for (j = j2t; i__4 < 0 ? j >= i__2 : j <= i__2; j += i__4) {

/*              create nonzero element a(j-ka,j+1) outside the band   
                and store it in WORK(j-m) */

		work[j - m] *= ab_ref(1, j + 1);
		ab_ref(1, j + 1) = work[*n + j - m] * ab_ref(1, j + 1);
/* L90: */
	    }

/*           generate rotations in 1st set to annihilate elements which   
             have been created outside the band */

	    if (nrt > 0) {
		dlargv_(&nrt, &ab_ref(1, j2t), &inca, &work[j2t - m], &ka1, &
			work[*n + j2t - m], &ka1);
	    }
	    if (nr > 0) {

/*              apply rotations in 1st set from the right */

		i__4 = *ka - 1;
		for (l = 1; l <= i__4; ++l) {
		    dlartv_(&nr, &ab_ref(ka1 - l, j2), &inca, &ab_ref(*ka - l,
			     j2 + 1), &inca, &work[*n + j2 - m], &work[j2 - m]
			    , &ka1);
/* L100: */
		}

/*              apply rotations in 1st set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(ka1, j2), &ab_ref(ka1, j2 + 1), &ab_ref(*
			ka, j2 + 1), &inca, &work[*n + j2 - m], &work[j2 - m],
			 &ka1);

	    }

/*           start applying rotations in 1st set from the left */

	    i__4 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__4; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j2 + ka1 - l), &inca, &ab_ref(l 
			    + 1, j2 + ka1 - l), &inca, &work[*n + j2 - m], &
			    work[j2 - m], &ka1);
		}
/* L110: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 1st set */

		i__4 = j1;
		i__2 = ka1;
		for (j = j2; i__2 < 0 ? j >= i__4 : j <= i__4; j += i__2) {
		    i__1 = *n - m;
		    drot_(&i__1, &x_ref(m + 1, j), &c__1, &x_ref(m + 1, j + 1)
			    , &c__1, &work[*n + j - m], &work[j - m]);
/* L120: */
		}
	    }
/* L130: */
	}

	if (update) {
	    if (i2 <= *n && kbt > 0) {

/*              create nonzero element a(i-kbt,i-kbt+ka+1) outside the   
                band and store it in WORK(i-kbt) */

		work[i__ - kbt] = -bb_ref(kb1 - kbt, i__) * ra1;
	    }
	}

	for (k = *kb; k >= 1; --k) {
	    if (update) {
/* Computing MAX */
		i__3 = 2, i__2 = k - i0 + 1;
		j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
	    } else {
/* Computing MAX */
		i__3 = 1, i__2 = k - i0 + 1;
		j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
	    }

/*           finish applying rotations in 2nd set from the left */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (*n - j2 + *ka + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j2 - l + 1), &inca, &ab_ref(l + 
			    1, j2 - l + 1), &inca, &work[*n + j2 - *ka], &
			    work[j2 - *ka], &ka1);
		}
/* L140: */
	    }
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    i__3 = j2;
	    i__2 = -ka1;
	    for (j = j1; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
		work[j] = work[j - *ka];
		work[*n + j] = work[*n + j - *ka];
/* L150: */
	    }
	    i__2 = j1;
	    i__3 = ka1;
	    for (j = j2; i__3 < 0 ? j >= i__2 : j <= i__2; j += i__3) {

/*              create nonzero element a(j-ka,j+1) outside the band   
                and store it in WORK(j) */

		work[j] *= ab_ref(1, j + 1);
		ab_ref(1, j + 1) = work[*n + j] * ab_ref(1, j + 1);
/* L160: */
	    }
	    if (update) {
		if (i__ - k < *n - *ka && k <= kbt) {
		    work[i__ - k + *ka] = work[i__ - k];
		}
	    }
/* L170: */
	}

	for (k = *kb; k >= 1; --k) {
/* Computing MAX */
	    i__3 = 1, i__2 = k - i0 + 1;
	    j2 = i__ - k - 1 + max(i__3,i__2) * ka1;
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    if (nr > 0) {

/*              generate rotations in 2nd set to annihilate elements   
                which have been created outside the band */

		dlargv_(&nr, &ab_ref(1, j2), &inca, &work[j2], &ka1, &work[*n 
			+ j2], &ka1);

/*              apply rotations in 2nd set from the right */

		i__3 = *ka - 1;
		for (l = 1; l <= i__3; ++l) {
		    dlartv_(&nr, &ab_ref(ka1 - l, j2), &inca, &ab_ref(*ka - l,
			     j2 + 1), &inca, &work[*n + j2], &work[j2], &ka1);
/* L180: */
		}

/*              apply rotations in 2nd set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(ka1, j2), &ab_ref(ka1, j2 + 1), &ab_ref(*
			ka, j2 + 1), &inca, &work[*n + j2], &work[j2], &ka1);

	    }

/*           start applying rotations in 2nd set from the left */

	    i__3 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__3; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j2 + ka1 - l), &inca, &ab_ref(l 
			    + 1, j2 + ka1 - l), &inca, &work[*n + j2], &work[
			    j2], &ka1);
		}
/* L190: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 2nd set */

		i__3 = j1;
		i__2 = ka1;
		for (j = j2; i__2 < 0 ? j >= i__3 : j <= i__3; j += i__2) {
		    i__4 = *n - m;
		    drot_(&i__4, &x_ref(m + 1, j), &c__1, &x_ref(m + 1, j + 1)
			    , &c__1, &work[*n + j], &work[j]);
/* L200: */
		}
	    }
/* L210: */
	}

	i__2 = *kb - 1;
	for (k = 1; k <= i__2; ++k) {
/* Computing MAX */
	    i__3 = 1, i__4 = k - i0 + 2;
	    j2 = i__ - k - 1 + max(i__3,i__4) * ka1;

/*           finish applying rotations in 1st set from the left */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j2 + ka1 - l), &inca, &ab_ref(l 
			    + 1, j2 + ka1 - l), &inca, &work[*n + j2 - m], &
			    work[j2 - m], &ka1);
		}
/* L220: */
	    }
/* L230: */
	}

	if (*kb > 1) {
	    i__2 = i__ - *kb + (*ka << 1) + 1;
	    for (j = *n - 1; j >= i__2; --j) {
		work[*n + j - m] = work[*n + j - *ka - m];
		work[j - m] = work[j - *ka - m];
/* L240: */
	    }
	}

    } else {

/*        Transform A, working with the lower triangle */

	if (update) {

/*           Form  inv(S(i))**T * A * inv(S(i)) */

	    bii = bb_ref(1, i__);
	    i__2 = i1;
	    for (j = i__; j <= i__2; ++j) {
		ab_ref(j - i__ + 1, i__) = ab_ref(j - i__ + 1, i__) / bii;
/* L250: */
	    }
/* Computing MAX */
	    i__2 = 1, i__3 = i__ - *ka;
	    i__4 = i__;
	    for (j = max(i__2,i__3); j <= i__4; ++j) {
		ab_ref(i__ - j + 1, j) = ab_ref(i__ - j + 1, j) / bii;
/* L260: */
	    }
	    i__4 = i__ - 1;
	    for (k = i__ - kbt; k <= i__4; ++k) {
		i__2 = k;
		for (j = i__ - kbt; j <= i__2; ++j) {
		    ab_ref(k - j + 1, j) = ab_ref(k - j + 1, j) - bb_ref(i__ 
			    - j + 1, j) * ab_ref(i__ - k + 1, k) - bb_ref(i__ 
			    - k + 1, k) * ab_ref(i__ - j + 1, j) + ab_ref(1, 
			    i__) * bb_ref(i__ - j + 1, j) * bb_ref(i__ - k + 
			    1, k);
/* L270: */
		}
/* Computing MAX */
		i__2 = 1, i__3 = i__ - *ka;
		i__1 = i__ - kbt - 1;
		for (j = max(i__2,i__3); j <= i__1; ++j) {
		    ab_ref(k - j + 1, j) = ab_ref(k - j + 1, j) - bb_ref(i__ 
			    - k + 1, k) * ab_ref(i__ - j + 1, j);
/* L280: */
		}
/* L290: */
	    }
	    i__4 = i1;
	    for (j = i__; j <= i__4; ++j) {
/* Computing MAX */
		i__1 = j - *ka, i__2 = i__ - kbt;
		i__3 = i__ - 1;
		for (k = max(i__1,i__2); k <= i__3; ++k) {
		    ab_ref(j - k + 1, k) = ab_ref(j - k + 1, k) - bb_ref(i__ 
			    - k + 1, k) * ab_ref(j - i__ + 1, i__);
/* L300: */
		}
/* L310: */
	    }

	    if (wantx) {

/*              post-multiply X by inv(S(i)) */

		i__4 = *n - m;
		d__1 = 1. / bii;
		dscal_(&i__4, &d__1, &x_ref(m + 1, i__), &c__1);
		if (kbt > 0) {
		    i__4 = *n - m;
		    i__3 = *ldbb - 1;
		    dger_(&i__4, &kbt, &c_b20, &x_ref(m + 1, i__), &c__1, &
			    bb_ref(kbt + 1, i__ - kbt), &i__3, &x_ref(m + 1, 
			    i__ - kbt), ldx);
		}
	    }

/*           store a(i1,i) in RA1 for use in next loop over K */

	    ra1 = ab_ref(i1 - i__ + 1, i__);
	}

/*        Generate and apply vectors of rotations to chase all the   
          existing bulges KA positions down toward the bottom of the   
          band */

	i__4 = *kb - 1;
	for (k = 1; k <= i__4; ++k) {
	    if (update) {

/*              Determine the rotations which would annihilate the bulge   
                which has in theory just been created */

		if (i__ - k + *ka < *n && i__ - k > 1) {

/*                 generate rotation to annihilate a(i-k+ka+1,i) */

		    dlartg_(&ab_ref(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_ref(k + 1, i__ - k) * ra1;
		    work[i__ - k] = work[*n + i__ - k + *ka - m] * t - work[
			    i__ - k + *ka - m] * ab_ref(ka1, i__ - k);
		    ab_ref(ka1, i__ - k) = work[i__ - k + *ka - m] * t + work[
			    *n + i__ - k + *ka - m] * ab_ref(ka1, i__ - k);
		    ra1 = ra;
		}
	    }
/* Computing MAX */
	    i__3 = 1, i__1 = k - i0 + 2;
	    j2 = i__ - k - 1 + max(i__3,i__1) * ka1;
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    if (update) {
/* Computing MAX */
		i__3 = j2, i__1 = i__ + (*ka << 1) - k + 1;
		j2t = max(i__3,i__1);
	    } else {
		j2t = j2;
	    }
	    nrt = (*n - j2t + *ka) / ka1;
	    i__3 = j1;
	    i__1 = ka1;
	    for (j = j2t; i__1 < 0 ? j >= i__3 : j <= i__3; j += i__1) {

/*              create nonzero element a(j+1,j-ka) outside the band   
                and store it in WORK(j-m) */

		work[j - m] *= ab_ref(ka1, j - *ka + 1);
		ab_ref(ka1, j - *ka + 1) = work[*n + j - m] * ab_ref(ka1, j - 
			*ka + 1);
/* L320: */
	    }

/*           generate rotations in 1st set to annihilate elements which   
             have been created outside the band */

	    if (nrt > 0) {
		dlargv_(&nrt, &ab_ref(ka1, j2t - *ka), &inca, &work[j2t - m], 
			&ka1, &work[*n + j2t - m], &ka1);
	    }
	    if (nr > 0) {

/*              apply rotations in 1st set from the left */

		i__1 = *ka - 1;
		for (l = 1; l <= i__1; ++l) {
		    dlartv_(&nr, &ab_ref(l + 1, j2 - l), &inca, &ab_ref(l + 2,
			     j2 - l), &inca, &work[*n + j2 - m], &work[j2 - m]
			    , &ka1);
/* L330: */
		}

/*              apply rotations in 1st set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(1, j2), &ab_ref(1, j2 + 1), &ab_ref(2, 
			j2), &inca, &work[*n + j2 - m], &work[j2 - m], &ka1);

	    }

/*           start applying rotations in 1st set from the right */

	    i__1 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__1; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j2), &inca, &ab_ref(
			    ka1 - l, j2 + 1), &inca, &work[*n + j2 - m], &
			    work[j2 - m], &ka1);
		}
/* L340: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 1st set */

		i__1 = j1;
		i__3 = ka1;
		for (j = j2; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
		    i__2 = *n - m;
		    drot_(&i__2, &x_ref(m + 1, j), &c__1, &x_ref(m + 1, j + 1)
			    , &c__1, &work[*n + j - m], &work[j - m]);
/* L350: */
		}
	    }
/* L360: */
	}

	if (update) {
	    if (i2 <= *n && kbt > 0) {

/*              create nonzero element a(i-kbt+ka+1,i-kbt) outside the   
                band and store it in WORK(i-kbt) */

		work[i__ - kbt] = -bb_ref(kbt + 1, i__ - kbt) * ra1;
	    }
	}

	for (k = *kb; k >= 1; --k) {
	    if (update) {
/* Computing MAX */
		i__4 = 2, i__3 = k - i0 + 1;
		j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
	    } else {
/* Computing MAX */
		i__4 = 1, i__3 = k - i0 + 1;
		j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
	    }

/*           finish applying rotations in 2nd set from the right */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (*n - j2 + *ka + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j2 - *ka), &inca, &
			    ab_ref(ka1 - l, j2 - *ka + 1), &inca, &work[*n + 
			    j2 - *ka], &work[j2 - *ka], &ka1);
		}
/* L370: */
	    }
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    i__4 = j2;
	    i__3 = -ka1;
	    for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
		work[j] = work[j - *ka];
		work[*n + j] = work[*n + j - *ka];
/* L380: */
	    }
	    i__3 = j1;
	    i__4 = ka1;
	    for (j = j2; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {

/*              create nonzero element a(j+1,j-ka) outside the band   
                and store it in WORK(j) */

		work[j] *= ab_ref(ka1, j - *ka + 1);
		ab_ref(ka1, j - *ka + 1) = work[*n + j] * ab_ref(ka1, j - *ka 
			+ 1);
/* L390: */
	    }
	    if (update) {
		if (i__ - k < *n - *ka && k <= kbt) {
		    work[i__ - k + *ka] = work[i__ - k];
		}
	    }
/* L400: */
	}

	for (k = *kb; k >= 1; --k) {
/* Computing MAX */
	    i__4 = 1, i__3 = k - i0 + 1;
	    j2 = i__ - k - 1 + max(i__4,i__3) * ka1;
	    nr = (*n - j2 + *ka) / ka1;
	    j1 = j2 + (nr - 1) * ka1;
	    if (nr > 0) {

/*              generate rotations in 2nd set to annihilate elements   
                which have been created outside the band */

		dlargv_(&nr, &ab_ref(ka1, j2 - *ka), &inca, &work[j2], &ka1, &
			work[*n + j2], &ka1);

/*              apply rotations in 2nd set from the left */

		i__4 = *ka - 1;
		for (l = 1; l <= i__4; ++l) {
		    dlartv_(&nr, &ab_ref(l + 1, j2 - l), &inca, &ab_ref(l + 2,
			     j2 - l), &inca, &work[*n + j2], &work[j2], &ka1);
/* L410: */
		}

/*              apply rotations in 2nd set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(1, j2), &ab_ref(1, j2 + 1), &ab_ref(2, 
			j2), &inca, &work[*n + j2], &work[j2], &ka1);

	    }

/*           start applying rotations in 2nd set from the right */

	    i__4 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__4; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j2), &inca, &ab_ref(
			    ka1 - l, j2 + 1), &inca, &work[*n + j2], &work[j2]
			    , &ka1);
		}
/* L420: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 2nd set */

		i__4 = j1;
		i__3 = ka1;
		for (j = j2; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
		    i__1 = *n - m;
		    drot_(&i__1, &x_ref(m + 1, j), &c__1, &x_ref(m + 1, j + 1)
			    , &c__1, &work[*n + j], &work[j]);
/* L430: */
		}
	    }
/* L440: */
	}

	i__3 = *kb - 1;
	for (k = 1; k <= i__3; ++k) {
/* Computing MAX */
	    i__4 = 1, i__1 = k - i0 + 2;
	    j2 = i__ - k - 1 + max(i__4,i__1) * ka1;

/*           finish applying rotations in 1st set from the right */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (*n - j2 + l) / ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j2), &inca, &ab_ref(
			    ka1 - l, j2 + 1), &inca, &work[*n + j2 - m], &
			    work[j2 - m], &ka1);
		}
/* L450: */
	    }
/* L460: */
	}

	if (*kb > 1) {
	    i__3 = i__ - *kb + (*ka << 1) + 1;
	    for (j = *n - 1; j >= i__3; --j) {
		work[*n + j - m] = work[*n + j - *ka - m];
		work[j - m] = work[j - *ka - m];
/* L470: */
	    }
	}

    }

    goto L10;

L480:

/*     **************************** Phase 2 *****************************   

       The logical structure of this phase is:   

       UPDATE = .TRUE.   
       DO I = 1, M   
          use S(i) to update A and create a new bulge   
          apply rotations to push all bulges KA positions upward   
       END DO   
       UPDATE = .FALSE.   
       DO I = M - KA - 1, 2, -1   
          apply rotations to push all bulges KA positions upward   
       END DO   

       To avoid duplicating code, the two loops are merged. */

    update = TRUE_;
    i__ = 0;
L490:
    if (update) {
	++i__;
/* Computing MIN */
	i__3 = *kb, i__4 = m - i__;
	kbt = min(i__3,i__4);
	i0 = i__ + 1;
/* Computing MAX */
	i__3 = 1, i__4 = i__ - *ka;
	i1 = max(i__3,i__4);
	i2 = i__ + kbt - ka1;
	if (i__ > m) {
	    update = FALSE_;
	    --i__;
	    i0 = m + 1;
	    if (*ka == 0) {
		return 0;
	    }
	    goto L490;
	}
    } else {
	i__ -= *ka;
	if (i__ < 2) {
	    return 0;
	}
    }

    if (i__ < m - kbt) {
	nx = m;
    } else {
	nx = *n;
    }

    if (upper) {

/*        Transform A, working with the upper triangle */

	if (update) {

/*           Form  inv(S(i))**T * A * inv(S(i)) */

	    bii = bb_ref(kb1, i__);
	    i__3 = i__;
	    for (j = i1; j <= i__3; ++j) {
		ab_ref(j - i__ + ka1, i__) = ab_ref(j - i__ + ka1, i__) / bii;
/* L500: */
	    }
/* Computing MIN */
	    i__4 = *n, i__1 = i__ + *ka;
	    i__3 = min(i__4,i__1);
	    for (j = i__; j <= i__3; ++j) {
		ab_ref(i__ - j + ka1, j) = ab_ref(i__ - j + ka1, j) / bii;
/* L510: */
	    }
	    i__3 = i__ + kbt;
	    for (k = i__ + 1; k <= i__3; ++k) {
		i__4 = i__ + kbt;
		for (j = k; j <= i__4; ++j) {
		    ab_ref(k - j + ka1, j) = ab_ref(k - j + ka1, j) - bb_ref(
			    i__ - j + kb1, j) * ab_ref(i__ - k + ka1, k) - 
			    bb_ref(i__ - k + kb1, k) * ab_ref(i__ - j + ka1, 
			    j) + ab_ref(ka1, i__) * bb_ref(i__ - j + kb1, j) *
			     bb_ref(i__ - k + kb1, k);
/* L520: */
		}
/* Computing MIN */
		i__1 = *n, i__2 = i__ + *ka;
		i__4 = min(i__1,i__2);
		for (j = i__ + kbt + 1; j <= i__4; ++j) {
		    ab_ref(k - j + ka1, j) = ab_ref(k - j + ka1, j) - bb_ref(
			    i__ - k + kb1, k) * ab_ref(i__ - j + ka1, j);
/* L530: */
		}
/* L540: */
	    }
	    i__3 = i__;
	    for (j = i1; j <= i__3; ++j) {
/* Computing MIN */
		i__1 = j + *ka, i__2 = i__ + kbt;
		i__4 = min(i__1,i__2);
		for (k = i__ + 1; k <= i__4; ++k) {
		    ab_ref(j - k + ka1, k) = ab_ref(j - k + ka1, k) - bb_ref(
			    i__ - k + kb1, k) * ab_ref(j - i__ + ka1, i__);
/* L550: */
		}
/* L560: */
	    }

	    if (wantx) {

/*              post-multiply X by inv(S(i)) */

		d__1 = 1. / bii;
		dscal_(&nx, &d__1, &x_ref(1, i__), &c__1);
		if (kbt > 0) {
		    i__3 = *ldbb - 1;
		    dger_(&nx, &kbt, &c_b20, &x_ref(1, i__), &c__1, &bb_ref(*
			    kb, i__ + 1), &i__3, &x_ref(1, i__ + 1), ldx);
		}
	    }

/*           store a(i1,i) in RA1 for use in next loop over K */

	    ra1 = ab_ref(i1 - i__ + ka1, i__);
	}

/*        Generate and apply vectors of rotations to chase all the   
          existing bulges KA positions up toward the top of the band */

	i__3 = *kb - 1;
	for (k = 1; k <= i__3; ++k) {
	    if (update) {

/*              Determine the rotations which would annihilate the bulge   
                which has in theory just been created */

		if (i__ + k - ka1 > 0 && i__ + k < m) {

/*                 generate rotation to annihilate a(i+k-ka-1,i) */

		    dlartg_(&ab_ref(k + 1, i__), &ra1, &work[*n + i__ + k - *
			    ka], &work[i__ + k - *ka], &ra);

/*                 create nonzero element a(i+k-ka-1,i+k) outside the   
                   band and store it in WORK(m-kb+i+k) */

		    t = -bb_ref(kb1 - k, i__ + k) * ra1;
		    work[m - *kb + i__ + k] = work[*n + i__ + k - *ka] * t - 
			    work[i__ + k - *ka] * ab_ref(1, i__ + k);
		    ab_ref(1, i__ + k) = work[i__ + k - *ka] * t + work[*n + 
			    i__ + k - *ka] * ab_ref(1, i__ + k);
		    ra1 = ra;
		}
	    }
/* Computing MAX */
	    i__4 = 1, i__1 = k + i0 - m + 1;
	    j2 = i__ + k + 1 - max(i__4,i__1) * ka1;
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    if (update) {
/* Computing MIN */
		i__4 = j2, i__1 = i__ - (*ka << 1) + k - 1;
		j2t = min(i__4,i__1);
	    } else {
		j2t = j2;
	    }
	    nrt = (j2t + *ka - 1) / ka1;
	    i__4 = j2t;
	    i__1 = ka1;
	    for (j = j1; i__1 < 0 ? j >= i__4 : j <= i__4; j += i__1) {

/*              create nonzero element a(j-1,j+ka) outside the band   
                and store it in WORK(j) */

		work[j] *= ab_ref(1, j + *ka - 1);
		ab_ref(1, j + *ka - 1) = work[*n + j] * ab_ref(1, j + *ka - 1)
			;
/* L570: */
	    }

/*           generate rotations in 1st set to annihilate elements which   
             have been created outside the band */

	    if (nrt > 0) {
		dlargv_(&nrt, &ab_ref(1, j1 + *ka), &inca, &work[j1], &ka1, &
			work[*n + j1], &ka1);
	    }
	    if (nr > 0) {

/*              apply rotations in 1st set from the left */

		i__1 = *ka - 1;
		for (l = 1; l <= i__1; ++l) {
		    dlartv_(&nr, &ab_ref(ka1 - l, j1 + l), &inca, &ab_ref(*ka 
			    - l, j1 + l), &inca, &work[*n + j1], &work[j1], &
			    ka1);
/* L580: */
		}

/*              apply rotations in 1st set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(ka1, j1), &ab_ref(ka1, j1 - 1), &ab_ref(*
			ka, j1), &inca, &work[*n + j1], &work[j1], &ka1);

	    }

/*           start applying rotations in 1st set from the right */

	    i__1 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__1; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j1t), &inca, &ab_ref(l + 1, j1t 
			    - 1), &inca, &work[*n + j1t], &work[j1t], &ka1);
		}
/* L590: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 1st set */

		i__1 = j2;
		i__4 = ka1;
		for (j = j1; i__4 < 0 ? j >= i__1 : j <= i__1; j += i__4) {
		    drot_(&nx, &x_ref(1, j), &c__1, &x_ref(1, j - 1), &c__1, &
			    work[*n + j], &work[j]);
/* L600: */
		}
	    }
/* L610: */
	}

	if (update) {
	    if (i2 > 0 && kbt > 0) {

/*              create nonzero element a(i+kbt-ka-1,i+kbt) outside the   
                band and store it in WORK(m-kb+i+kbt) */

		work[m - *kb + i__ + kbt] = -bb_ref(kb1 - kbt, i__ + kbt) * 
			ra1;
	    }
	}

	for (k = *kb; k >= 1; --k) {
	    if (update) {
/* Computing MAX */
		i__3 = 2, i__4 = k + i0 - m;
		j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
	    } else {
/* Computing MAX */
		i__3 = 1, i__4 = k + i0 - m;
		j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
	    }

/*           finish applying rotations in 2nd set from the right */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (j2 + *ka + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j1t + *ka), &inca, &ab_ref(l + 1,
			     j1t + *ka - 1), &inca, &work[*n + m - *kb + j1t 
			    + *ka], &work[m - *kb + j1t + *ka], &ka1);
		}
/* L620: */
	    }
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    i__3 = j2;
	    i__4 = ka1;
	    for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
		work[m - *kb + j] = work[m - *kb + j + *ka];
		work[*n + m - *kb + j] = work[*n + m - *kb + j + *ka];
/* L630: */
	    }
	    i__4 = j2;
	    i__3 = ka1;
	    for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {

/*              create nonzero element a(j-1,j+ka) outside the band   
                and store it in WORK(m-kb+j) */

		work[m - *kb + j] *= ab_ref(1, j + *ka - 1);
		ab_ref(1, j + *ka - 1) = work[*n + m - *kb + j] * ab_ref(1, j 
			+ *ka - 1);
/* L640: */
	    }
	    if (update) {
		if (i__ + k > ka1 && k <= kbt) {
		    work[m - *kb + i__ + k - *ka] = work[m - *kb + i__ + k];
		}
	    }
/* L650: */
	}

	for (k = *kb; k >= 1; --k) {
/* Computing MAX */
	    i__3 = 1, i__4 = k + i0 - m;
	    j2 = i__ + k + 1 - max(i__3,i__4) * ka1;
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    if (nr > 0) {

/*              generate rotations in 2nd set to annihilate elements   
                which have been created outside the band */

		dlargv_(&nr, &ab_ref(1, j1 + *ka), &inca, &work[m - *kb + j1],
			 &ka1, &work[*n + m - *kb + j1], &ka1);

/*              apply rotations in 2nd set from the left */

		i__3 = *ka - 1;
		for (l = 1; l <= i__3; ++l) {
		    dlartv_(&nr, &ab_ref(ka1 - l, j1 + l), &inca, &ab_ref(*ka 
			    - l, j1 + l), &inca, &work[*n + m - *kb + j1], &
			    work[m - *kb + j1], &ka1);
/* L660: */
		}

/*              apply rotations in 2nd set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(ka1, j1), &ab_ref(ka1, j1 - 1), &ab_ref(*
			ka, j1), &inca, &work[*n + m - *kb + j1], &work[m - *
			kb + j1], &ka1);

	    }

/*           start applying rotations in 2nd set from the right */

	    i__3 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__3; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j1t), &inca, &ab_ref(l + 1, j1t 
			    - 1), &inca, &work[*n + m - *kb + j1t], &work[m - 
			    *kb + j1t], &ka1);
		}
/* L670: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 2nd set */

		i__3 = j2;
		i__4 = ka1;
		for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
		    drot_(&nx, &x_ref(1, j), &c__1, &x_ref(1, j - 1), &c__1, &
			    work[*n + m - *kb + j], &work[m - *kb + j]);
/* L680: */
		}
	    }
/* L690: */
	}

	i__4 = *kb - 1;
	for (k = 1; k <= i__4; ++k) {
/* Computing MAX */
	    i__3 = 1, i__1 = k + i0 - m + 1;
	    j2 = i__ + k + 1 - max(i__3,i__1) * ka1;

/*           finish applying rotations in 1st set from the right */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(l, j1t), &inca, &ab_ref(l + 1, j1t 
			    - 1), &inca, &work[*n + j1t], &work[j1t], &ka1);
		}
/* L700: */
	    }
/* L710: */
	}

	if (*kb > 1) {
/* Computing MIN */
	    i__3 = i__ + *kb;
	    i__4 = min(i__3,m) - (*ka << 1) - 1;
	    for (j = 2; j <= i__4; ++j) {
		work[*n + j] = work[*n + j + *ka];
		work[j] = work[j + *ka];
/* L720: */
	    }
	}

    } else {

/*        Transform A, working with the lower triangle */

	if (update) {

/*           Form  inv(S(i))**T * A * inv(S(i)) */

	    bii = bb_ref(1, i__);
	    i__4 = i__;
	    for (j = i1; j <= i__4; ++j) {
		ab_ref(i__ - j + 1, j) = ab_ref(i__ - j + 1, j) / bii;
/* L730: */
	    }
/* Computing MIN */
	    i__3 = *n, i__1 = i__ + *ka;
	    i__4 = min(i__3,i__1);
	    for (j = i__; j <= i__4; ++j) {
		ab_ref(j - i__ + 1, i__) = ab_ref(j - i__ + 1, i__) / bii;
/* L740: */
	    }
	    i__4 = i__ + kbt;
	    for (k = i__ + 1; k <= i__4; ++k) {
		i__3 = i__ + kbt;
		for (j = k; j <= i__3; ++j) {
		    ab_ref(j - k + 1, k) = ab_ref(j - k + 1, k) - bb_ref(j - 
			    i__ + 1, i__) * ab_ref(k - i__ + 1, i__) - bb_ref(
			    k - i__ + 1, i__) * ab_ref(j - i__ + 1, i__) + 
			    ab_ref(1, i__) * bb_ref(j - i__ + 1, i__) * 
			    bb_ref(k - i__ + 1, i__);
/* L750: */
		}
/* Computing MIN */
		i__1 = *n, i__2 = i__ + *ka;
		i__3 = min(i__1,i__2);
		for (j = i__ + kbt + 1; j <= i__3; ++j) {
		    ab_ref(j - k + 1, k) = ab_ref(j - k + 1, k) - bb_ref(k - 
			    i__ + 1, i__) * ab_ref(j - i__ + 1, i__);
/* L760: */
		}
/* L770: */
	    }
	    i__4 = i__;
	    for (j = i1; j <= i__4; ++j) {
/* Computing MIN */
		i__1 = j + *ka, i__2 = i__ + kbt;
		i__3 = min(i__1,i__2);
		for (k = i__ + 1; k <= i__3; ++k) {
		    ab_ref(k - j + 1, j) = ab_ref(k - j + 1, j) - bb_ref(k - 
			    i__ + 1, i__) * ab_ref(i__ - j + 1, j);
/* L780: */
		}
/* L790: */
	    }

	    if (wantx) {

/*              post-multiply X by inv(S(i)) */

		d__1 = 1. / bii;
		dscal_(&nx, &d__1, &x_ref(1, i__), &c__1);
		if (kbt > 0) {
		    dger_(&nx, &kbt, &c_b20, &x_ref(1, i__), &c__1, &bb_ref(2,
			     i__), &c__1, &x_ref(1, i__ + 1), ldx);
		}
	    }

/*           store a(i,i1) in RA1 for use in next loop over K */

	    ra1 = ab_ref(i__ - i1 + 1, i1);
	}

/*        Generate and apply vectors of rotations to chase all the   
          existing bulges KA positions up toward the top of the band */

	i__4 = *kb - 1;
	for (k = 1; k <= i__4; ++k) {
	    if (update) {

/*              Determine the rotations which would annihilate the bulge   
                which has in theory just been created */

		if (i__ + k - ka1 > 0 && i__ + k < m) {

/*                 generate rotation to annihilate a(i,i+k-ka-1) */

		    dlartg_(&ab_ref(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_ref(k + 1, i__) * ra1;
		    work[m - *kb + i__ + k] = work[*n + i__ + k - *ka] * t - 
			    work[i__ + k - *ka] * ab_ref(ka1, i__ + k - *ka);
		    ab_ref(ka1, i__ + k - *ka) = work[i__ + k - *ka] * t + 
			    work[*n + i__ + k - *ka] * ab_ref(ka1, i__ + k - *
			    ka);
		    ra1 = ra;
		}
	    }
/* Computing MAX */
	    i__3 = 1, i__1 = k + i0 - m + 1;
	    j2 = i__ + k + 1 - max(i__3,i__1) * ka1;
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    if (update) {
/* Computing MIN */
		i__3 = j2, i__1 = i__ - (*ka << 1) + k - 1;
		j2t = min(i__3,i__1);
	    } else {
		j2t = j2;
	    }
	    nrt = (j2t + *ka - 1) / ka1;
	    i__3 = j2t;
	    i__1 = ka1;
	    for (j = j1; i__1 < 0 ? j >= i__3 : j <= i__3; j += i__1) {

/*              create nonzero element a(j+ka,j-1) outside the band   
                and store it in WORK(j) */

		work[j] *= ab_ref(ka1, j - 1);
		ab_ref(ka1, j - 1) = work[*n + j] * ab_ref(ka1, j - 1);
/* L800: */
	    }

/*           generate rotations in 1st set to annihilate elements which   
             have been created outside the band */

	    if (nrt > 0) {
		dlargv_(&nrt, &ab_ref(ka1, j1), &inca, &work[j1], &ka1, &work[
			*n + j1], &ka1);
	    }
	    if (nr > 0) {

/*              apply rotations in 1st set from the right */

		i__1 = *ka - 1;
		for (l = 1; l <= i__1; ++l) {
		    dlartv_(&nr, &ab_ref(l + 1, j1), &inca, &ab_ref(l + 2, j1 
			    - 1), &inca, &work[*n + j1], &work[j1], &ka1);
/* L810: */
		}

/*              apply rotations in 1st set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(1, j1), &ab_ref(1, j1 - 1), &ab_ref(2, 
			j1 - 1), &inca, &work[*n + j1], &work[j1], &ka1);

	    }

/*           start applying rotations in 1st set from the left */

	    i__1 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__1; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j1t - ka1 + l), &inca, 
			    &ab_ref(ka1 - l, j1t - ka1 + l), &inca, &work[*n 
			    + j1t], &work[j1t], &ka1);
		}
/* L820: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 1st set */

		i__1 = j2;
		i__3 = ka1;
		for (j = j1; i__3 < 0 ? j >= i__1 : j <= i__1; j += i__3) {
		    drot_(&nx, &x_ref(1, j), &c__1, &x_ref(1, j - 1), &c__1, &
			    work[*n + j], &work[j]);
/* L830: */
		}
	    }
/* L840: */
	}

	if (update) {
	    if (i2 > 0 && kbt > 0) {

/*              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_ref(kbt + 1, i__) * ra1;
	    }
	}

	for (k = *kb; k >= 1; --k) {
	    if (update) {
/* Computing MAX */
		i__4 = 2, i__3 = k + i0 - m;
		j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
	    } else {
/* Computing MAX */
		i__4 = 1, i__3 = k + i0 - m;
		j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
	    }

/*           finish applying rotations in 2nd set from the left */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (j2 + *ka + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j1t + l - 1), &inca, &
			    ab_ref(ka1 - l, j1t + l - 1), &inca, &work[*n + m 
			    - *kb + j1t + *ka], &work[m - *kb + j1t + *ka], &
			    ka1);
		}
/* L850: */
	    }
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    i__4 = j2;
	    i__3 = ka1;
	    for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
		work[m - *kb + j] = work[m - *kb + j + *ka];
		work[*n + m - *kb + j] = work[*n + m - *kb + j + *ka];
/* L860: */
	    }
	    i__3 = j2;
	    i__4 = ka1;
	    for (j = j1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {

/*              create nonzero element a(j+ka,j-1) outside the band   
                and store it in WORK(m-kb+j) */

		work[m - *kb + j] *= ab_ref(ka1, j - 1);
		ab_ref(ka1, j - 1) = work[*n + m - *kb + j] * ab_ref(ka1, j - 
			1);
/* L870: */
	    }
	    if (update) {
		if (i__ + k > ka1 && k <= kbt) {
		    work[m - *kb + i__ + k - *ka] = work[m - *kb + i__ + k];
		}
	    }
/* L880: */
	}

	for (k = *kb; k >= 1; --k) {
/* Computing MAX */
	    i__4 = 1, i__3 = k + i0 - m;
	    j2 = i__ + k + 1 - max(i__4,i__3) * ka1;
	    nr = (j2 + *ka - 1) / ka1;
	    j1 = j2 - (nr - 1) * ka1;
	    if (nr > 0) {

/*              generate rotations in 2nd set to annihilate elements   
                which have been created outside the band */

		dlargv_(&nr, &ab_ref(ka1, j1), &inca, &work[m - *kb + j1], &
			ka1, &work[*n + m - *kb + j1], &ka1);

/*              apply rotations in 2nd set from the right */

		i__4 = *ka - 1;
		for (l = 1; l <= i__4; ++l) {
		    dlartv_(&nr, &ab_ref(l + 1, j1), &inca, &ab_ref(l + 2, j1 
			    - 1), &inca, &work[*n + m - *kb + j1], &work[m - *
			    kb + j1], &ka1);
/* L890: */
		}

/*              apply rotations in 2nd set from both sides to diagonal   
                blocks */

		dlar2v_(&nr, &ab_ref(1, j1), &ab_ref(1, j1 - 1), &ab_ref(2, 
			j1 - 1), &inca, &work[*n + m - *kb + j1], &work[m - *
			kb + j1], &ka1);

	    }

/*           start applying rotations in 2nd set from the left */

	    i__4 = *kb - k + 1;
	    for (l = *ka - 1; l >= i__4; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j1t - ka1 + l), &inca, 
			    &ab_ref(ka1 - l, j1t - ka1 + l), &inca, &work[*n 
			    + m - *kb + j1t], &work[m - *kb + j1t], &ka1);
		}
/* L900: */
	    }

	    if (wantx) {

/*              post-multiply X by product of rotations in 2nd set */

		i__4 = j2;
		i__3 = ka1;
		for (j = j1; i__3 < 0 ? j >= i__4 : j <= i__4; j += i__3) {
		    drot_(&nx, &x_ref(1, j), &c__1, &x_ref(1, j - 1), &c__1, &
			    work[*n + m - *kb + j], &work[m - *kb + j]);
/* L910: */
		}
	    }
/* L920: */
	}

	i__3 = *kb - 1;
	for (k = 1; k <= i__3; ++k) {
/* Computing MAX */
	    i__4 = 1, i__1 = k + i0 - m + 1;
	    j2 = i__ + k + 1 - max(i__4,i__1) * ka1;

/*           finish applying rotations in 1st set from the left */

	    for (l = *kb - k; l >= 1; --l) {
		nrt = (j2 + l - 1) / ka1;
		j1t = j2 - (nrt - 1) * ka1;
		if (nrt > 0) {
		    dlartv_(&nrt, &ab_ref(ka1 - l + 1, j1t - ka1 + l), &inca, 
			    &ab_ref(ka1 - l, j1t - ka1 + l), &inca, &work[*n 
			    + j1t], &work[j1t], &ka1);
		}
/* L930: */
	    }
/* L940: */
	}

	if (*kb > 1) {
/* Computing MIN */
	    i__4 = i__ + *kb;
	    i__3 = min(i__4,m) - (*ka << 1) - 1;
	    for (j = 2; j <= i__3; ++j) {
		work[*n + j] = work[*n + j + *ka];
		work[j] = work[j + *ka];
/* L950: */
	    }
	}

    }

    goto L490;

/*     End of DSBGST */

} /* dsbgst_ */

#undef bb_ref
#undef ab_ref
#undef x_ref


