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 ap = qr. the householder transformation for column k, k = 1,2,...,min(m,n), is of the form
t
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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | m | a positive integer input variable set to the number of rows of a. |
||
| integer, | intent(in) | :: | n | a positive integer input variable set to the number of columns of a. |
||
| real(kind=wp), | intent(inout) | :: | a(Lda,n) | 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). |
||
| integer, | intent(in) | :: | Lda | a positive integer input variable not less than m which specifies the leading dimension of the array a. |
||
| logical, | intent(in) | :: | Pivot | 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. |
||
| integer, | intent(out) | :: | Ipvt(Lipvt) | an integer output array of length lipvt. ipvt defines the permutation matrix p such that ap = qr. column j of p is column ipvt(j) of the identity matrix. if pivot is false, ipvt is not referenced. |
||
| integer, | intent(in) | :: | Lipvt | 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. |
||
| real(kind=wp), | intent(out) | :: | Rdiag(n) | an output array of length n which contains the diagonal elements of r. |
||
| real(kind=wp), | intent(out) | :: | Acnorm(n) | 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. |
||
| real(kind=wp), | intent(inout) | :: | Wa(n) | a work array of length n. if pivot is false, then wa can coincide with rdiag. |
subroutine qrfac(m, n, a, Lda, Pivot, Ipvt, Lipvt, Rdiag, Acnorm, Wa)
implicit none
integer, intent(in) :: m !! a positive integer input variable set to the number
!! of rows of a.
integer, intent(in) :: n !! a positive integer input variable set to the number
!! of columns of a.
integer, intent(in) :: Lda !! a positive integer input variable not less than m
!! which specifies the leading dimension of the array a.
integer, intent(in) :: Lipvt !! 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.
integer, intent(out) :: Ipvt(Lipvt) !! 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.
logical, intent(in) :: Pivot !! 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.
real(wp), intent(inout) :: a(Lda, n) !! 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).
real(wp), intent(out) :: Rdiag(n) !! an output array of length n which contains the
!! diagonal elements of r.
real(wp), intent(out) :: Acnorm(n) !! 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.
real(wp), intent(inout) :: Wa(n) !! a work array of length n. if pivot is false, then wa
!! can coincide with rdiag.
integer :: i, j, jp1, k, kmax, minmn
real(wp) :: ajnorm, sum, temp
real(wp), parameter :: p05 = 5.0e-2_wp
! 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 (Pivot) then
! 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) then
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
end if
end if
! 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) then
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) then
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. (.not. Pivot .or. Rdiag(k) == zero)) then
temp = a(j, k)/Rdiag(k)
Rdiag(k) = Rdiag(k)*sqrt(max(zero, one - temp**2))
if (p05*(Rdiag(k)/Wa(k))**2 <= epsmch) then
Rdiag(k) = enorm(m - j, a(jp1, k))
Wa(k) = Rdiag(k)
end if
end if
end do
end if
end if
Rdiag(j) = -ajnorm
end do
end subroutine qrfac