SUBROUTINE Q7RFH(IERR, IPIVOT, N, NN, NOPIVK, P, Q, R, RLEN, W)
C
C *** COMPUTE QR FACTORIZATION VIA HOUSEHOLDER TRANSFORMATIONS
C *** WITH COLUMN PIVOTING ***
C
C *** PARAMETER DECLARATIONS ***
C
INTEGER IERR, N, NN, NOPIVK, P, RLEN
INTEGER IPIVOT(P)
REAL Q(NN,P), R(RLEN), W(P)
C DIMENSION R(P*(P+1)/2)
C
C---------------------------- DESCRIPTION ----------------------------
C
C THIS ROUTINE COMPUTES A QR FACTORIZATION (VIA HOUSEHOLDER TRANS-
C FORMATIONS) OF THE MATRIX A THAT ON INPUT IS STORED IN Q.
C IF NOPIVK ALLOWS IT, THIS ROUTINE DOES COLUMN PIVOTING -- IF
C K .GT. NOPIVK, THEN ORIGINAL COLUMN K IS ELIGIBLE FOR PIVOTING.
C THE Q AND R RETURNED ARE SUCH THAT COLUMN I OF Q*R EQUALS
C COLUMN IPIVOT(I) OF THE ORIGINAL MATRIX A. THE UPPER TRIANGULAR
C MATRIX R IS STORED COMPACTLY BY COLUMNS, I.E., THE OUTPUT VECTOR R
C CONTAINS R(1,1), R(1,2), R(2,2), R(1,3), R(2,3), ..., R(P,P) (IN
C THAT ORDER). IF ALL GOES WELL, THEN THIS ROUTINE SETS IERR = 0.
C BUT IF (PERMUTED) COLUMN K OF A IS LINEARLY DEPENDENT ON
C (PERMUTED) COLUMNS 1,2,...,K-1, THEN IERR IS SET TO K AND THE R
C MATRIX RETURNED HAS R(I,J) = 0 FOR I .GE. K AND J .GE. K.
C THE ORIGINAL MATRIX A IS AN N BY P MATRIX. NN IS THE LEAD
C DIMENSION OF THE ARRAY Q AND MUST SATISFY NN .GE. N. NO
C PARAMETER CHECKING IS DONE.
C PIVOTING IS DONE AS THOUGH ALL COLUMNS OF Q WERE FIRST
C SCALED TO HAVE THE SAME NORM. IF COLUMN K IS ELIGIBLE FOR
C PIVOTING AND ITS (SCALED) NORM**2 LOSS IS MORE THAN THE
C MINIMUM SUCH LOSS (OVER COLUMNS K THRU P), THEN COLUMN K IS
C SWAPPED WITH THE COLUMN OF LEAST NORM**2 LOSS.
C
C CODED BY DAVID M. GAY (FALL 1979, SPRING 1984).
C
C-------------------------- LOCAL VARIABLES --------------------------
C
INTEGER I, II, J, K, KK, KM1, KP1, NK1
REAL AK, QKK, S, SINGTL, T, T1, WK
REAL D7TPR, R7MDC, V2NRM
EXTERNAL D7TPR, R7MDC, V2AXY, V7SCL, V7SCP, V7SWP, V2NRM
C/+
REAL SQRT
C/
REAL BIG, BIGRT, MEPS10, ONE, TEN, TINY, TINYRT,
1 WTOL, ZERO
C/6
C DATA ONE/1.0E+0/, TEN/1.E+1/, WTOL/0.75E+0/, ZERO/0.0E+0/
C/7
PARAMETER (ONE=1.0E+0, TEN=1.E+1, WTOL=0.75E+0, ZERO=0.0E+0)
SAVE BIGRT, MEPS10, TINY, TINYRT
C/
DATA BIGRT/0.0E+0/, MEPS10/0.0E+0/, TINY/0.E+0/, TINYRT/0.E+0/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
IERR = 0
IF (MEPS10 .GT. ZERO) GO TO 10
BIGRT = R7MDC(5)
MEPS10 = TEN * R7MDC(3)
TINYRT = R7MDC(2)
TINY = R7MDC(1)
BIG = R7MDC(6)
IF (TINY*BIG .LT. ONE) TINY = ONE / BIG
10 SINGTL = FLOAT(MAX0(N,P)) * MEPS10
C
C *** INITIALIZE W, IPIVOT, AND DIAG(R) ***
C
J = 0
DO 40 I = 1, P
IPIVOT(I) = I
T = V2NRM(N, Q(1,I))
IF (T .GT. ZERO) GO TO 20
W(I) = ONE
GO TO 30
20 W(I) = ZERO
30 J = J + I
R(J) = T
40 CONTINUE
C
C *** MAIN LOOP ***
C
KK = 0
NK1 = N + 1
DO 130 K = 1, P
IF (NK1 .LE. 1) GO TO 999
NK1 = NK1 - 1
KK = KK + K
KP1 = K + 1
IF (K .LE. NOPIVK) GO TO 60
IF (K .GE. P) GO TO 60
C
C *** FIND COLUMN WITH MINIMUM WEIGHT LOSS ***
C
T = W(K)
IF (T .LE. ZERO) GO TO 60
J = K
DO 50 I = KP1, P
IF (W(I) .GE. T) GO TO 50
T = W(I)
J = I
50 CONTINUE
IF (J .EQ. K) GO TO 60
C
C *** INTERCHANGE COLUMNS K AND J ***
C
I = IPIVOT(K)
IPIVOT(K) = IPIVOT(J)
IPIVOT(J) = I
W(J) = W(K)
W(K) = T
I = J*(J+1)/2
T1 = R(I)
R(I) = R(KK)
R(KK) = T1
CALL V7SWP(N, Q(1,K), Q(1,J))
IF (K .LE. 1) GO TO 60
I = I - J + 1
J = KK - K + 1
CALL V7SWP(K-1, R(I), R(J))
C
C *** COLUMN K OF Q SHOULD BE NEARLY ORTHOGONAL TO THE PREVIOUS
C *** COLUMNS. NORMALIZE IT, TEST FOR SINGULARITY, AND DECIDE
C *** WHETHER TO REORTHOGONALIZE IT.
C
60 AK = R(KK)
IF (AK .LE. ZERO) GO TO 140
WK = W(K)
C
C *** SET T TO THE NORM OF (Q(K,K),...,Q(N,K))
C *** AND CHECK FOR SINGULARITY.
C
IF (WK .LT. WTOL) GO TO 70
T = V2NRM(NK1, Q(K,K))
IF (T / AK .LE. SINGTL) GO TO 140
GO TO 80
70 T = SQRT(ONE - WK)
IF (T .LE. SINGTL) GO TO 140
T = T * AK
C
C *** DETERMINE HOUSEHOLDER TRANSFORMATION ***
C
80 QKK = Q(K,K)
IF (T .LE. TINYRT) GO TO 90
IF (T .GE. BIGRT) GO TO 90
IF (QKK .LT. ZERO) T = -T
QKK = QKK + T
S = SQRT(T * QKK)
GO TO 110
90 S = SQRT(T)
IF (QKK .LT. ZERO) GO TO 100
QKK = QKK + T
S = S * SQRT(QKK)
GO TO 110
100 T = -T
QKK = QKK + T
S = S * SQRT(-QKK)
110 Q(K,K) = QKK
C
C *** SCALE (Q(K,K),...,Q(N,K)) TO HAVE NORM SQRT(2) ***
C
IF (S .LE. TINY) GO TO 140
CALL V7SCL(NK1, Q(K,K), ONE/S, Q(K,K))
C
R(KK) = -T
C
C *** COMPUTE R(K,I) FOR I = K+1,...,P AND UPDATE Q ***
C
IF (K .GE. P) GO TO 999
J = KK + K
II = KK
DO 120 I = KP1, P
II = II + I
CALL V2AXY(NK1, Q(K,I), - D7TPR(NK1,Q(K,K),Q(K,I)),
1 Q(K,K), Q(K,I))
T = Q(K,I)
R(J) = T
J = J + I
T1 = R(II)
IF (T1 .GT. ZERO) W(I) = W(I) + (T/T1)**2
120 CONTINUE
130 CONTINUE
C
C *** SINGULAR Q ***
C
140 IERR = K
KM1 = K - 1
J = KK
DO 150 I = K, P
CALL V7SCP(I-KM1, R(J), ZERO)
J = J + I
150 CONTINUE
C
999 RETURN
C *** LAST CARD OF Q7RFH FOLLOWS ***
END
.