SUBROUTINE PTLS(C, LDC, M, N, L, RANK, THETA, X, LDX, Q, INUL, * WRK, IWRK, LWRK, TOL1, TOL2, IERR, IWARN) C C CONTRIBUTOR: S. Van Huffel (ESAT Laboratory, KU Leuven). C C REVISIONS: 1988, February 15. C C .. Scalar Arguments .. INTEGER LDC, M, N, L, RANK, LDX, IERR, IWARN DOUBLE PRECISION THETA, TOL1, TOL2 C .. Array Arguments .. INTEGER IWRK(N+L) DOUBLE PRECISION C(LDC,*), X(LDX,*), Q(*), WRK(*) LOGICAL INUL(N+L), LWRK(N+L) C .. External Subroutines/Functions .. DOUBLE PRECISION DDOT EXTERNAL DAXPY,DDOT,DCOPY,DQRDC, BIDIAG, CANCEL, HOUSH, INIT, QRQL C .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN C .. Local Scalars .. INTEGER P, NL, PP1, J1, MC, NJ, K, II, I1, N1, MNL1 DOUBLE PRECISION TEMP, HH, INPROD LOGICAL ZERO C .. Local Arrays .. INTEGER DUMAR1(1) DOUBLE PRECISION DUMAR2(1,1) C .. Executable Statements .. C IERR = 0 IWARN = 0 IF (M .LT. 1) IERR = 1 IF (N .LT. 1) IERR = 2 IF (L .LT. 1) IERR = 3 IF (LDC .LT. MAX(M,N+L)) IERR = 4 IF (LDX .LT. N) IERR = 5 IF (RANK .GT. MIN(M,N)) IERR = 6 IF ((RANK .LT. 0) .AND. (THETA .LT. 0.0D0)) IERR = 7 IF (TOL1 .LT. 0.0D0) IERR = 8 IF (TOL2 .LT. 0.0D0) IERR = 9 IF (IERR .NE. 0) RETURN C C Initializations. C NL = N + L MNL1 = M + NL + 1 P = MIN(M,NL) DO 1 I = 1, P INUL(I) = .FALSE. LWRK(I) = .FALSE. 1 CONTINUE PP1 = P + 1 DO 2 I = PP1, NL INUL(I) = .TRUE. LWRK(I) = .FALSE. 2 CONTINUE C DO 3 J = 1, L DO 3 I = 1, N X(I,J) = 0.0D0 3 CONTINUE C C Subroutine PTLS solves a set of linear equations by a Total Least C Squares Approximation, based on the Partial SVD. C C Step 1: Bidiagonalization phase C ----------------------- C 1.a): If M >= 5*(N+L)/3, transform C into upper triangular form R C by Householder transformations. C MC = M IF (3*M .GE. 5*NL) THEN CALL DQRDC(C, LDC, M, NL, Q, DUMAR1, WRK, 0) DO 4 J = 1, NL J1 = J + 1 DO 4 I = J1, NL C(I,J) = 0.0D0 4 CONTINUE MC = NL END IF C C 1.b): Transform C (or R) into bidiagonal form Q using Householder C transformations. C CALL BIDIAG(C, LDC, MC, NL, Q, Q(P+1), WRK) C C Store the Householder transformations performed onto the rows of C C in the last storage locations of the work array WRK. C MC = MIN(NL-2,M) IF (MC .GT. 0) THEN K = MNL1 DO 5 II = 1, MC J = MC - II + 1 NJ = NL - J CALL DCOPY(NJ, C(J,J+1), LDC, WRK(K), 1) K = K + NJ 5 CONTINUE END IF C C 1.c): Initialize the right singular base matrix V with the identi- C ty matrix (V overwrites C). C CALL INIT(C, LDC, NL, NL) C C 1.d): If M < N+L, bring the bidiagonal Q to M by M by cancelling C its last superdiagonal element using Givens rotations. C PP1 = P IF (M .LT. NL) THEN PP1 = P + 1 CALL CANCEL(DUMAR2, 1, C, LDC, Q, Q(PP1), P, PP1, PP1, PP1, * TOL2, .FALSE., .TRUE.) END IF C C REPEAT C C Compute the Householder matrix Q and matrices F and Y such that C F is nonsingular. C Step 2: Partial diagonalization phase. C ----------------------------- C Diagonalize the bidiagonal Q partially until convergence to C the desired right singular subspace. C 6 CALL QRQL(DUMAR2, 1, C, LDC, P, PP1, RANK, THETA, Q, Q(P+1), * INUL, TOL1, TOL2, .FALSE., .TRUE., IERR, IWARN) C IF (IERR.NE.0) RETURN C C Step 3: Back transformation phase. C ------------------------- C Apply the Householder transformations (stored in WRK) perfor- C med onto the rows of C during the bidiagonalization phase, to C the selected base vectors in the right singular base matrix C of C. C DO 8 I = 1, NL IF (INUL(I) .AND. (.NOT. LWRK(I))) THEN K = MNL1 DO 7 II = 1, MC J = MC - II + 1 NJ = NL - J J1 = J + 1 IF (ABS(WRK(K)) .GT. TOL2) THEN TEMP = -DDOT(NJ, WRK(K), 1, C(J1,I), 1)/WRK(K) CALL DAXPY(NJ, TEMP, WRK(K), 1, C(J1,I), 1) K = K + NJ END IF 7 CONTINUE LWRK(I) = .TRUE. END IF 8 CONTINUE IF (RANK.LE.0) RETURN C C Step 4: Compute matrices F and Y using Householder transf. Q. C ------------------------ K = 0 DO 9 I = 1, NL IF (INUL(I)) THEN K = K + 1 IWRK(K) = I END IF 9 CONTINUE C IF (K .LT. L) THEN C C Rank TLS approximation is larger than min(M,N). C IERR = 3 RETURN END IF N1 = N + 1 ZERO = .FALSE. I = NL C C WHILE ((K > 1) .AND. (I > N) .AND. (.NOT.ZERO)) DO 10 IF ((K.GT.1) .AND. (I.GT.N) .AND. (.NOT.ZERO)) THEN DO 11 J = 1, K WRK(J) = C(I,IWRK(J)) 11 CONTINUE C C Compute Householder transformation. C CALL HOUSH(WRK, K, K, TOL2, ZERO, TEMP) IF (.NOT.ZERO) THEN C C Apply Householder transformation onto the selected base C vectors. C DO 13 I1 = 1, I INPROD = 0.0D0 DO 12 J = 1, K INPROD = INPROD + WRK(J) * C(I1,IWRK(J)) 12 CONTINUE HH = INPROD * TEMP DO 13 J = 1, K J1 = IWRK(J) C(I1,J1) = C(I1,J1) - WRK(J) * HH 13 CONTINUE C K = K - 1 END IF I = I - 1 GOTO 10 END IF C END WHILE 10 C IF (.NOT.ZERO) K = N1 - RANK C C If F singular, lower the rank of the TLS approximation . C IF (ABS(C(N1,IWRK(K))) .LE. TOL2) THEN RANK = RANK - 1 IWARN = 2 THETA = -1.0D0 GO TO 6 END IF C UNTIL ((.NOT.ZERO) .AND. (F nonsingular)) C C Step 5: Compute TLS solution C -------------------- C Solve X F = -Y by forward elimination (F is upper triangular). NJ = IWRK(K) CALL DAXPY(N, -1.0D0/C(N1,NJ), C(1,NJ), 1, X, 1) DO 15 J = 2, L J1 = J - 1 NJ = IWRK(K+J1) TEMP = C(N+J,NJ) DO 14 I = 1, N X(I,J) = -(C(I,NJ) + DDOT(J1,C(N1,NJ),1,X(I,1),LDX))/TEMP 14 CONTINUE 15 CONTINUE RETURN C *** Last line of PTLS *********************************************** END