lu1pq1 Subroutine

  subroutine lu1pq1( m, n, len, iperm, loc, inv, num )

    integer(ip),   intent(in)    :: m, n
    integer(ip),   intent(in)    :: len(m)
    integer(ip),   intent(out)   :: iperm(m), loc(n), inv(m)
    integer(ip),   intent(out)   :: num(n) ! workspace

    !------------------------------------------------------------------
    ! lu1pq1  constructs a permutation  iperm  from the array  len.
    !
    ! On entry:
    ! len(i)  holds the number of nonzeros in the i-th row (say)
    !         of an m by n matrix.
    ! num(*)  can be anything (workspace).
    !
    ! On exit:
    ! iperm   contains a list of row numbers in the order
    !         rows of length 0,  rows of length 1,..., rows of length n.
    ! loc(nz) points to the first row containing  nz  nonzeros,
    !         nz = 1, n.
    ! inv(i)  points to the position of row i within iperm(*).
    !
    ! 10 Jan 2010: First f90 version.
    ! 12 Dec 2011: Declare intent and local variables.
    !------------------------------------------------------------------

    integer(ip)        :: i, l, nz, nzero

    ! Count the number of rows of each length.

    nzero    = 0
    num(1:n) = 0
    loc(1:n) = 0

    do i = 1, m
       nz      = len(i)
       if (nz == 0) then
          nzero   = nzero   + 1
       else
          num(nz) = num(nz) + 1
       end if
    end do

    ! Set starting locations for each length.

    l      = nzero + 1
    do nz  = 1, n
       loc(nz) = l
       l       = l + num(nz)
       num(nz) = 0
    end do

    ! Form the list.

    nzero  = 0
    do i = 1, m
       nz = len(i)
       if (nz == 0) then
          nzero = nzero + 1
          iperm(nzero) = i
       else
          l        = loc(nz) + num(nz)
          iperm(l) = i
          num(nz)  = num(nz) + 1
       end if
    end do

    ! Define the inverse of iperm.

    do l = 1, m
       i      = iperm(l)
       inv(i) = l
    end do

  end subroutine lu1pq1