Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE cp_min_heap
9 : USE kinds, ONLY: int_4,&
10 : int_8
11 : #include "../base/base_uses.f90"
12 :
13 : IMPLICIT NONE
14 : PRIVATE
15 :
16 : PUBLIC :: cp_heap_type, keyt, valt
17 : PUBLIC :: cp_heap_pop, cp_heap_reset_node, cp_heap_fill
18 : PUBLIC :: cp_heap_new, cp_heap_release
19 : PUBLIC :: cp_heap_get_first, cp_heap_reset_first
20 :
21 : ! Sets the types
22 : INTEGER, PARAMETER :: keyt = int_4
23 : INTEGER, PARAMETER :: valt = int_8
24 :
25 : TYPE cp_heap_node
26 : INTEGER(KIND=keyt) :: key = -1_keyt
27 : INTEGER(KIND=valt) :: value = -1_valt
28 : END TYPE cp_heap_node
29 :
30 : TYPE cp_heap_node_e
31 : TYPE(cp_heap_node) :: node = cp_heap_node()
32 : END TYPE cp_heap_node_e
33 :
34 : TYPE cp_heap_type
35 : INTEGER :: n = -1
36 : INTEGER, DIMENSION(:), POINTER :: index => NULL()
37 : TYPE(cp_heap_node_e), DIMENSION(:), POINTER :: nodes => NULL()
38 : END TYPE cp_heap_type
39 :
40 : CONTAINS
41 :
42 : ! Lookup functions
43 :
44 : ! **************************************************************************************************
45 : !> \brief ...
46 : !> \param n ...
47 : !> \return ...
48 : ! **************************************************************************************************
49 1323895 : ELEMENTAL FUNCTION get_parent(n) RESULT(parent)
50 : INTEGER, INTENT(IN) :: n
51 : INTEGER :: parent
52 :
53 1323895 : parent = INT(n/2)
54 1323895 : END FUNCTION get_parent
55 :
56 : ! **************************************************************************************************
57 : !> \brief ...
58 : !> \param n ...
59 : !> \return ...
60 : ! **************************************************************************************************
61 615323 : ELEMENTAL FUNCTION get_left_child(n) RESULT(child)
62 : INTEGER, INTENT(IN) :: n
63 : INTEGER :: child
64 :
65 615323 : child = 2*n
66 615323 : END FUNCTION get_left_child
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param n ...
71 : !> \return ...
72 : ! **************************************************************************************************
73 0 : ELEMENTAL FUNCTION get_right_child(n) RESULT(child)
74 : INTEGER, INTENT(IN) :: n
75 : INTEGER :: child
76 :
77 0 : child = 2*n + 1
78 0 : END FUNCTION get_right_child
79 :
80 : ! **************************************************************************************************
81 : !> \brief ...
82 : !> \param heap ...
83 : !> \param n ...
84 : !> \return ...
85 : ! **************************************************************************************************
86 1230685 : ELEMENTAL FUNCTION get_value(heap, n) RESULT(value)
87 : TYPE(cp_heap_type), INTENT(IN) :: heap
88 : INTEGER, INTENT(IN) :: n
89 : INTEGER(KIND=valt) :: value
90 :
91 1230685 : value = heap%nodes(n)%node%value
92 1230685 : END FUNCTION get_value
93 :
94 : ! **************************************************************************************************
95 : !> \brief ...
96 : !> \param heap ...
97 : !> \param key ...
98 : !> \return ...
99 : ! **************************************************************************************************
100 0 : ELEMENTAL FUNCTION get_value_by_key(heap, key) RESULT(value)
101 : TYPE(cp_heap_type), INTENT(IN) :: heap
102 : INTEGER(KIND=keyt), INTENT(IN) :: key
103 : INTEGER(KIND=valt) :: value
104 :
105 : INTEGER :: n
106 :
107 0 : n = heap%index(key)
108 0 : value = get_value(heap, n)
109 0 : END FUNCTION get_value_by_key
110 :
111 : ! Initialization functions
112 :
113 : ! **************************************************************************************************
114 : !> \brief ...
115 : !> \param heap ...
116 : !> \param n ...
117 : ! **************************************************************************************************
118 30445 : SUBROUTINE cp_heap_new(heap, n)
119 : TYPE(cp_heap_type), INTENT(OUT) :: heap
120 : INTEGER, INTENT(IN) :: n
121 :
122 30445 : heap%n = n
123 91335 : ALLOCATE (heap%index(n))
124 149434 : ALLOCATE (heap%nodes(n))
125 30445 : END SUBROUTINE cp_heap_new
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param heap ...
130 : ! **************************************************************************************************
131 30445 : SUBROUTINE cp_heap_release(heap)
132 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
133 :
134 30445 : DEALLOCATE (heap%index)
135 30445 : DEALLOCATE (heap%nodes)
136 30445 : heap%n = 0
137 30445 : END SUBROUTINE cp_heap_release
138 :
139 : ! **************************************************************************************************
140 : !> \brief Fill heap with given values
141 : !> \param heap ...
142 : !> \param values ...
143 : ! **************************************************************************************************
144 30445 : SUBROUTINE cp_heap_fill(heap, values)
145 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
146 : INTEGER(KIND=valt), DIMENSION(:), INTENT(IN) :: values
147 :
148 : INTEGER :: first, i, n
149 :
150 30445 : n = SIZE(values)
151 30445 : CPASSERT(heap%n >= n)
152 :
153 88544 : DO i = 1, n
154 58099 : heap%index(i) = i
155 58099 : heap%nodes(i)%node%key = i
156 88544 : heap%nodes(i)%node%value = values(i)
157 : END DO
158 : ! Sort from the last full subtree
159 30445 : first = get_parent(n)
160 58088 : DO i = first, 1, -1
161 58088 : CALL bubble_down(heap, i)
162 : END DO
163 :
164 30445 : END SUBROUTINE cp_heap_fill
165 :
166 : ! **************************************************************************************************
167 : !> \brief Returns the first heap element without removing it.
168 : !> \param heap ...
169 : !> \param key ...
170 : !> \param value ...
171 : !> \param found ...
172 : ! **************************************************************************************************
173 650456 : SUBROUTINE cp_heap_get_first(heap, key, value, found)
174 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
175 : INTEGER(KIND=keyt), INTENT(OUT) :: key
176 : INTEGER(KIND=valt), INTENT(OUT) :: value
177 : LOGICAL, INTENT(OUT) :: found
178 :
179 650456 : IF (heap%n .LT. 1) THEN
180 0 : found = .FALSE.
181 : ELSE
182 650456 : found = .TRUE.
183 650456 : key = heap%nodes(1)%node%key
184 650456 : value = heap%nodes(1)%node%value
185 : END IF
186 650456 : END SUBROUTINE cp_heap_get_first
187 :
188 : ! **************************************************************************************************
189 : !> \brief Returns and removes the first heap element and rebalances
190 : !> the heap.
191 : !> \param heap ...
192 : !> \param key ...
193 : !> \param value ...
194 : !> \param found ...
195 : ! **************************************************************************************************
196 85 : SUBROUTINE cp_heap_pop(heap, key, value, found)
197 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
198 : INTEGER(KIND=keyt), INTENT(OUT) :: key
199 : INTEGER(KIND=valt), INTENT(OUT) :: value
200 : LOGICAL, INTENT(OUT) :: found
201 :
202 : !
203 :
204 85 : CALL cp_heap_get_first(heap, key, value, found)
205 85 : IF (found) THEN
206 85 : IF (heap%n .GT. 1) THEN
207 48 : CALL cp_heap_copy_node(heap, 1, heap%n)
208 48 : heap%n = heap%n - 1
209 48 : CALL bubble_down(heap, 1)
210 : ELSE
211 37 : heap%n = heap%n - 1
212 : END IF
213 : END IF
214 85 : END SUBROUTINE cp_heap_pop
215 :
216 : ! **************************************************************************************************
217 : !> \brief Changes the value of the heap element with given key and
218 : !> rebalances the heap.
219 : !> \param heap ...
220 : !> \param key ...
221 : !> \param value ...
222 : ! **************************************************************************************************
223 102 : SUBROUTINE cp_heap_reset_node(heap, key, value)
224 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
225 : INTEGER(KIND=keyt), INTENT(IN) :: key
226 : INTEGER(KIND=valt), INTENT(IN) :: value
227 :
228 : INTEGER :: n, new_pos
229 :
230 51 : CPASSERT(heap%n > 0)
231 :
232 51 : n = heap%index(key)
233 51 : CPASSERT(heap%nodes(n)%node%key == key)
234 51 : heap%nodes(n)%node%value = value
235 51 : CALL bubble_up(heap, n, new_pos)
236 51 : CALL bubble_down(heap, new_pos)
237 51 : END SUBROUTINE cp_heap_reset_node
238 :
239 : ! **************************************************************************************************
240 : !> \brief Changes the value of the minimum heap element and rebalances the heap.
241 : !> \param heap ...
242 : !> \param value ...
243 : ! **************************************************************************************************
244 650371 : SUBROUTINE cp_heap_reset_first(heap, value)
245 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
246 : INTEGER(KIND=valt), INTENT(IN) :: value
247 :
248 650371 : CPASSERT(heap%n > 0)
249 650371 : heap%nodes(1)%node%value = value
250 650371 : CALL bubble_down(heap, 1)
251 650371 : END SUBROUTINE cp_heap_reset_first
252 :
253 : ! **************************************************************************************************
254 : !> \brief ...
255 : !> \param heap ...
256 : !> \param e1 ...
257 : !> \param e2 ...
258 : ! **************************************************************************************************
259 615330 : PURE SUBROUTINE cp_heap_swap(heap, e1, e2)
260 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
261 : INTEGER, INTENT(IN) :: e1, e2
262 :
263 : INTEGER(KIND=keyt) :: key1, key2
264 : TYPE(cp_heap_node) :: tmp_node
265 :
266 615330 : key1 = heap%nodes(e1)%node%key
267 615330 : key2 = heap%nodes(e2)%node%key
268 :
269 615330 : tmp_node = heap%nodes(e1)%node
270 615330 : heap%nodes(e1)%node = heap%nodes(e2)%node
271 615330 : heap%nodes(e2)%node = tmp_node
272 :
273 615330 : heap%index(key1) = e2
274 615330 : heap%index(key2) = e1
275 615330 : END SUBROUTINE cp_heap_swap
276 :
277 : ! **************************************************************************************************
278 : !> \brief Sets node e1 to e2
279 : !> \param heap ...
280 : !> \param e1 ...
281 : !> \param e2 ...
282 : ! **************************************************************************************************
283 48 : PURE SUBROUTINE cp_heap_copy_node(heap, e1, e2)
284 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
285 : INTEGER, INTENT(IN) :: e1, e2
286 :
287 : INTEGER(KIND=keyt) :: key1, key2
288 :
289 48 : key1 = heap%nodes(e1)%node%key
290 48 : key2 = heap%nodes(e2)%node%key
291 :
292 48 : heap%nodes(e1)%node = heap%nodes(e2)%node
293 :
294 48 : heap%index(key1) = 0
295 48 : heap%index(key2) = e1
296 48 : END SUBROUTINE cp_heap_copy_node
297 :
298 : ! **************************************************************************************************
299 : !> \brief Balances a heap by bubbling down from the given element.
300 : !> \param heap ...
301 : !> \param first ...
302 : ! **************************************************************************************************
303 678113 : SUBROUTINE bubble_down(heap, first)
304 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
305 : INTEGER, INTENT(IN) :: first
306 :
307 : INTEGER :: e, left_child, right_child, smallest
308 : INTEGER(kind=valt) :: left_child_value, min_value, &
309 : right_child_value
310 : LOGICAL :: all_done
311 :
312 : !
313 678113 : CPASSERT(0 < first .AND. first <= heap%n)
314 :
315 678113 : e = first
316 678113 : all_done = .FALSE.
317 : ! Check whether we are finished, i.e,. whether the element to
318 : ! bubble down is childless.
319 1293436 : DO WHILE (e .LE. get_parent(heap%n) .AND. .NOT. all_done)
320 : ! Determines which node (current, left, or right child) has the
321 : ! smallest value.
322 615323 : smallest = e
323 615323 : min_value = get_value(heap, e)
324 615323 : left_child = get_left_child(e)
325 615323 : IF (left_child .LE. heap%n) THEN
326 615323 : left_child_value = get_value(heap, left_child)
327 615323 : IF (left_child_value .LT. min_value) THEN
328 365138 : min_value = left_child_value
329 365138 : smallest = left_child
330 : END IF
331 : END IF
332 615323 : right_child = left_child + 1
333 615323 : IF (right_child .LE. heap%n) THEN
334 11 : right_child_value = get_value(heap, right_child)
335 11 : IF (right_child_value .LT. min_value) THEN
336 0 : min_value = right_child_value
337 0 : smallest = right_child
338 : END IF
339 : END IF
340 : !
341 615323 : CALL cp_heap_swap(heap, e, smallest)
342 1293436 : IF (smallest .EQ. e) THEN
343 : all_done = .TRUE.
344 : ELSE
345 365138 : e = smallest
346 : END IF
347 : END DO
348 678113 : END SUBROUTINE bubble_down
349 :
350 : ! **************************************************************************************************
351 : !> \brief Balances a heap by bubbling up from the given element.
352 : !> \param heap ...
353 : !> \param first ...
354 : !> \param new_pos ...
355 : ! **************************************************************************************************
356 51 : SUBROUTINE bubble_up(heap, first, new_pos)
357 : TYPE(cp_heap_type), INTENT(INOUT) :: heap
358 : INTEGER, INTENT(IN) :: first
359 : INTEGER, INTENT(OUT) :: new_pos
360 :
361 : INTEGER :: e, parent
362 : INTEGER(kind=valt) :: my_value, parent_value
363 : LOGICAL :: all_done
364 :
365 51 : CPASSERT(0 < first .AND. first <= heap%n)
366 :
367 51 : e = first
368 51 : all_done = .FALSE.
369 51 : IF (e .GT. 1) THEN
370 14 : my_value = get_value(heap, e)
371 : END IF
372 : ! Check whether we are finished, i.e,. whether the element to
373 : ! bubble up is an orphan.
374 51 : new_pos = e
375 65 : DO WHILE (e .GT. 1 .AND. .NOT. all_done)
376 : ! Switches the parent and the current element if the current
377 : ! element's value is greater than the parent's value.
378 14 : parent = get_parent(e)
379 14 : parent_value = get_value(heap, parent)
380 65 : IF (my_value .LT. parent_value) THEN
381 7 : CALL cp_heap_swap(heap, e, parent)
382 7 : e = parent
383 : ELSE
384 : all_done = .TRUE.
385 : END IF
386 : END DO
387 51 : new_pos = e
388 51 : END SUBROUTINE bubble_up
389 :
390 0 : END MODULE cp_min_heap
|
