map_most_common Subroutine

   SUBROUTINE map_most_common(array, most_common_map, nmost_common, &
                              most_common_elements, size_limit, max_val)
      !! Makes a lookup table from the most common elements.
      !!
      !! Lookup table
      !! The lookup table is indexed by the most common array values
      !! (i.e., block sizes).  Its values are the order of their frequency.

      INTEGER, DIMENSION(:), INTENT(IN)                  :: array
         !! Array for which to find the most common elements.
      INTEGER(KIND=int_4), DIMENSION(:), POINTER         :: most_common_map
         !! Ranking of the most common elements in array
      INTEGER, INTENT(IN)                                :: nmost_common
         !! The number of most common elements
      INTEGER, DIMENSION(:), INTENT(OUT)                 :: most_common_elements
         !! The most common elements in array
      INTEGER, INTENT(IN)                                :: size_limit
         !! Limit maximum size to this value
      INTEGER, INTENT(OUT)                               :: max_val

      INTEGER                                            :: i, max_val_l, nmc
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: permutation, size_counts

!   ---------------------------------------------------------------------------

      IF (SIZE(array) .GT. 0) THEN
         max_val = MAXVAL(array)
         max_val_l = MIN(MIN(size_limit, max_val), INT(HUGE(most_common_map)))
      ELSE
         max_val = 0
         max_val_l = 0
      END IF
      ! Count the frequency of all block sizes up to max_val_l.
      ALLOCATE (size_counts(0:max_val_l))
      ALLOCATE (permutation(0:max_val_l))
      size_counts = 0
      permutation = 0
      DO i = 1, SIZE(array)
         ! Counts are decreased to easily get a reverse sort order.
         IF (array(i) .LE. max_val_l) &
            size_counts(array(i)) = size_counts(array(i)) - 1
      END DO
      IF (SIZE(array) .GT. 0) THEN
         CALL sort(size_counts, max_val_l + 1, permutation)
      END IF
      ! Limiting nmc to max_val_l prevents out-of-bounds.
      nmc = MIN(nmost_common, max_val_l)
      ! Determine the biggest block size and allocate the map.
      ALLOCATE (most_common_map(0:max_val_l))
      ! Create the mapping from block size to order.
      most_common_map = nmost_common + 1
      DO i = 1, nmc
         most_common_map(permutation(i - 1) - 1) = i
      END DO
      ! Copy the most common elements
      most_common_elements(:) = 0
      most_common_elements(1:nmc) = permutation(0:nmc - 1) - 1
      DEALLOCATE (size_counts)
      DEALLOCATE (permutation)
   END SUBROUTINE map_most_common