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