HeapSortByKey Subroutine

  subroutine HeapSortByKey(keys,vals)
    !! Sort keys(1:n) ascending in place, permuting vals(1:n) alongside. Heap sort:
    !! O(n log n), no extra storage. Keys are global side ids, distinct within one
    !! neighbour's face list, so stability does not matter.
    implicit none
    integer,intent(inout) :: keys(1:)
    integer,intent(inout) :: vals(1:)
    ! Local
    integer :: n,i,last

    n = size(keys)
    do i = n/2,1,-1
      call SiftDown(i,n)
    enddo
    do last = n,2,-1
      call SwapEntries(1,last)
      call SiftDown(1,last-1)
    enddo

  contains

    subroutine SiftDown(first,last)
      implicit none
      integer,intent(in) :: first,last
      integer :: root,child

      root = first
      do while(2*root <= last)
        child = 2*root
        if(child < last) then
          if(keys(child+1) > keys(child)) child = child+1
        endif
        if(keys(root) >= keys(child)) return
        call SwapEntries(root,child)
        root = child
      enddo

    endsubroutine SiftDown

    subroutine SwapEntries(i1,i2)
      implicit none
      integer,intent(in) :: i1,i2
      integer :: t

      t = keys(i1); keys(i1) = keys(i2); keys(i2) = t
      t = vals(i1); vals(i1) = vals(i2); vals(i2) = t

    endsubroutine SwapEntries

  endsubroutine HeapSortByKey