lu6U Subroutine

  subroutine lu6U  ( 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)
    real(rp),      intent(in)    :: v(m)

    integer(ip),   intent(inout) :: luparm(30)
    real(rp),      intent(inout) :: parmlu(30)

    integer(ip),   intent(out)   :: inform
    real(rp),      intent(out)   :: w(n)

    !------------------------------------------------------------------
    ! lu6U   solves   U w = v.          v  is not altered.
    !
    ! 15 Dec 2002: First version derived from lu6sol.
    ! 15 Dec 2002: Current version.
    ! 13 Dec 2011: First f90 version.
    !------------------------------------------------------------------

    integer(ip)            :: i, j, k, klast, l, l1, l2, l3, nrank, nrank1
    real(rp)               :: resid, small, t


    nrank  = luparm(16)
    small  = parmlu(3)
    inform = 0
    nrank1 = nrank + 1
    resid  = zero

    ! Find the first nonzero in v(1:nrank), counting backwards.

    do klast = nrank, 1, -1
       i     = p(klast)
       if (abs(v(i)) > small) exit
    end do

    do k = klast + 1, n
       j    = q(k)
       w(j) = zero
    end do

    ! Do the back-substitution, using rows 1:klast of U.

    do k  = klast, 1, -1
       i  = p(k)
       t  = v(i)
       l1 = locr(i)
       l2 = l1 + 1
       l3 = l1 + lenr(i) - 1

       !***** This loop could be coded specially.
       do l = l2, l3
          j = indr(l)
          t = t - a(l)*w(j)
       end do

       j  = q(k)
       if (abs(t) <= small) then
          w(j) = zero
       else
          w(j) = t/a(l1)
       end if
    end do

    ! Compute residual for overdetermined systems.

    do k = nrank1, m
       i     = p(k)
       resid = resid + abs(v(i))
    end do

    ! Exit.

    if (resid > zero) inform = 1
    luparm(10) = inform
    parmlu(20) = resid

  end subroutine lu6U