C C ________________________________________________________ C | | C | QR FACTOR A GENERAL MATRIX WITH COLUMN PIVOTING | C | | C | INPUT: | C | | C | A --ARRAY CONTAINING MATRIX | C | (LENGTH AT LEAST 3+N+MN + MIN(M,N)) | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | M --NUMBER OF ROWS IN COEFFICIENT MATRIX | C | | C | N --NUMBER OF COLUMNS IN COEFF. MATRIX | C | | C | OUTPUT: | C | | C | A --FACTORED MATRIX | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C | PACKAGE SUBROUTINES: RPACK | C |________________________________________________________| C SUBROUTINE QR(A,LA,M,N) REAL A(1),R,S,T,U,V,W INTEGER B,C,D,E,F,G,H,I,J,K,L,LA,M,N,O,P,Q,Z IF ( LA .GT. M ) CALL RPACK(A,LA,M,N) U = 1. 10 U = .5*U T = 1 + U IF ( T .GT. 1. ) GOTO 10 U = SQRT(40.*U) O = M + 1 I = 2 + N*O F = I IF ( M .LT. N ) F = 2 + M*O D = I - M J = 2 + N H = I + 1 C = H L = H + N G = H V = 0. C --------------------------------------- C |*** COMPUTE 2-NORM OF EACH COLUMN ***| C --------------------------------------- 20 S = 0. R = S K = I - M 30 T = A(I-J) A(I) = T I = I - 1 IF ( T .NE. 0. ) GOTO 40 IF ( I .GT. K ) GOTO 30 GOTO 70 40 S = ABS(T) R = 1. IF ( I .EQ. K ) GOTO 70 50 T = A(I-J) A(I) = T I = I - 1 IF ( ABS(T) .GT. S ) GOTO 60 R = R + (T/S)**2 IF ( I .GT. K ) GOTO 50 GOTO 70 60 R = 1 + R*(S/T)**2 S = ABS(T) IF ( I .GT. K ) GOTO 50 70 W = S*SQRT(R) A(G) = W A(L) = U*W IF ( W .LT. V ) GOTO 80 V = W P = G 80 L = L - 1 G = I I = I - 1 J = J - 1 K = I - M IF ( K .GT. 2 ) GOTO 20 A(1) = 3230 A(2) = M A(3) = N K = 4 G = M C ----------------------------- C |*** START FACTORIZATION ***| C ----------------------------- 90 E = K + G L = P - M J = P - E H = H + 1 A(H+J/O) = A(H) A(H) = (P-3)/O C ----------------------------- C |*** INTERCHANGE COLUMNS ***| C ----------------------------- DO 100 I = L,P T = A(I) Q = I - J A(I) = A(Q) 100 A(Q) = T S = 0. T = A(E) IF ( T .EQ. 0. ) GOTO 240 IF ( K .EQ. F ) RETURN E = E - 1 DO 110 I = K,E 110 S = S + (A(I)/T)**2 IF ( S .EQ. 0. ) GOTO 240 S = T*SQRT(S) T = A(K) A(K) = S IF ( T .GE. 0. ) A(K) = -S R = 1./SQRT(S*(S+ABS(T))) I = E C ---------------------------------- C |*** STORE HOUSEHOLDER MATRIX ***| C ---------------------------------- 120 A(I+1) = R*A(I) I = I - 1 IF ( I .GT. K ) GOTO 120 IF ( T .LT. 0. ) S = -S A(K+1) = R*(T+S) IF ( K .GT. D ) RETURN J = K E = -1 Z = H G = G - 1 V = 0. 130 IF ( J .GT. D ) GOTO 230 J = J + O E = E + O L = J + G Z = Z + 1 B = L + 1 S = A(B) T = 0. IF ( S .EQ. 0. ) GOTO 160 DO 140 I = J,L 140 T = T + A(I)*A(I-E) C ------------------------------ C |*** UPDATE FACTORIZATION ***| C ------------------------------ DO 150 I = J,L 150 A(I) = A(I) - T*A(I-E) T = S*SQRT(ABS(1.-(A(J)/S)**2)) A(B) = T IF ( T .LT. A(Z) ) GOTO 170 160 IF ( T .LT. V ) GOTO 130 V = T P = B GOTO 130 170 I = J + 1 C ------------------------------- C |*** COMPUTE COLUMN 2-NORM ***| C ------------------------------- S = 0. R = S DO 180 Q = I,L 180 IF ( A(Q) .NE. 0. ) GOTO 190 GOTO 220 190 S = ABS(A(Q)) DO 210 I = Q,L T = ABS(A(I)) IF ( T .GT. S ) GOTO 200 R = R + (T/S)**2 GOTO 210 200 R = 1 + R*(S/T)**2 S = T 210 CONTINUE 220 T = S*SQRT(R) A(B) = T A(Z) = U*T IF ( T .LT. V ) GOTO 130 V = T P = B GOTO 130 230 K = K + O + 1 IF ( K .LT. F ) GOTO 90 IF ( M .GE. N ) GOTO 90 RETURN 240 A(H) = 0. A(1) = -3230 RETURN END