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 