WPS/metgrid/src/minheap_module.f90

   1 ! Implements a heap using an array; top of the heap is the item
   2 !   with minimum key value
   3 module minheap_module
   4
   5    use datatype_module
   6    use module_debug
   7
   8    ! Maximum heap size -- maybe make this magically dynamic somehow?
   9    integer, parameter :: HEAPSIZE = 10000
  10
  11    ! Type of item to be stored in the heap
  12    type heap_object
  13       type (data_node), pointer :: object
  14    end type heap_object
  15
  16    ! The heap itself
  17    type (heap_object), allocatable, dimension(:) :: heap
  18
  19    ! Index of last item in the heap
  20    integer :: end_of_heap
  21
  22    contains
  23
  24
  25    ! Initialize the heap; current functionality can be had without
  26    !   the need for init function, but we may want more things later
  27    subroutine init_heap()
  28
  29       implicit none
  30
  31       end_of_heap = 0
  32       allocate(heap(HEAPSIZE))
  33
  34    end subroutine init_heap
  35
  36
  37    subroutine heap_destroy()
  38
  39       implicit none
  40
  41       deallocate(heap)
  42
  43    end subroutine heap_destroy
  44
  45
  46    subroutine add_to_heap(x)
  47
  48       implicit none
  49
  50       ! Arguments
  51       type (data_node), pointer :: x
  52
  53       ! Local variables
  54       integer :: idx, parent
  55
  56       call mprintf((end_of_heap == HEAPSIZE),ERROR, 'add_to_heap(): Maximum heap size exceeded')
  57
  58       end_of_heap = end_of_heap + 1
  59       idx = end_of_heap
  60       heap(idx)%object => x
  61       heap(idx)%object%heap_index = idx
  62
  63       do while (idx > 1)
  64          parent = floor(real(idx)/2.)
  65          if (heap(idx)%object%last_used < heap(parent)%object%last_used) then
  66             heap(idx)%object => heap(parent)%object
  67             heap(idx)%object%heap_index = idx
  68             heap(parent)%object => x
  69             heap(parent)%object%heap_index = parent
  70             idx = parent
  71          else
  72             idx = 1
  73          end if
  74       end do
  75
  76    end subroutine add_to_heap
  77
  78
  79    subroutine remove_index(idx)
  80
  81       implicit none
  82
  83       ! Arguments
  84       integer, intent(in) :: idx
  85
  86       ! Local variables
  87       integer :: indx, left, right
  88       type (data_node), pointer :: temp
  89
  90       heap(idx)%object => heap(end_of_heap)%object
  91       heap(idx)%object%heap_index = idx
  92       end_of_heap = end_of_heap - 1
  93
  94       indx = idx
  95
  96       do while (indx <= end_of_heap)
  97          left = indx*2
  98          right = indx*2+1
  99          if (right <= end_of_heap) then
 100             if (heap(right)%object%last_used < heap(left)%object%last_used) then
 101                if (heap(right)%object%last_used < heap(indx)%object%last_used) then
 102                   temp => heap(indx)%object
 103                   heap(indx)%object => heap(right)%object
 104                   heap(indx)%object%heap_index = indx
 105                   heap(right)%object => temp
 106                   heap(right)%object%heap_index = right
 107                   indx = right
 108                else
 109                   indx = end_of_heap + 1
 110                end if
 111             else
 112                if (heap(left)%object%last_used < heap(indx)%object%last_used) then
 113                   temp => heap(indx)%object
 114                   heap(indx)%object => heap(left)%object
 115                   heap(indx)%object%heap_index = indx
 116                   heap(left)%object => temp
 117                   heap(left)%object%heap_index = left
 118                   indx = left
 119                else
 120                   indx = end_of_heap + 1
 121                end if
 122             end if
 123          else if (left <= end_of_heap) then
 124             if (heap(left)%object%last_used < heap(indx)%object%last_used) then
 125                temp => heap(indx)%object
 126                heap(indx)%object => heap(left)%object
 127                heap(indx)%object%heap_index = indx
 128                heap(left)%object => temp
 129                heap(left)%object%heap_index = left
 130                indx = left
 131             else
 132                indx = end_of_heap + 1
 133             end if
 134          else
 135             indx = end_of_heap + 1
 136          end if
 137       end do
 138
 139    end subroutine remove_index
 140
 141
 142    subroutine get_min(x)
 143
 144       implicit none
 145
 146       ! Arguments
 147       type (data_node), pointer :: x
 148
 149       ! Local variables
 150       integer :: idx, left, right
 151       type (data_node), pointer :: temp
 152
 153       call mprintf((end_of_heap <= 0),ERROR, 'get_min(): No items left in the heap.')
 154
 155       x => heap(1)%object
 156
 157       heap(1)%object => heap(end_of_heap)%object
 158       heap(1)%object%heap_index = 1
 159       end_of_heap = end_of_heap - 1
 160       idx = 1
 161
 162       do while (idx <= end_of_heap)
 163          left = idx*2
 164          right = idx*2+1
 165          if (right <= end_of_heap) then
 166             if (heap(right)%object%last_used < heap(left)%object%last_used) then
 167                if (heap(right)%object%last_used < heap(idx)%object%last_used) then
 168                   temp => heap(idx)%object
 169                   heap(idx)%object => heap(right)%object
 170                   heap(idx)%object%heap_index = idx
 171                   heap(right)%object => temp
 172                   heap(right)%object%heap_index = right
 173                   idx = right
 174                else
 175                   idx = end_of_heap + 1
 176                end if
 177             else
 178                if (heap(left)%object%last_used < heap(idx)%object%last_used) then
 179                   temp => heap(idx)%object
 180                   heap(idx)%object => heap(left)%object
 181                   heap(idx)%object%heap_index = idx
 182                   heap(left)%object => temp
 183                   heap(left)%object%heap_index = left
 184                   idx = left
 185                else
 186                   idx = end_of_heap + 1
 187                end if
 188             end if
 189          else if (left <= end_of_heap) then
 190             if (heap(left)%object%last_used < heap(idx)%object%last_used) then
 191                temp => heap(idx)%object
 192                heap(idx)%object => heap(left)%object
 193                heap(idx)%object%heap_index = idx
 194                heap(left)%object => temp
 195                heap(left)%object%heap_index = left
 196                idx = left
 197             else
 198                idx = end_of_heap + 1
 199             end if
 200          else
 201             idx = end_of_heap + 1
 202          end if
 203       end do
 204
 205    end subroutine get_min
 206
 207 end module minheap_module