subroutine lu6Ut ( inform, m, n, v, w, &
lena, luparm, parmlu, a, indr, p, q, lenr, locr )
integer(ip), intent(in) :: m, n, lena
integer(ip), intent(in) :: indr(lena), p(m), q(n), lenr(m), locr(m)
real(rp), intent(in) :: a(lena)
integer(ip), intent(inout) :: luparm(30)
real(rp), intent(inout) :: parmlu(30), w(n)
integer(ip), intent(out) :: inform
real(rp), intent(out) :: v(m)
!------------------------------------------------------------------
! lu6Ut solves U'v = w. w is destroyed.
!
! 15 Dec 2002: First version derived from lu6sol.
! 15 Dec 2002: Current version.
! 13 Dec 2011: First f90 version.
!------------------------------------------------------------------
integer(ip) :: i, j, k, l, l1, l2, nrank, nrank1
real(rp) :: resid, small, t
nrank = luparm(16)
small = parmlu(3)
inform = 0
nrank1 = nrank + 1
resid = zero
do k = nrank1, m
i = p(k)
v(i) = zero
end do
! Do the forward-substitution, skipping columns of U(transpose)
! when the associated element of w(*) is negligible.
do k = 1, nrank
i = p(k)
j = q(k)
t = w(j)
if (abs(t) <= small) then
v(i) = zero
cycle
end if
l1 = locr(i)
t = t/a(l1)
v(i) = t
l2 = l1 + lenr(i) - 1
l1 = l1 + 1
!***** This loop could be coded specially.
do l = l1, l2
j = indr(l)
w(j) = w(j) - t*a(l)
end do
end do
! Compute residual for overdetermined systems.
do k = nrank1, n
j = q(k)
resid = resid + abs(w(j))
end do
! Exit.
if (resid > zero) inform = 1
luparm(10) = inform
parmlu(20) = resid
end subroutine lu6Ut