SUBROUTINE QRFAC(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
! This subroutine uses householder transformations with column
! pivoting (optional) to compute a QR factorization of the
! M by N matrix A. That is, QRFAC determines an orthogonal
! matrix Q, a permutation matrix P, and an upper trapezoidal
! matrix R with diagonal elements of nonincreasing magnitude,
! such that A*P = Q*R. The householder transformation for
! column K, K = 1,2,...,min(M,N), is of the form
!
! i - (1/u(k))*u*u
!
! where U has zeros in the first K-1 positions. The form of
! this transformation and the method of pivoting first
! appeared in the corresponding Linpack subroutine.
!
! The subroutine statement is
!
! subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
!
! where
!
! M is a positive integer input variable set to the number
! of rows of A.
!
! N is a positive integer input variable set to the number
! of columns of A.
!
! A is an M by N array. On input A contains the matrix for
! which the QR factorization is to be computed. On output
! the strict upper trapezoidal part of A contains the strict
! upper trapezoidal part of R, and the lower trapezoidal
! part of A contains a factored form of Q (the non-trivial
! elements of the U vectors described above).
!
! LDA is a positive integer input variable not less than M
! which specifies the leading dimension of the array A.
!
! PIVOT is a logical input variable. If PIVOT is set true,
! then column pivoting is enforced. If PIVOT is set false,
! then no column pivoting is done.
!
! IPVT is an integer output array of length LIPVT. IPVT
! defines the permutation matrix P such that A*P = Q*R.
! Column J of P is column IPVT(J) of the identity matrix.
! If PIVOT is false, IPVT is not referenced.
!
! LIPVT is a positive integer input variable. If PIVOT is false,
! then LIPVT may be as small as 1. If PIVOT is true, then
! LIPVT must be at least N.
!
! RDIAG is an output array of length N which contains the
! diagonal elements of R.
!
! ACNORM is an output array of length N which contains the
! norms of the corresponding columns of the input matrix A.
! If this information is not needed, then ACNORM can coincide
! with RDIAG.
!
! WA is a work array of length N. If PIVOT is false, then WA
! can coincide with RDIAG.
!
! Argonne National Laboratory. Minpack project. March 1980.
! Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More
IMPLICIT NONE
INTEGER m,n,lda,lipvt,i,j,jp1,k,kmax,minmn
INTEGER ipvt(lipvt)
LOGICAL pivot
real(dp) a(lda,n),rdiag(n),acnorm(n),wa(n)
real(dp) ajnorm,epsmch,one,p05,sum,temp,zero
one = 1.0D0
p05 = 0.05D0
zero = 0.0D0
! Machine precision
epsmch = spmpar(1)
! Compute the initial column norms and initialize several arrays.
do j = 1, n
acnorm(j) = enorm(m,a(1,j))
rdiag(j) = acnorm(j)
wa(j) = rdiag(j)
if (pivot) ipvt(j) = j
end do
! Reduce a to r with householder transformations.
minmn = min(m,n)
do j = 1, minmn
if (.not.pivot) goto 40
! Bring the column of largest norm into the pivot position.
kmax = j
do k = j, n
if (rdiag(k) > rdiag(kmax)) kmax = k
end do
if (kmax == j) goto 40
do i = 1, m
temp = a(i,j)
a(i,j) = a(i,kmax)
a(i,kmax) = temp
end do
rdiag(kmax) = rdiag(j)
wa(kmax) = wa(j)
k = ipvt(j)
ipvt(j) = ipvt(kmax)
ipvt(kmax) = k
40 continue
! Compute the householder transformation to reduce the
! j-th column of a to a multiple of the j-th unit vector.
ajnorm = enorm(m-j+1,a(j,j))
if (ajnorm == zero) goto 100
if (a(j,j) < zero) ajnorm = -ajnorm
do i = j, m
a(i,j) = a(i,j)/ajnorm
end do
a(j,j) = a(j,j) + one
! Apply the transformation to the remaining columns
! and update the norms.
jp1 = j + 1
if (n < jp1) goto 100
do k = jp1, n
sum = zero
do i = j, m
sum = sum + a(i,j)*a(i,k)
end do
temp = sum/a(j,j)
do i = j, m
a(i,k) = a(i,k) - temp*a(i,j)
end do
if ((.not.pivot).or.(rdiag(k) == zero)) goto 80
temp = a(j,k)/rdiag(k)
rdiag(k) = rdiag(k)*sqrt(max(zero,one-temp**2))
if ((p05*(rdiag(k)/wa(k))**2) > epsmch) goto 80
rdiag(k) = enorm(m-j,a(jp1,k))
wa(k) = rdiag(k)
80 continue
end do
100 continue
rdiag(j) = -ajnorm
end do
return
end SUBROUTINE QRFAC