subroutine lu1or2( n, numa, lena, a, inum, jnum, lenc, locc )
integer(ip), intent(in) :: n, numa, lena
integer(ip), intent(in) :: lenc(n)
integer(ip), intent(inout) :: inum(lena), jnum(lena)
real(rp), intent(inout) :: a(lena)
integer(ip), intent(out) :: locc(n)
!------------------------------------------------------------------
! lu1or2 sorts a list of matrix elements a(i,j) into column
! order, given numa entries a(i,j), i, j in the parallel
! arrays a, inum, jnum respectively. The matrix is assumed
! to have n columns and an arbitrary number of rows.
!
! On entry, lenc(*) must contain the length of each column.
!
! On exit, a(*) and inum(*) are sorted, jnum(*) = 0, and
! locc(j) points to the start of column j.
!
! lu1or2 is derived from mc20ad, a routine in the Harwell
! Subroutine Library, author J. K. Reid.
! xx Feb 1985: Original version.
! 17 Oct 2000: a, inum, jnum now have size lena to allow nelem = 0.
!
! 10 Jan 2010: First f90 version.
! 12 Dec 2011: Declare intent and local variables.
!------------------------------------------------------------------
integer(ip) :: i, ice, icep, j, ja, jb, jce, jcep, l
real(rp) :: ace, acep
! Set loc(j) to point to the beginning of column j.
l = 1
do j = 1, n
locc(j) = l
l = l + lenc(j)
end do
! Sort the elements into column order.
! The algorithm is an in-place sort and is of order numa.
do i = 1, numa
! Establish the current entry.
jce = jnum(i)
if (jce == 0) cycle
ace = a(i)
ice = inum(i)
jnum(i) = 0
! Chain from current entry.
do j = 1, numa
! The current entry is not in the correct position.
! Determine where to store it.
l = locc(jce)
locc(jce) = locc(jce) + 1
! Save the contents of that location.
acep = a(l)
icep = inum(l)
jcep = jnum(l)
! Store current entry.
a(l) = ace
inum(l) = ice
jnum(l) = 0
! If next current entry needs to be processed,
! copy it into current entry.
if (jcep == 0) exit
ace = acep
ice = icep
jce = jcep
end do
end do
! Reset loc(j) to point to the start of column j.
ja = 1
do j = 1, n
jb = locc(j)
locc(j) = ja
ja = jb
end do
end subroutine lu1or2