SUBROUTINE : QRQL 1 PURPOSE: The subroutine QRQL diagonalizes the given bidiagonal !Q(1) E(2) 0 ... 0 ! ! 0 Q(2) E(3) . ! J = ! . . ! (1) ! . E(p)! ! 0 ... Q(p)! with p = min(M,N) partially, using QR or QL iterations and in such a way that J is splitted into unreduced subbidiagonals Jj whose singular values are either all larger than a given bound or are all smaller than or equal to this bound. The Givens rotations performed on J on the left and right and corres- ponding to each QR or QL iteration step, are postmultiplied into the arrays U and V, if desired. 2 SPECIFICATION: SUBROUTINE QRQL(U, LDU, V, LDV, M, N, RANK, THETA, Q, E, INUL, TOL1, TOL2, WANTU, WANTV, IERR, IWARN) INTEGER LDU, LDV, M, N, RANK, IERR, IWARN DOUBLE PRECISION THETA, TOL1, TOL2 DOUBLE PRECISION U(LDU,min(M,N)), V(LDV,min(M,N)), Q(min(M,N)), E(min(M,N)) LOGICAL WANTU, WANTV LOGICAL INUL(min(M,N)) 3 ARGUMENT LIST: 3.1 ARGUMENTS IN U - DOUBLE PRECISION array of DIMENSION (LDU,min(M,N)). The leading M x p part (with p=min(M,N)) of this array con- tains either the leading M x p part of the identity matrix or a left transformation matrix applied to the original matrix of the problem. NOTE that U is changed by the routine if WANTU = .TRUE. and that U is not referenced if WANTU = .FALSE. LDU - INTEGER. LDU is the leading dimension of the array U (LDU >= M). V - DOUBLE PRECISION array of DIMENSION (LDV,min(M,N)). The leading N x p part (with p=min(M,N)) of this array con- tains either the leading N x p part of the identity matrix or a right transformation matrix applied to the original matrix of the problem. NOTE that V is changed by the routine if WANTV = .TRUE. and that V is not referenced if WANTV = .FALSE. LDV - INTEGER. LDV is the leading dimension of the array V (LDV >= N). M - INTEGER. M is the number of rows of the matrix U. N - INTEGER. N is the number of rows of the matrix V. RANK - INTEGER. On entry, there are 2 possibilities : i) RANK < 0: the rank is to be computed by the routine. ii) RANK >= 0: specifies the desired rank of J. If RANK > min(M,N), then the routine call is an empty statement. NOTE that RANK may be overwritten if J has multiple singular values. THETA - DOUBLE PRECISION. On entry, there are 2 possibilities depending on RANK. i) RANK < 0: THETA specifies an upper bound on the smallest singular values of J. THETA >= 0. ii) RANK >= 0: THETA is an initial estimate for computing an upper bound such that precisely RANK singular values of J are > this bound. If not available, assign a negative value (< 0) to THETA. NOTE that THETA is overwritten. Q - DOUBLE PRECISION array of DIMENSION (min(M,N)). Contains the diagonal entries of the bidiagonal J. NOTE that this array is overwritten. E - DOUBLE PRECISION array of DIMENSION (min(M,N)). E(i), i=2,...,p, with p=min(M,N), contain the superdiagonal entries of the bidiagonal J. E(1) = 0. NOTE that this array is overwritten. INUL - LOGICAL array of DIMENSION (min(M,N)). All elements of INUL must be initialized: INUL(i) = .TRUE. if the i-th column of U and/or V contains already a computed base vector of the desi- red singular subspace of the original matrix All other elements of INUL must be set to .FALSE. NOTE that this array is overwritten. 3.2 ARGUMENTS OUT U - DOUBLE PRECISION array of DIMENSION (LDU,min(M,N)). If WANTU = .TRUE., the Givens rotations on the left, corres- ponding to each QR or QL iteration step performed, have been postmultiplied into U. If WANTU = .FALSE., then U is not referenced. V - DOUBLE PRECISION array of DIMENSION (LDU,min(M,N)). If WANTV = .TRUE., the Givens rotations on the right, corres- ponding to each QR or QL iteration step performed, have been postmultiplied into V. If WANTV = .FALSE., then V is not referenced. RANK - INTEGER. If not specified by the user, then RANK is computed by the routine as the number of singular values > THETA. If specified by the user, then the specified RANK is changed by the routine if the RANK-th and the (RANK+1)-th singular value of A are considered to be equal. THETA - DOUBLE PRECISION. If RANK >= 0, then THETA specifies the computed upper bound such that precisely RANK singular values of A are > THETA + TOL1. Q - DOUBLE PRECISION array of DIMENSION (min(M,N)). Contains the diagonal entries of the transformed bidiagonal matrix J. E - DOUBLE PRECISION array of DIMENSION (min(M,N)). E(i), i=2,...,p, contain the superdiagonal entries of the transformed bidiagonal J. E(1) is unchanged. INUL - LOGICAL array of DIMENSION (min(M,N)). The indices of the elements of INUL with value .TRUE. indi- cate the indices of the diagonal entries of J which belong to those subbidiagonals whose singular values are all <= THETA. 3.4 TOLERANCES TOL1 - DOUBLE PRECISION. This parameter defines the multiplicity of singular values by considering all singular values within an interval of length TOL1 as coinciding. TOL1 is used in checking how many singu- lar values are <= THETA. Also in computing an appropriate upper bound THETA by a bisection method, TOL1 is used as stop criterion defining the minimal subinterval length. According to the numerical properties of the SVD, TOL1 must be >= !!J!! x EPS where !!J!! denotes the L2-norm and EPS is the machine precision. TOL2 - DOUBLE PRECISION. Working precision for the computations. This parameter specifies that elements Q(i) and E(i), which are <= TOL2 in absolute value, are considered to be zero. 3.5 MODE PARAMETERS WANTU - LOGICAL. If WANTU = .TRUE., the Givens rotations are postmultiplied on the left into U. WANTV - LOGICAL. If WANTV = .TRUE., the Givens rotations are postmultiplied on the right into V. 3.6 ERROR INDICATORS IERR - INTEGER. On return, IERR contains 0 unless the routine has failed. IWARN - INTEGER. On return, IWARN contains 0 unless RANK has been lowered by the routine. 4 ERROR INDICATORS and WARNINGS: Errors detected by the routine. IERR = 0: successful completion. 10: maximum number of QR/QL iteration steps (50) exceeded. Warnings given by the routine. IWARN = 0: no warnings. 1: the rank of the bidiagonal J, specified by the user, has been lowered because a singular value of multi- plicity > 1 was found. 5 EXTERNAL SUBROUTINES and FUNCTIONS: ESTIM, CANCEL, NSINGV, QRSTEP, QLSTEP. 6 METHOD DESCRIPTION: If the upper bound THETA is not given, THETA is computed such that precisely p - RANK singular values of J are <= THETA, using a bisec- tion method. (p = min(M,N)). QRQL then proceeds as follows. The unreduced subbidiagonals of J(j), where J(j) is the transformed bidiagonal after the j-th iteration step, are classified into the following three classes: - C1 contains the subbidiag. with all sing. values > THETA, - C2 contains the subbidiag. with all sing. values <= THETA, - C3 contains the subbidiag. with sing. values > THETA and also singular values <= THETA. If C3 is empty, then the partial diagonalization is completed, and RANK is the sum of the dimension of the elements of C1. As long as C3 is not empty, QR or QL iterations are performed on each subbidiagonal of C3, until this subbidiagonal has been splitted into two subbidiagonals. Then these two submatrices are classified and the iterations are restarted. If the upper left diagonal element of the subbidiagonal is larger than its lower right diagonal element, then QL iterations are per- formed, else QR iterations are used. The shift is equal to the mini- mal diagonal element in absolute value. However, if this element is > THETA, the shift is set to zero. 7 REFERENCES: [1] S. Van Huffel, J. Vandewalle and A. Haegemans, An efficient and reliable algorithm for computing the singular subspace of a matrix associated with its smallest singular values. J. Comput. and Applied Math., 19 (1987), 313 - 330. *********************************************************************** 1988, February 15.