*DECK DQRFAC
SUBROUTINE DQRFAC (M, N, A, LDA, PIVOT, IPVT, LIPVT, SIGMA,
+ ACNORM, WA)
C***BEGIN PROLOGUE DQRFAC
C***SUBSIDIARY
C***PURPOSE Subsidiary to DNLS1, DNLS1E, DNSQ and DNSQE
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (QRFAC-S, DQRFAC-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
C
C **** Double Precision version of QRFAC ****
C
C This subroutine uses Householder transformations with column
C pivoting (optional) to compute a QR factorization of the
C M by N matrix A. That is, DQRFAC determines an orthogonal
C matrix Q, a permutation matrix P, and an upper trapezoidal
C matrix R with diagonal elements of nonincreasing magnitude,
C such that A*P = Q*R. The Householder transformation for
C column K, K = 1,2,...,MIN(M,N), is of the form
C
C T
C I - (1/U(K))*U*U
C
C where U has zeros in the first K-1 positions. The form of
C this transformation and the method of pivoting first
C appeared in the corresponding LINPACK subroutine.
C
C The subroutine statement is
C
C SUBROUTINE DQRFAC(M,N,A,LDA,PIVOT,IPVT,LIPVT,SIGMA,ACNORM,WA)
C
C where
C
C M is a positive integer input variable set to the number
C of rows of A.
C
C N is a positive integer input variable set to the number
C of columns of A.
C
C A is an M by N array. On input A contains the matrix for
C which the QR factorization is to be computed. On output
C the strict upper trapezoidal part of A contains the strict
C upper trapezoidal part of R, and the lower trapezoidal
C part of A contains a factored form of Q (the non-trivial
C elements of the U vectors described above).
C
C LDA is a positive integer input variable not less than M
C which specifies the leading dimension of the array A.
C
C PIVOT is a logical input variable. If pivot is set .TRUE.,
C then column pivoting is enforced. If pivot is set .FALSE.,
C then no column pivoting is done.
C
C IPVT is an integer output array of length LIPVT. IPVT
C defines the permutation matrix P such that A*P = Q*R.
C Column J of P is column IPVT(J) of the identity matrix.
C If pivot is .FALSE., IPVT is not referenced.
C
C LIPVT is a positive integer input variable. If PIVOT is
C .FALSE., then LIPVT may be as small as 1. If PIVOT is
C .TRUE., then LIPVT must be at least N.
C
C SIGMA is an output array of length N which contains the
C diagonal elements of R.
C
C ACNORM is an output array of length N which contains the
C norms of the corresponding columns of the input matrix A.
C If this information is not needed, then ACNORM can coincide
C with SIGMA.
C
C WA is a work array of length N. If pivot is .FALSE., then WA
C can coincide with SIGMA.
C
C***SEE ALSO DNLS1, DNLS1E, DNSQ, DNSQE
C***ROUTINES CALLED D1MACH, DENORM
C***REVISION HISTORY (YYMMDD)
C 800301 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE DQRFAC
INTEGER M,N,LDA,LIPVT
INTEGER IPVT(*)
LOGICAL PIVOT
SAVE ONE, P05, ZERO
DOUBLE PRECISION A(LDA,*),SIGMA(*),ACNORM(*),WA(*)
INTEGER I,J,JP1,K,KMAX,MINMN
DOUBLE PRECISION AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
DOUBLE PRECISION D1MACH,DENORM
DATA ONE,P05,ZERO /1.0D0,5.0D-2,0.0D0/
C***FIRST EXECUTABLE STATEMENT DQRFAC
EPSMCH = D1MACH(4)
C
C COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
C
DO 10 J = 1, N
ACNORM(J) = DENORM(M,A(1,J))
SIGMA(J) = ACNORM(J)
WA(J) = SIGMA(J)
IF (PIVOT) IPVT(J) = J
10 CONTINUE
C
C REDUCE A TO R WITH HOUSEHOLDER TRANSFORMATIONS.
C
MINMN = MIN(M,N)
DO 110 J = 1, MINMN
IF (.NOT.PIVOT) GO TO 40
C
C BRING THE COLUMN OF LARGEST NORM INTO THE PIVOT POSITION.
C
KMAX = J
DO 20 K = J, N
IF (SIGMA(K) .GT. SIGMA(KMAX)) KMAX = K
20 CONTINUE
IF (KMAX .EQ. J) GO TO 40
DO 30 I = 1, M
TEMP = A(I,J)
A(I,J) = A(I,KMAX)
A(I,KMAX) = TEMP
30 CONTINUE
SIGMA(KMAX) = SIGMA(J)
WA(KMAX) = WA(J)
K = IPVT(J)
IPVT(J) = IPVT(KMAX)
IPVT(KMAX) = K
40 CONTINUE
C
C COMPUTE THE HOUSEHOLDER TRANSFORMATION TO REDUCE THE
C J-TH COLUMN OF A TO A MULTIPLE OF THE J-TH UNIT VECTOR.
C
AJNORM = DENORM(M-J+1,A(J,J))
IF (AJNORM .EQ. ZERO) GO TO 100
IF (A(J,J) .LT. ZERO) AJNORM = -AJNORM
DO 50 I = J, M
A(I,J) = A(I,J)/AJNORM
50 CONTINUE
A(J,J) = A(J,J) + ONE
C
C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS
C AND UPDATE THE NORMS.
C
JP1 = J + 1
IF (N .LT. JP1) GO TO 100
DO 90 K = JP1, N
SUM = ZERO
DO 60 I = J, M
SUM = SUM + A(I,J)*A(I,K)
60 CONTINUE
TEMP = SUM/A(J,J)
DO 70 I = J, M
A(I,K) = A(I,K) - TEMP*A(I,J)
70 CONTINUE
IF (.NOT.PIVOT .OR. SIGMA(K) .EQ. ZERO) GO TO 80
TEMP = A(J,K)/SIGMA(K)
SIGMA(K) = SIGMA(K)*SQRT(MAX(ZERO,ONE-TEMP**2))
IF (P05*(SIGMA(K)/WA(K))**2 .GT. EPSMCH) GO TO 80
SIGMA(K) = DENORM(M-J,A(JP1,K))
WA(K) = SIGMA(K)
80 CONTINUE
90 CONTINUE
100 CONTINUE
SIGMA(J) = -AJNORM
110 CONTINUE
RETURN
C
C LAST CARD OF SUBROUTINE DQRFAC.
C
END