Branch data Line data Source code
1 : : !!****m* ABINIT/m_mergesort
2 : : !!
3 : : !! NAME
4 : : !! m_mergesort
5 : : !!
6 : : !! FUNCTION
7 : : !! Module for sorting integer arrays using merge sorting algorithm
8 : : !!
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2010-2026 ABINIT group (hexu)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public Licence, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : : #include "abi_common.h"
24 : :
25 : :
26 : : module m_mergesort
27 : : !!***
28 : : use defs_basis
29 : : use m_errors
30 : : use m_abicore
31 : : use m_mathfuncs, only: array_lessthan, array_morethan, array_le
32 : : implicit none
33 : : private
34 : : public :: MergeSort2D
35 : : public :: MergeSort
36 : :
37 : : interface MergeSort
38 : : procedure MergeSort_dp
39 : : procedure MergeSort_int
40 : : end interface MergeSort
41 : :
42 : : contains
43 : :
44 : : !======== Merge sort algorithm for 1d array ===============
45 : : ! Modified from https://rosettacode.org/wiki/Sorting_algorithms/Merge_sort#Fortran
46 : : ! The original code is public domain.
47 : :
48 [ # # # # : 0 : subroutine merge_int(A, B, C)
# # ]
49 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
50 : : integer, target, intent(in) :: A(:), B(:)
51 : : integer, target, intent(inout) :: C(:)
52 : : integer :: i, j, k
53 : :
54 [ # # ]: 0 : if (size(A) + size(B) > size(C)) then
55 : 0 : stop (1)
56 : : end if
57 : :
58 : 0 : i = 1; j = 1
59 [ # # ]: 0 : do k = 1, size(C)
60 [ # # # # ]: 0 : if (i <= size(A) .and. j <= size(B)) then
61 [ # # ]: 0 : if (A(i) <= B(j)) then
62 : 0 : C(k) = A(i)
63 : 0 : i = i + 1
64 : : else
65 : 0 : C(k) = B(j)
66 : 0 : j = j + 1
67 : : end if
68 [ # # ]: 0 : else if (i <= size(A)) then
69 : 0 : C(k) = A(i)
70 : 0 : i = i + 1
71 [ # # ]: 0 : else if (j <= size(B)) then
72 : 0 : C(k) = B(j)
73 : 0 : j = j + 1
74 : : end if
75 : : end do
76 : 0 : end subroutine merge_int
77 : :
78 [ # # # # : 0 : subroutine merge_with_order_int(A, B, C, orderA, orderB, orderC)
# # # # #
# # # ]
79 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
80 : : integer, target, intent(in) :: A(:), B(:), orderA(:), orderB(:)
81 : : integer, target, intent(inout) :: C(:), orderC(:)
82 : : integer :: i, j, k
83 : :
84 [ # # ]: 0 : if (size(A) + size(B) > size(C)) then
85 : 0 : stop (1)
86 : : end if
87 : :
88 : 0 : i = 1; j = 1
89 [ # # ]: 0 : do k = 1, size(C)
90 [ # # # # ]: 0 : if (i <= size(A) .and. j <= size(B)) then
91 [ # # ]: 0 : if (A(i) <= B(j)) then
92 : 0 : C(k) = A(i)
93 : 0 : orderC(k) = orderA(i)
94 : 0 : i = i + 1
95 : : else
96 : 0 : C(k) = B(j)
97 : 0 : orderC(k) =orderB(j)
98 : 0 : j = j + 1
99 : : end if
100 [ # # ]: 0 : else if (i <= size(A)) then
101 : 0 : C(k) = A(i)
102 : 0 : orderC(k) = orderA(i)
103 : 0 : i = i + 1
104 [ # # ]: 0 : else if (j <= size(B)) then
105 : 0 : C(k) = B(j)
106 : 0 : orderC(k) = orderB(j)
107 : 0 : j = j + 1
108 : : end if
109 : : end do
110 : 0 : end subroutine merge_with_order_int
111 : :
112 : :
113 : 12286 : subroutine swap_int(x, y)
114 : : integer, intent(inout) :: x, y
115 : : integer :: tmp
116 : 12286 : tmp = x; x = y; y = tmp
117 : 12286 : end subroutine swap_int
118 : :
119 [ # # # # ]: 0 : recursive subroutine MergeSort_no_init_order_int(A, work, order, worder)
120 : : integer, intent(inout) :: A(:)
121 : : integer, intent(inout) :: work(:)
122 : : integer, optional, intent(inout):: order(size(A)), worder(size(work))
123 : : integer :: half, n
124 : : logical :: ordered
125 : 0 : n=size(A)
126 : 0 : ordered=.False.
127 [ # # ]: 0 : if (present(order)) then
128 : 0 : ordered=.True.
129 : : end if
130 : :
131 : 0 : half = (size(A) + 1) / 2
132 [ # # ]: 0 : if (size(A) < 2) then
133 : : continue
134 [ # # ]: 0 : else if (size(A) == 2) then
135 [ # # ]: 0 : if (A(1) > A(2)) then
136 : 0 : call swap_int(A(1), A(2))
137 [ # # ]: 0 : if(ordered) call swap_int(order(1), order(2))
138 : : end if
139 : : else
140 [ # # ]: 0 : if(ordered) then
141 [ # # ]: 0 : call MergeSort_no_init_order_int(A( : half), work, order(:half), worder)
142 [ # # ]: 0 : call MergeSort_no_init_order_int(A(half + 1 :), work, order(half+1:), worder)
143 [ # # ]: 0 : if (A(half) > A(half + 1)) then
144 [ # # ]: 0 : work(1 : half) = A(1 : half)
145 [ # # ]: 0 : worder(1:half) = order(1: half)
146 : : call merge_with_order_int(work(1 : half), A(half + 1:), A &
147 [ # # ]: 0 : &, worder(1 : half), order(half + 1:), order)
148 : : endif
149 : : else
150 : 0 : call MergeSort_no_init_order_int(A( : half), work)
151 : 0 : call MergeSort_no_init_order_int(A(half + 1 :), work)
152 [ # # ]: 0 : if (A(half) > A(half + 1)) then
153 [ # # ]: 0 : work(1 : half) = A(1 : half)
154 : 0 : call merge_int(work(1 : half), A(half + 1:), A)
155 : : endif
156 : : endif
157 : : end if
158 [ # # # # ]: 0 : end subroutine MergeSort_No_Init_Order_Int
159 : :
160 [ # # # # ]: 0 : subroutine MergeSort_int(A, work, order, worder)
161 : : integer, intent(inout) :: A(:)
162 : : integer, intent(inout) :: work(:)
163 : : integer, optional, intent(inout):: order(size(A)), worder(size(work))
164 : : integer :: i
165 [ # # ]: 0 : if (present(order)) then
166 [ # # ]: 0 : do i =1 , size(A)
167 : 0 : order(i) = i
168 : : end do
169 : : end if
170 : 0 : call MergeSort_no_init_order_int(A, work, order, worder)
171 [ # # # # ]: 0 : end subroutine MergeSort_int
172 : :
173 : :
174 [ # # # # : 0 : subroutine merge_dp(A, B, C)
# # ]
175 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
176 : : real(dp), target, intent(in) :: A(:), B(:)
177 : : real(dp), target, intent(inout) :: C(:)
178 : : integer :: i, j, k
179 [ # # ]: 0 : if (size(A) + size(B) > size(C)) then
180 : 0 : stop (1)
181 : : end if
182 : 0 : i = 1; j = 1
183 [ # # ]: 0 : do k = 1, size(C)
184 [ # # # # ]: 0 : if (i <= size(A) .and. j <= size(B)) then
185 [ # # ]: 0 : if (A(i) <= B(j)) then
186 : 0 : C(k) = A(i)
187 : 0 : i = i + 1
188 : : else
189 : 0 : C(k) = B(j)
190 : 0 : j = j + 1
191 : : end if
192 [ # # ]: 0 : else if (i <= size(A)) then
193 : 0 : C(k) = A(i)
194 : 0 : i = i + 1
195 [ # # ]: 0 : else if (j <= size(B)) then
196 : 0 : C(k) = B(j)
197 : 0 : j = j + 1
198 : : end if
199 : : end do
200 : 0 : end subroutine merge_dp
201 : :
202 : :
203 : :
204 [ - + - + : 18216 : subroutine merge_with_order_dp(A, B, C, orderA, orderB, orderC)
- + - + -
+ - + ]
205 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
206 : : real(dp), target, intent(in) :: A(:), B(:)
207 : : integer, target, intent(in) :: orderA(:), orderB(:)
208 : : real(dp), target, intent(inout) :: C(:)
209 : : integer, target, intent(inout) :: orderC(:)
210 : : integer :: i, j, k
211 : :
212 [ - + ]: 18216 : if (size(A) + size(B) > size(C)) then
213 : 0 : stop (1)
214 : : end if
215 : 18216 : i = 1; j = 1
216 [ + + ]: 305330 : do k = 1, size(C)
217 [ + + + + ]: 305330 : if (i <= size(A) .and. j <= size(B)) then
218 [ + + ]: 247959 : if (A(i) <= B(j)) then
219 : 129364 : C(k) = A(i)
220 : 129364 : orderC(k) = orderA(i)
221 : 129364 : i = i + 1
222 : : else
223 : 118595 : C(k) = B(j)
224 : 118595 : orderC(k) =orderB(j)
225 : 118595 : j = j + 1
226 : : end if
227 [ + + ]: 39155 : else if (i <= size(A)) then
228 : 18775 : C(k) = A(i)
229 : 18775 : orderC(k) = orderA(i)
230 : 18775 : i = i + 1
231 [ + - ]: 20380 : else if (j <= size(B)) then
232 : 20380 : C(k) = B(j)
233 : 20380 : orderC(k) = orderB(j)
234 : 20380 : j = j + 1
235 : : end if
236 : : end do
237 : 18216 : end subroutine merge_with_order_dp
238 : :
239 : :
240 : :
241 : 8944 : subroutine swap_dp(x, y)
242 : : real(dp), intent(inout) :: x, y
243 : : real(dp) :: tmp
244 : 8944 : tmp = x; x = y; y = tmp
245 : : end subroutine swap_dp
246 : :
247 : :
248 [ - + - + ]: 130332 : recursive subroutine MergeSort_no_init_order_dp(A, work, order, worder)
249 : : real(dp), intent(inout) :: A(:)
250 : : real(dp), intent(inout) :: work(:)
251 : : integer, optional, intent(inout):: order(size(A)), worder(size(work))
252 : : integer :: half, n
253 : : logical :: ordered
254 : 43444 : n=size(A)
255 : 43444 : ordered=.False.
256 [ + - ]: 43444 : if (present(order)) then
257 : 43444 : ordered=.True.
258 : : end if
259 : :
260 : 43444 : half = (size(A) + 1) / 2
261 [ + + ]: 43444 : if (size(A) < 2) then
262 : : continue
263 [ + + ]: 39972 : else if (size(A) == 2) then
264 [ + + ]: 18336 : if (A(1) > A(2)) then
265 : 8944 : call swap_dp(A(1), A(2))
266 [ + - ]: 8944 : if(ordered) call swap_int(order(1), order(2))
267 : : end if
268 : : else
269 [ + - ]: 21636 : if(ordered) then
270 [ + - ]: 21636 : call MergeSort_no_init_order_dp(A( : half), work, order(:half), worder)
271 [ + - ]: 21636 : call MergeSort_no_init_order_dp(A(half + 1 :), work, order(half+1:), worder)
272 [ + + ]: 21636 : if (A(half) > A(half + 1)) then
273 [ + + ]: 166355 : work(1 : half) = A(1 : half)
274 [ + + ]: 166355 : worder(1:half) = order(1: half)
275 : : call merge_with_order_dp(work(1 : half), A(half + 1:), A &
276 [ + - ]: 18216 : &, worder(1 : half), order(half + 1:), order)
277 : : endif
278 : : else
279 : 0 : call MergeSort_no_init_order_dp(A( : half), work)
280 : 0 : call MergeSort_no_init_order_dp(A(half + 1 :), work)
281 [ # # ]: 0 : if (A(half) > A(half + 1)) then
282 [ # # ]: 0 : work(1 : half) = A(1 : half)
283 : 0 : call merge_dp(work(1 : half), A(half + 1:), A)
284 : : endif
285 : : endif
286 : : end if
287 [ + - + - ]: 86888 : end subroutine MergeSort_No_Init_Order_Dp
288 : :
289 [ - + - + ]: 516 : subroutine MergeSort_dp(A, work, order, worder)
290 : : real(dp), intent(inout) :: A(:)
291 : : real(dp), intent(inout) :: work(:)
292 : : integer, optional, intent(inout):: order(size(A)), worder(size(work))
293 : : integer :: i
294 [ + - ]: 172 : if (present(order)) then
295 [ + + ]: 40316 : do i =1 , size(A)
296 : 40316 : order(i) = i
297 : : end do
298 : : end if
299 : 172 : call MergeSort_no_init_order_dp(A, work, order, worder)
300 [ + - + - ]: 344 : end subroutine MergeSort_dp
301 : :
302 : :
303 : :
304 : : !== Merge sort for 2D integer array=================================
305 : :
306 : :
307 [ # # # # : 0 : subroutine merge2D(A, B, C)
# # ]
308 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
309 : : integer, target, intent(in) :: A(:, :), B(:, :)
310 : : integer, target, intent(inout) :: C(:,:)
311 : : integer :: i, j, k
312 : :
313 [ # # ]: 0 : if (size(A, 2) + size(B, 2) > size(C, 2)) then
314 : 0 : stop (1)
315 : : end if
316 : :
317 : 0 : i = 1; j = 1
318 [ # # ]: 0 : do k = 1, size(C, 2)
319 [ # # # # ]: 0 : if (i <= size(A, 2) .and. j <= size(B,2)) then
320 [ # # ]: 0 : if (array_le (A(:, i) , B(:,j), size(A, 1))) then
321 [ # # # # : 0 : C(:,k) = A(:,i)
# # ]
322 : 0 : i = i + 1
323 : : else
324 [ # # # # : 0 : C(:,k) = B(:, j)
# # ]
325 : 0 : j = j + 1
326 : : end if
327 [ # # ]: 0 : else if (i <= size(A, 2)) then
328 [ # # # # : 0 : C(:,k) = A(:, i)
# # ]
329 : 0 : i = i + 1
330 [ # # ]: 0 : else if (j <= size(B, 2)) then
331 [ # # # # : 0 : C(:,k) = B(:,j)
# # ]
332 : 0 : j = j + 1
333 : : end if
334 : : end do
335 : 0 : end subroutine merge2D
336 : :
337 [ - + - + : 46137 : subroutine merge2D_with_order(A, B, C, orderA, orderB, orderC)
- + - + -
+ - + ]
338 : : ! The targe attribute is necessary, because A .or. B might overlap with C.
339 : : integer, target, intent(in) :: A(:, :), B(:,:), orderA(:), orderB(:)
340 : : integer, target, intent(inout) :: C(:,:), orderC(:)
341 : : integer :: i, j, k
342 : :
343 [ - + ]: 46137 : if (size(A, 2) + size(B,2) > size(C,2)) then
344 : 0 : stop (1)
345 : : end if
346 : :
347 : 15999535 : i = 1; j = 1
348 [ + + ]: 15999535 : do k = 1, size(C, 2)
349 [ + + + + ]: 15999535 : if (i <= size(A, 2) .and. j <= size(B, 2)) then
350 [ + + ]: 15325708 : if (array_le (A(:, i) , B(:,j), size(A, 1))) then
351 [ + - + + : 54500663 : C(:,k) = A(:,i)
+ + ]
352 : 7689063 : orderC(k) = orderA(i)
353 : 7689063 : i = i + 1
354 : : else
355 [ + - + + : 54220061 : C(:,k) = B(:,j)
+ + ]
356 : 7636645 : orderC(k) =orderB(j)
357 : 15325708 : j = j + 1
358 : : end if
359 [ + + ]: 627690 : else if (i <= size(A,2)) then
360 [ + - + + : 2104029 : C(:,k) = A(:,i)
+ + ]
361 : 296613 : orderC(k) = orderA(i)
362 : 296613 : i = i + 1
363 [ + - ]: 331077 : else if (j <= size(B, 2)) then
364 [ + - + + : 2255633 : C(:,k) = B(:,j)
+ + ]
365 : 331077 : orderC(k) = orderB(j)
366 : 331077 : j = j + 1
367 : : end if
368 : : end do
369 : 46137 : end subroutine merge2D_with_order
370 : :
371 : :
372 [ - + - + ]: 3342 : subroutine swap2D(x, y)
373 : : integer, intent(inout) :: x(:), y(:)
374 : : integer :: tmp, i
375 [ + + ]: 13804 : do i =1, size(x)
376 : 13804 : tmp = x(i); x(i) = y(i); y(i) = tmp
377 : : end do
378 : 3342 : end subroutine swap2D
379 : :
380 [ - + - + ]: 2607704 : recursive subroutine MergeSort2D_no_init_order(A, work, order, worder)
381 : : integer, intent(inout) :: A(:, :)
382 : : integer, intent(inout) :: work(:, :)
383 : : integer, optional, intent(inout):: order(size(A, 2)), worder(size(work, 2))
384 : : integer :: half, n
385 : : logical :: ordered
386 : 1303852 : n=size(A,2)
387 : 1303852 : ordered=.False.
388 [ + - ]: 1303852 : if (present(order)) then
389 : 1303852 : ordered=.True.
390 : : end if
391 : :
392 : 1303852 : half = (n + 1) / 2
393 [ + + ]: 1303852 : if (n < 2) then
394 : : continue
395 [ + + ]: 1278480 : else if (n == 2) then
396 [ + + ]: 626560 : if (array_morethan(A(:, 1) , A(:,2), size(A, 1))) then
397 : 3342 : call swap2D(A(:,1), A(:,2))
398 [ + - ]: 3342 : if(ordered) call swap_int(order(1), order(2))
399 : : end if
400 : : else
401 [ + - ]: 651920 : if(ordered) then
402 [ + - ]: 651920 : call MergeSort2D_no_init_order(A(:, : half), work, order(:half), worder)
403 [ + - ]: 651920 : call MergeSort2D_no_init_order(A(:, half + 1 :), work, order(half+1:), worder)
404 [ + + ]: 651920 : if (array_morethan(A(:,half) , A(:,half + 1), size(A, 1))) then
405 [ + + + + ]: 32341321 : work(:, 1 : half) = A(:, 1 : half)
406 [ + + ]: 8031813 : worder(1:half) = order(1: half)
407 : : call merge2D_with_order(work(:, 1 : half), A(:, half + 1:), A &
408 [ + - ]: 46137 : &, worder(1 : half), order(half + 1:), order)
409 : : endif
410 : : else
411 : 0 : call MergeSort2D_no_init_order(A(:, : half), work)
412 : 0 : call MergeSort2D_no_init_order(A(:, half + 1 :), work)
413 [ # # ]: 0 : if (array_morethan(A(:,half) , A(:,half + 1), size(A, 1))) then
414 [ # # # # ]: 0 : work(:, 1 : half) = A(:,1 : half)
415 : 0 : call merge2D(work(:, 1 : half), A(:, half + 1:), A)
416 : : endif
417 : : endif
418 : : end if
419 [ + - ]: 1303852 : end subroutine MergeSort2D_No_Init_Order
420 : :
421 [ - + - + ]: 36 : subroutine MergeSort2D(A, work, order, worder)
422 : : integer, intent(inout) :: A(:,:)
423 : : integer, intent(inout) :: work(:, :)
424 : : integer, optional, intent(inout):: order(size(A, 2)), worder(size(work, 2))
425 : :
426 : : integer :: i
427 [ + - ]: 12 : if(present(order)) then
428 [ + + ]: 1278504 : do i =1 , size(A, 2)
429 : 1278504 : order(i) = i
430 : : end do
431 : : end if
432 : 12 : call MergeSort2D_no_init_order(A, work, order, worder)
433 [ + - + - ]: 24 : end subroutine MergeSort2D
434 : :
435 : :
436 : :
437 : : end module m_mergesort
|
