2 SUBROUTINE QRFAC (M, N, A, LDA, PIVOT, IPVT, LIPVT, SIGMA, ACNORM,
4 C***BEGIN PROLOGUE QRFAC
6 C***PURPOSE Subsidiary to SNLS1, SNLS1E, SNSQ and SNSQE
8 C***TYPE SINGLE PRECISION (QRFAC-S, DQRFAC-D)
12 C This subroutine uses Householder transformations with column
13 C pivoting (optional) to compute a QR factorization of the
14 C M by N matrix A. That is, QRFAC determines an orthogonal
15 C matrix Q, a permutation matrix P, and an upper trapezoidal
16 C matrix R with diagonal elements of nonincreasing magnitude,
17 C such that A*P = Q*R. The Householder transformation for
18 C column K, K = 1,2,...,MIN(M,N), is of the form
23 C where U has zeros in the first K-1 positions. The form of
24 C this transformation and the method of pivoting first
25 C appeared in the corresponding LINPACK subroutine.
27 C The subroutine statement is
29 C SUBROUTINE QRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
33 C M is a positive integer input variable set to the number
36 C N is a positive integer input variable set to the number
39 C A is an M by N array. On input A contains the matrix for
40 C which the QR factorization is to be computed. On output
41 C the strict upper trapezoidal part of A contains the strict
42 C upper trapezoidal part of R, and the lower trapezoidal
43 C part of A contains a factored form of Q (the non-trivial
44 C elements of the U vectors described above).
46 C LDA is a positive integer input variable not less than M
47 C which specifies the leading dimension of the array A.
49 C PIVOT is a logical input variable. If pivot is set .TRUE.,
50 C then column pivoting is enforced. If pivot is set .FALSE.,
51 C then no column pivoting is done.
53 C IPVT is an integer output array of length LIPVT. IPVT
54 C defines the permutation matrix P such that A*P = Q*R.
55 C Column J of P is column IPVT(J) of the identity matrix.
56 C If pivot is .FALSE., IPVT is not referenced.
58 C LIPVT is a positive integer input variable. If PIVOT is
59 C .FALSE., then LIPVT may be as small as 1. If PIVOT is
60 C .TRUE., then LIPVT must be at least N.
62 C SIGMA is an output array of length N which contains the
63 C diagonal elements of R.
65 C ACNORM is an output array of length N which contains the
66 C norms of the corresponding columns of the input matrix A.
67 C If this information is not needed, then ACNORM can coincide
70 C WA is a work array of length N. If pivot is .FALSE., then WA
71 C can coincide with SIGMA.
73 C***SEE ALSO SNLS1, SNLS1E, SNSQ, SNSQE
74 C***ROUTINES CALLED ENORM, R1MACH
75 C***REVISION HISTORY (YYMMDD)
77 C 890531 Changed all specific intrinsics to generic. (WRB)
78 C 890831 Modified array declarations. (WRB)
79 C 891214 Prologue converted to Version 4.0 format. (BAB)
80 C 900326 Removed duplicate information from DESCRIPTION section.
82 C 900328 Added TYPE section. (WRB)
83 C***END PROLOGUE QRFAC
87 REAL A(LDA,*),SIGMA(*),ACNORM(*),WA(*)
88 INTEGER I,J,JP1,K,KMAX,MINMN
89 REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
92 DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/
93 C***FIRST EXECUTABLE STATEMENT QRFAC
96 C COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
99 ACNORM(J) = ENORM(M,A(1,J))
102 IF (PIVOT) IPVT(J) = J
105 C REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
109 IF (.NOT.PIVOT) GO TO 40
111 C BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
115 IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
117 IF (KMAX .EQ. J) GO TO 40
123 SIGMA(KMAX) = SIGMA(J)
130 C COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
131 C J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
133 AJNORM = ENORM(M-J+1,A(J,J))
134 IF (AJNORM .EQ. ZERO) GO TO 100
135 IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
137 A(I,J) = A(I,J)/AJNORM
139 A(J,J) = A(J,J) + ONE
141 C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
142 C AND UPDATE THE NORMS.
145 IF (N .LT. JP1) GO TO 100
149 SUM = SUM + A(I,J)*A(I,K)
153 A(I,K) = A(I,K) - TEMP*A(I,J)
155 IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
156 TEMP = A(J,K)/SIGMA(K)
157 SIGMA(K) = SIGMA(K)*SQRT(MAX(ZERO,ONE-TEMP**2))
158 IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
159 SIGMA(K) = ENORM(M-J,A(JP1,K))
168 C LAST CARD OF SUBROUTINE QRFAC.