SUBROUTINE CANCEL(U, LDU, V, LDV, Q, E, M, N, I, K, TOL, * WANTU, WANTV) C C PURPOSE: C C Either, subroutine CANCEL separates a zero singular value of a C subbidiagonal matrix of order k, k <= p, of the bidiagonal 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 with p = min(M,N), by annihilating one or two superdiagonal C elements E(i) and/or E(i+1). C Or, CANCEL annihilates the element E(M+1) of the bidiagonal matrix C C !Q(1) E(2) 0 ... 0 0 ! C ! 0 Q(2) E(3) . . ! C J = ! . . . ! C ! . E(M) . ! C ! 0 ... Q(M) E(M+1)! C C ARGUMENT LIST: C C U - DOUBLE PRECISION array of DIMENSION (LDU,p). C On entry, U contains the M by p (p = min(M,N)) left trans- C formation matrix. C On return, the Givens rotations S on the left, annihilating C E(i+1), 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,s). C On entry, V contains the N by s (s = min(M+1,N)) right trans- C formation matrix. C On return, the Givens rotations T on the right, annihilating C E(i), have been postmultiplied into V. C NOTE: V is not referenced if WANTV = .FALSE. . C LDV - INTEGER. C LDV is the leading dimension of the array V (LDV >= N). C Q - DOUBLE PRECISION array of DIMENSION (p). C On entry, Q contains the diagonal entries of the bidiagonal J. C p = min(M,N). C On return, Q contains the diagonal elements of the transformed C bidiagonal S' J T. C E - DOUBLE PRECISION array of DIMENSION (s). C On entry, E(i), i=2,...,s, contain the superdiagonal entries C of the bidiagonal J. s = min(M+1,N), E(1) = 0.0D0. C On return, E contains the superdiagonal elements of the trans- C formed bidiagonal S' J T. 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 Either, I is the index of the negligible diagonal entry Q(I) C of the bidiagonal J, i,e. abs(Q(I)) <= TOL, I <= p. C Or, I = M + 1 if E(M+1) is to be annihilated. C K - INTEGER. C Either, K is the index of the last diagonal entry of the con- C sidered subbidiagonal of J, i.e. abs(E(K+1)) <= TOL, K <= p. C Or, K = M + 1 if E(M+1) is to be annihilated. C TOL - DOUBLE PRECISION. C Specifies that matrix elements Q(i), which are <= TOL in C absolute value, are considered to be zero. C WANTU - LOGICAL. C Logical indicating the need for postmultiplying the Givens C rotations S on the left into U. C WANTV - LOGICAL. C Logical indicating the need for postmultiplying the Givens C rotations T on the right into V. C C EXTERNAL SUBROUTINES and FUNCTIONS: C C DROT from BLAS. C C METHOD DESCRIPTION: C C Let the considered subbidiagonal be C C !Q(1) E(2) 0 ... 0 ! C ! 0 Q(2) E(3) ... ! C ! . ... ! C ! Q(i-1) E(i) . ! C Jk = ! Q(i) E(i+1) . ! C ! Q(i+1) . ! C ! . .. ! C ! . E(k)! C ! 0 ... ... Q(k)! C C A zero singular value of Jk manifests itself by a zero diagonal C entry Q(i) or in practice, a negligible value of Q(i). C We call Q(i) negligible if abs(Q(i)) <= TOL. C When such a negligible diagonal element Q(i) in Jk is present, C the subbidiagonal Jk is splitted by the routine CANCEL into 2 or C 3 unreduced subbidiagonals by annihilating E(i+1) (if i1) C using Givens rotations T on the right until Jk is reduced to the C form : C C !Q(1) E(2) 0 ... 0 ! C ! 0 . ... ! C ! . ... ! C ! Q(i-1) 0 . ! C S' Jk T = ! 0 0 . ! C ! Q(i+1) . ! C ! . .. ! C ! . E(k)! C ! 0 ... ... Q(k)! C C For more details, see [1, pp.11.12-11.14]. C The case of the annihilation of E(M+1) can be treated by the same C process. This may be seen by augmenting the matrix J with an extra C row of zeros, i.e. by introducing Q(M+1) = 0. C C REFERENCES: C C [1] J.J. Dongarra, J.R. Bunch, C.B. Moler and G.W. Stewart, C LINPACK User's Guide. SIAM, Philadelphia (1979). 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 TOL LOGICAL WANTU, WANTV C .. Array Arguments .. DOUBLE PRECISION U(LDU,*), V(LDV,*), Q(*), E(*) C .. External Subroutines/Functions .. EXTERNAL DROT C .. Intrinsic Functions .. INTRINSIC ABS, SQRT C .. Local Scalars .. INTEGER I1, L, L1 DOUBLE PRECISION C, S, F, G, H C .. Executable Statements .. C IF (I .LE. M) Q(I) = 0.0D0 C C Annihilate E(I+1) (if I < K). C IF (I .LT. K) THEN C = 0.0D0 S = 1.0D0 I1 = I + 1 DO 1 L = I1, K G = E(L) F = S * G E(L) = C * G IF (ABS(F) .LE. TOL) GO TO 2 G = Q(L) H = SQRT(F**2 + G**2) Q(L) = H C = G/H S = -F/H IF (WANTU) CALL DROT(M, U(1,I), 1, U(1,L), 1, C, S) 1 CONTINUE END IF C C Annihilate E(I) (if I > 1). C 2 IF (I .GT. 1) THEN I1 = I - 1 F = E(I) E(I) = 0.0D0 DO 3 L1 = 1, I1 IF (ABS(F) .LE. TOL) RETURN L = I - L1 G = Q(L) IF (ABS(G) .LE. TOL) THEN G = 0.0D0 H = ABS(F) ELSE H = SQRT(F**2 + G**2) END IF Q(L) = H C = G/H S = -F/H G = E(L) F = S * G E(L) = C * G IF (WANTV) CALL DROT(N, V(1,I), 1, V(1,L), 1, C, S) 3 CONTINUE E(1) = 0.0D0 END IF RETURN C *** Last line of CANCEL ********************************************* END