+++ /dev/null
-*DECK QRFAC
- SUBROUTINE QRFAC (M, N, A, LDA, PIVOT, IPVT, LIPVT, SIGMA, ACNORM,
- + WA)
-C***BEGIN PROLOGUE QRFAC
-C***SUBSIDIARY
-C***PURPOSE Subsidiary to SNLS1, SNLS1E, SNSQ and SNSQE
-C***LIBRARY SLATEC
-C***TYPE SINGLE PRECISION (QRFAC-S, DQRFAC-D)
-C***AUTHOR (UNKNOWN)
-C***DESCRIPTION
-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, QRFAC 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 QRFAC(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 SNLS1, SNLS1E, SNSQ, SNSQE
-C***ROUTINES CALLED ENORM, R1MACH
-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 QRFAC
- INTEGER M,N,LDA,LIPVT
- INTEGER IPVT(*)
- LOGICAL PIVOT
- REAL A(LDA,*),SIGMA(*),ACNORM(*),WA(*)
- INTEGER I,J,JP1,K,KMAX,MINMN
- REAL AJNORM,EPSMCH,ONE,P05,SUM,TEMP,ZERO
- REAL R1MACH,ENORM
- SAVE ONE, P05, ZERO
- DATA ONE,P05,ZERO /1.0E0,5.0E-2,0.0E0/
-C***FIRST EXECUTABLE STATEMENT QRFAC
- EPSMCH = R1MACH(4)
-C
-C COMPUTE THE INITIAL COLUMN NORMS AND INITIALIZE SEVERAL ARRAYS.
-C
- DO 10 J = 1, N
- ACNORM(J) = ENORM(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 = ENORM(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) = ENORM(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 QRFAC.
-C
- END
git://source.jalview.org / jabaws.git / blobdiff