*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