given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b, the problem is to determine an x which solves the system
a*x = b , d*x = 0 ,
in the least squares sense.
this subroutine completes the solution of the problem if it is provided with the necessary information from the qr factorization, with column pivoting, of a. that is, if ap = qr, where p is a permutation matrix, q has orthogonal columns, and r is an upper triangular matrix with diagonal elements of nonincreasing magnitude, then qrsolv expects the full upper triangle of r, the permutation matrix p, and the first n components of (q transpose)b. the system ax = b, d*x = 0, is then equivalent to
t t
r*z = q *b , p *d*p*z = 0 ,
where x = p*z. if this system does not have full rank, then a least squares solution is obtained. on output qrsolv also provides an upper triangular matrix s such that
t t t
p *(a *a + d*d)*p = s *s .
s is computed within qrsolv and may be of separate interest.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n | a positive integer input variable set to the order of r. |
||
| real(kind=wp), | intent(inout) | :: | r(Ldr,n) | an n by n array. on input the full upper triangle must contain the full upper triangle of the matrix r. on output the full upper triangle is unaltered, and the strict lower triangle contains the strict upper triangle (transposed) of the upper triangular matrix s. |
||
| integer, | intent(in) | :: | Ldr | a positive integer input variable not less than n which specifies the leading dimension of the array r. |
||
| integer, | intent(in) | :: | Ipvt(n) | an integer input array of length n which defines the permutation matrix p such that ap = qr. column j of p is column ipvt(j) of the identity matrix. |
||
| real(kind=wp), | intent(in) | :: | Diag(n) | an input array of length n which must contain the diagonal elements of the matrix d. |
||
| real(kind=wp), | intent(in) | :: | Qtb(n) | an input array of length n which must contain the first n elements of the vector (q transpose)*b. |
||
| real(kind=wp), | intent(out) | :: | x(n) | an output array of length n which contains the least squares solution of the system ax = b, dx = 0. |
||
| real(kind=wp), | intent(out) | :: | Sdiag(n) | an output array of length n which contains the diagonal elements of the upper triangular matrix s. |
||
| real(kind=wp), | intent(inout) | :: | Wa(n) | a work array of length n. |
subroutine qrsolv(n, r, Ldr, Ipvt, Diag, Qtb, x, Sdiag, Wa)
implicit none
integer, intent(in) :: n !! a positive integer input variable set to the order of r.
integer, intent(in) :: Ldr !! a positive integer input variable not less than n
!! which specifies the leading dimension of the array r.
integer, intent(in) :: Ipvt(n) !! an integer input array of length n which defines the
!! permutation matrix p such that a*p = q*r. column j of p
!! is column ipvt(j) of the identity matrix.
real(wp), intent(inout) :: r(Ldr, n) !! an n by n array. on input the full upper triangle
!! must contain the full upper triangle of the matrix r.
!! on output the full upper triangle is unaltered, and the
!! strict lower triangle contains the strict upper triangle
!! (transposed) of the upper triangular matrix s.
real(wp), intent(in) :: Diag(n) !! an input array of length n which must contain the
!! diagonal elements of the matrix d.
real(wp), intent(in) :: Qtb(n) !! an input array of length n which must contain the first
!! n elements of the vector (q transpose)*b.
real(wp), intent(out) :: x(n) !! an output array of length n which contains the least
!! squares solution of the system a*x = b, d*x = 0.
real(wp), intent(out) :: Sdiag(n) !! an output array of length n which contains the
!! diagonal elements of the upper triangular matrix s.
real(wp), intent(inout) :: Wa(n) !! a work array of length n.
integer :: i, j, jp1, k, kp1, l, nsing
real(wp) :: cos, cotan, qtbpj, sin, sum, tan, temp
real(wp), parameter :: p5 = 5.0e-1_wp
real(wp), parameter :: p25 = 2.5e-1_wp
! copy r and (q transpose)*b to preserve input and initialize s.
! in particular, save the diagonal elements of r in x.
do j = 1, n
do i = j, n
r(i, j) = r(j, i)
end do
x(j) = r(j, j)
Wa(j) = Qtb(j)
end do
! eliminate the diagonal matrix d using a givens rotation.
do j = 1, n
! prepare the row of d to be eliminated, locating the
! diagonal element using p from the qr factorization.
l = Ipvt(j)
if (Diag(l) /= zero) then
do k = j, n
Sdiag(k) = zero
end do
Sdiag(j) = Diag(l)
! the transformations to eliminate the row of d
! modify only a single element of (q transpose)*b
! beyond the first n, which is initially zero.
qtbpj = zero
do k = j, n
! determine a givens rotation which eliminates the
! appropriate element in the current row of d.
if (Sdiag(k) /= zero) then
if (abs(r(k, k)) >= abs(Sdiag(k))) then
tan = Sdiag(k)/r(k, k)
cos = p5/sqrt(p25 + p25*tan**2)
sin = cos*tan
else
cotan = r(k, k)/Sdiag(k)
sin = p5/sqrt(p25 + p25*cotan**2)
cos = sin*cotan
end if
! compute the modified diagonal element of r and
! the modified element of ((q transpose)*b,0).
r(k, k) = cos*r(k, k) + sin*Sdiag(k)
temp = cos*Wa(k) + sin*qtbpj
qtbpj = -sin*Wa(k) + cos*qtbpj
Wa(k) = temp
! accumulate the tranformation in the row of s.
kp1 = k + 1
if (n >= kp1) then
do i = kp1, n
temp = cos*r(i, k) + sin*Sdiag(i)
Sdiag(i) = -sin*r(i, k) + cos*Sdiag(i)
r(i, k) = temp
end do
end if
end if
end do
end if
! store the diagonal element of s and restore
! the corresponding diagonal element of r.
Sdiag(j) = r(j, j)
r(j, j) = x(j)
end do
! solve the triangular system for z. if the system is
! singular, then obtain a least squares solution.
nsing = n
do j = 1, n
if (Sdiag(j) == zero .and. nsing == n) nsing = j - 1
if (nsing < n) Wa(j) = zero
end do
if (nsing >= 1) then
do k = 1, nsing
j = nsing - k + 1
sum = zero
jp1 = j + 1
if (nsing >= jp1) then
do i = jp1, nsing
sum = sum + r(i, j)*Wa(i)
end do
end if
Wa(j) = (Wa(j) - sum)/Sdiag(j)
end do
end if
! permute the components of z back to components of x.
do j = 1, n
l = Ipvt(j)
x(l) = Wa(j)
end do
end subroutine qrsolv