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
|
