SUBROUTINE QLSTEP(U, LDU, V, LDV, Q, E, M, N, I, K, SHIFT,
     *                  WANTU, WANTV)
C
C     PURPOSE:
C
C     The subroutine QLSTEP performs one QL iteration step onto the
C     unreduced subbidiagonal Jk:
C
C              !Q(i) E(i+1)  0  ...    0  !
C              ! 0   Q(i+1) E(i+2)     .  !
C         Jk = ! .                     .  !
C              ! .                        !
C              ! .                    E(k)!
C              ! 0   ...              Q(k)!
C
C     with k <= p and i >= 1, p = min(M,N), of the bidiagonal J:
C
C              !Q(1) E(2)  0    ...   0  !
C              ! 0   Q(2) E(3)        .  !
C          J = ! .                    .  !
C              ! .                   E(p)!
C              ! 0   ...             Q(p)!
C
C     Hereby, Jk is transformed to  S'Jk T with S and T products of
C     Givens rotations. These Givens rotations S (resp.,T) will be post-
C     multiplied into U (resp.,V), if WANTU (resp.,WANTV) = .TRUE.
C
C     ARGUMENT LIST:
C
C     U - DOUBLE PRECISION array of DIMENSION (LDU,min(M,N)).
C         On entry, U may contain the M by p (p=min(M,N)) left transfor-
C         mation matrix.
C         On return, if WANTU = .TRUE., the Givens rotations S on the
C         left have been postmultiplied into U.
C         NOTE: U is not referenced if WANTU = .FALSE.
C     LDU - INTEGER.
C         LDU is the leading dimension of the array U (LDU >= M).
C     V - DOUBLE PRECISION array of DIMENSION (LDV,min(M,N)).
C         On entry, V may contain the N by p (p=min(M,N)) right trans-
C         formation matrix.
C         On return, if WANTV = .TRUE., the Givens rotations T on the
C         right have been postmultiplied into V.
C         NOTE: V is not referenced if WANTV is .false.
C     LDV - INTEGER.
C         LDV is the leading dimension of the array V (LDV >= N).
C     Q - DOUBLE PRECISION array of DIMENSION (min(M,N)).
C         On entry, Q contains the diagonal entries of the bidiagonal J.
C         On return, Q contains the diagonal entries of the transformed
C         matrix S' J T.
C     E - DOUBLE PRECISION array of DIMENSION (min(M,N)).
C         On entry, E contains the superdiagonal entries of J.
C         On return, E contains the superdiagonal entries of the trans-
C         formed matrix S' J T. E(k+1) = 0 if k < min(M,N).
C     M - INTEGER.
C         M is the number of rows of the matrix U.
C     N - INTEGER.
C         N is the number of rows of the matrix V.
C     I - INTEGER.
C         I is the index of the first diagonal entry of the considered
C         unreduced subbidiagonal Jk of J.
C     K - INTEGER.
C         K is the index of the last diagonal entry of the considered
C         unreduced subbidiagonal Jk of J.
C     SHIFT - DOUBLE PRECISION.
C         Value of the shift used in the QL iteration step.
C     WANTU - LOGICAL.
C         WANTU = .TRUE. if the Givens rotations S must be postmulti-
C         plied on the left into U, else .FALSE.
C     WANTV - LOGICAL.
C         WANTV = .TRUE. if the Givens rotations T must be postmulti-
C         plied  on the left into V, else .FALSE.
C
C     EXTERNAL SUBROUTINES AND FUNCTIONS:
C
C     DROT from BLAS.
C
C     METHOD DESCRIPTION:
C
C     QL iterations diagonalize the bidiagonal by zeroing the super-
C     diagonal elements of Jk from top to bottom.
C     The routine QLSTEP overwrites Jk with the bidiagonal matrix
C     S' Jk T where S and T are Givens rotations.
C     T is essentially the orthogonal matrix that would be obtained by
C     applying one implicit symmetric shift QL step onto the matrix
C     Jk'Jk. This step factors the matrix (Jk'Jk - shift*I) into a
C     product of an orthogonal matrix T and a lower triangular matrix.
C     See [1,Sec.8.2-8.3] and [2] for more details.
C
C     REFERENCES:
C     [1] G.H. Golub and C.F. Van Loan, Matrix Computations. The Johns
C         Hopkins University Press, Baltimore,Maryland (1983).
C     [2] H. Bowdler, R.S. Martin and J.H. Wilkinson, The QR and QL
C         algorithms for symmetric matrices. Numer. Math., 11 (1968),
C         293-306.
C
C     CONTRIBUTOR: S. Van Huffel (ESAT Laboratory, KU Leuven).
C
C     REVISIONS: 1988, February 15.
C
C     .. Scalar Arguments ..
      INTEGER LDU, LDV, M, N, I, K
      DOUBLE PRECISION SHIFT
      LOGICAL WANTU, WANTV
C     .. Array Arguments ..
      DOUBLE PRECISION U(LDU,*), V(LDV,*), Q(*), E(*)
C     .. External Subroutines/Functions ..
      EXTERNAL DROT
C     .. Intrinsic Functions ..
      INTRINSIC MIN, SQRT
C     .. Local Scalars ..
      INTEGER IK, J, LL, JJ
      DOUBLE PRECISION F, G, H, C, S, X, Y, Z
C     .. Executable Statements ..
C
      X = Q(K)
      G = SHIFT
      F = (X - G) * (X + G)/X
      C = 1.0D0
      S = 1.0D0
      IK = K - I
      DO 10 JJ = 1, IK
         J = K - JJ
         LL= J + 1
         G = E(LL)
         Y = Q(J)
         H = S * G
         G = C * G
         Z = SQRT(F**2 + H**2)
         IF (JJ. GT. 1) E(LL+1) = Z
         C = F/Z
         S = H/Z
         F = X * C + G * S
         G = -X * S + G * C
         H = Y * S
         Y = Y * C
         IF (WANTU) CALL DROT(M, U(1,LL), 1, U(1,J), 1, C, S)
         Z = SQRT(F**2 + H**2)
         Q(LL) = Z
         C = F/Z
         S = H/Z
         F = C * G + S * Y
         X = -S * G + C * Y
         IF (WANTV) CALL DROT(N, V(1,LL), 1, V(1,J), 1, C, S)
   10 CONTINUE
      E(I+1) = F
      Q(I) = X
      IF (K .LT. MIN(M,N)) E(K+1) = 0.0D0
      RETURN
C Last line of QLSTEP *************************************************
      END