SCALE-RM
scale_sort.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
11 module scale_sort
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  !-----------------------------------------------------------------------------
19  implicit none
20  private
21  !-----------------------------------------------------------------------------
22  !
23  !++ Public procedure
24  !
25  public :: sort_exec
26  public :: sort_uniq_int_sorted
27  public :: sort_quicksort
28  public :: sort_heapsort
29  public :: sort_quickselect
30  public :: sort_quickselect_arg
31  public :: sort_quickselect_desc
32  public :: sort_quickselect_desc_arg
33 
34  interface sort_exec
35  module procedure sort_exec_without_idx
36  module procedure sort_exec_with_idxs
37  module procedure sort_exec_with_idx
38  end interface sort_exec
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  private :: swap
53  private :: swap_i
54  private :: partition
55  private :: partition_arg
56  private :: partition_desc
57  private :: partition_desc_arg
58  private :: median_of_three
59  private :: median_of_three_arg
60  private :: sample_second_min
61  private :: sample_second_min_arg
62  private :: sample_second_max
63  private :: sample_second_max_arg
64 
65  !-----------------------------------------------------------------------------
66 contains
67  !-----------------------------------------------------------------------------
69 !OCL SERIAL
70  subroutine sort_exec_with_idxs( &
71  npoints, &
72  val, &
73  idx_i, idx_j, &
74  reverse )
75  !$acc routine seq
76  implicit none
77  integer, intent(in) :: npoints ! number of interpolation points
78  real(RP), intent(inout) :: val (npoints) ! value to sort
79  integer, intent(inout) :: idx_i(npoints) ! i-index
80  integer, intent(inout) :: idx_j(npoints) ! j-index
81 
82  logical, intent(in), optional :: reverse
83 
84  real(RP) :: sig
85  integer :: itmp
86  integer :: jtmp
87  real(RP) :: vtmp
88 
89  integer :: n1, n2
90  !---------------------------------------------------------------------------
91 
92  sig = 1.0_rp
93  if ( present(reverse) ) then
94  if ( reverse ) sig = -1.0_rp
95  end if
96 
97  do n1 = 1, npoints-1
98  do n2 = n1+1, npoints
99  if ( val(n1) * sig > val(n2) * sig ) then
100  itmp = idx_i(n1)
101  jtmp = idx_j(n1)
102  vtmp = val(n1)
103 
104  idx_i(n1) = idx_i(n2)
105  idx_j(n1) = idx_j(n2)
106  val(n1) = val(n2)
107 
108  idx_i(n2) = itmp
109  idx_j(n2) = jtmp
110  val(n2) = vtmp
111  endif
112  enddo
113  enddo
114 
115  return
116  end subroutine sort_exec_with_idxs
117 
118 !OCL SERIAL
119  subroutine sort_exec_with_idx( &
120  npoints, &
121  val, index, &
122  reverse )
123  !$acc routine seq
124  implicit none
125  integer, intent(in) :: npoints ! number of interpolation points
126  real(RP), intent(inout) :: val (npoints) ! value to sort
127  integer, intent(inout) :: index(npoints) ! index
128 
129  logical, intent(in), optional :: reverse
130 
131  real(RP) :: sig
132  integer :: itmp
133  real(RP) :: vtmp
134 
135  integer :: n1, n2
136  !---------------------------------------------------------------------------
137 
138  sig = 1.0_rp
139  if ( present(reverse) ) then
140  if ( reverse ) sig = -1.0_rp
141  end if
142 
143  do n1 = 1, npoints-1
144  do n2 = n1+1, npoints
145  if ( val(n1) * sig > val(n2) * sig ) then
146  itmp = index(n1)
147  vtmp = val(n1)
148 
149  index(n1) = index(n2)
150  val(n1) = val(n2)
151 
152  index(n2) = itmp
153  val(n2) = vtmp
154  endif
155  enddo
156  enddo
157 
158  return
159  end subroutine sort_exec_with_idx
160 
161 !OCL SERIAL
162  subroutine sort_exec_without_idx( &
163  npoints, &
164  val, &
165  reverse )
166  !$acc routine seq
167  implicit none
168  integer, intent(in) :: npoints ! number of interpolation points
169  real(RP), intent(inout) :: val (npoints) ! value to sort
170 
171  logical, intent(in), optional :: reverse
172 
173  real(RP) :: sig
174  real(RP) :: vtmp
175 
176  integer :: n1, n2
177  !---------------------------------------------------------------------------
178 
179  sig = 1.0_rp
180  if ( present(reverse) ) then
181  if ( reverse ) sig = -1.0_rp
182  end if
183 
184  do n1 = 1, npoints-1
185  do n2 = n1+1, npoints
186  if ( val(n1) * sig > val(n2) * sig ) then
187  vtmp = val(n1)
188 
189  val(n1) = val(n2)
190 
191  val(n2) = vtmp
192  endif
193  enddo
194  enddo
195 
196  return
197  end subroutine sort_exec_without_idx
198 
199  subroutine sort_uniq_int_sorted(n, ary, c)
200  integer(DP), intent(in) :: n
201  integer(DP), intent(inout) :: ary(n)
202  integer(DP), intent(out) :: c
203  integer(DP) :: i, j
204 
205  c = 1
206  j = ary(1)
207  do i = 2, n
208  if(j .ne. ary(i)) then
209  j = ary(i)
210  c = c + 1
211  ary(c) = ary(i)
212  end if
213  end do
214 
215  return
216  end subroutine sort_uniq_int_sorted
217 
218  recursive subroutine sort_quicksort(n, array)
219  implicit none
220  integer(DP), intent(in) :: n
221  integer(DP), intent(inout) :: array(n)
222  integer(DP) :: i, j, k, k1, kn2, kn, tmp
223  integer, parameter :: threshold = 4
224 
225  if(n <= 1) return
226 
227  k1 = array(1)
228  kn2 = array(n / 2)
229  kn = array(n)
230  if(k1 < kn2) then
231  if(kn2 < kn) then
232  k = kn2 ! 1, n / 2, n
233  else ! array(n) <= array(n / 2)
234  if(k1 < kn) then
235  k = kn ! 1, n, n / 2
236  else
237  k = k1 ! n, 1, n / 2
238  end if
239  end if
240  else ! array(n / 2) <= array(1)
241  if(k1 < kn) then
242  k = k1 ! n / 2, 1, n
243  else ! array(n) <= array(1)
244  if(kn2 < kn) then
245  k = kn ! n / 2, n, 1
246  else ! array(n) <= array(n / 2)
247  k = kn2 ! n, n / 2, 1
248  end if
249  end if
250  end if
251 
252  i = 1
253  j = n
254  do
255  do i = i, n, 1
256  if(array(i) .ge. k) exit
257  end do
258  do j = j, 1, -1
259  if(array(j) .le. k) exit
260  end do
261  if(i >= j) exit
262  tmp = array(i)
263  array(i) = array(j)
264  array(j) = tmp
265  i = i + 1
266  j = j - 1
267  end do
268 
269  if(i - 1 > threshold) then
270  call sort_quicksort(i - 1, array(1:(i - 1)))
271  else if(i - 1 > 1) then
272  call sort_heapsort(i - 1, array(1:(i - 1)))
273  end if
274  if(n - j > threshold) then
275  call sort_quicksort(n - j, array((j + 1):n))
276  else if(n - j > 1) then
277  call sort_heapsort(n - j, array((j + 1):n))
278  end if
279 
280  return
281  end subroutine sort_quicksort
282 
283  subroutine sort_heapsort(n, array)
284  implicit none
285  integer(DP), intent(in) :: n
286  integer(DP), intent(inout) :: array(n)
287 
288  integer(DP) :: i, j, k, l
289  integer(DP) :: t
290 
291  l = n / 2 + 1
292  k = n
293  do while(k .ne. 1)
294  if(l .gt. 1) then
295  l = l - 1
296  t = array(l)
297  else
298  t = array(k)
299  array(k) = array(1)
300  k = k - 1
301  if(k .eq. 1) then
302  array(1) = t
303  exit
304  end if
305  end if
306  i = l
307  j = l + l
308  do while(j .le. k)
309  if(j .lt. k)then
310  if(array(j) .lt. array(j + 1)) j = j + 1
311  end if
312  if(t .lt. array(j))then
313  array(i) = array(j)
314  i = j
315  j = j + j
316  else
317  j = k + 1
318  endif
319  enddo
320  array(i) = t
321  end do
322 
323  return
324  end subroutine sort_heapsort
325 
326  recursive subroutine sort_quickselect(A, left, right, K)
327  implicit none
328  real(rp), intent(inout) :: a(:)
329  integer, intent(in) :: left, right
330  integer, intent(in) :: k
331  integer :: middle, pivot
332 
333  if (left < right) then
334  if ((right-left)/k >= 2) then
335  call sample_second_min(a, left, right, k, pivot)
336  else
337  middle = (left + right) / 2
338  call median_of_three(a, left, middle, right, pivot)
339  end if
340  call partition(a, left, right, pivot)
341  if (k < pivot) then
342  call sort_quickselect(a, left, pivot-1, k)
343  else if (k > pivot) then
344  call sort_quickselect(a, pivot+1, right, k)
345  end if
346  end if
347 
348  return
349  end subroutine sort_quickselect
350 
351  recursive subroutine sort_quickselect_arg(A, X, left, right, K)
352  implicit none
353  real(rp), intent(in) :: a(:)
354  integer, intent(inout) :: x(:)
355  integer, intent(in) :: left, right
356  integer, intent(in) :: k
357  integer :: middle, pivot
358 
359  if (left < right) then
360  if ((right-left)/k >= 2) then
361  call sample_second_min_arg(a, x, left, right, k, pivot)
362  else
363  middle = (left + right) / 2
364  call median_of_three_arg(a, x, left, middle, right, pivot)
365  end if
366  call partition_arg(a, x, left, right, pivot)
367  if (k < pivot) then
368  call sort_quickselect_arg(a, x, left, pivot-1, k)
369  else if (k > pivot) then
370  call sort_quickselect_arg(a, x, pivot+1, right, k)
371  end if
372  end if
373 
374  return
375  end subroutine sort_quickselect_arg
376 
377  recursive subroutine sort_quickselect_desc(A, left, right, K)
378  implicit none
379  real(rp), intent(inout) :: a(:)
380  integer, intent(in) :: left, right
381  integer, intent(in) :: k
382  integer :: middle, pivot
383 
384  if (left < right) then
385  if ((right-left)/k >= 2) then
386  call sample_second_max(a, left, right, k, pivot)
387  else
388  middle = (left + right) / 2
389  call median_of_three(a, left, middle, right, pivot)
390  end if
391  call partition_desc(a, left, right, pivot)
392  if (k < pivot) then
393  call sort_quickselect_desc(a, left, pivot-1, k)
394  else if (k > pivot) then
395  call sort_quickselect_desc(a, pivot+1, right, k)
396  end if
397  end if
398 
399  return
400  end subroutine sort_quickselect_desc
401 
402  recursive subroutine sort_quickselect_desc_arg(A, X, left, right, K)
403  implicit none
404  real(rp), intent(in) :: a(:)
405  integer, intent(inout) :: x(:)
406  integer, intent(in) :: left, right
407  integer, intent(in) :: k
408  integer :: middle, pivot
409 
410  if (left < right) then
411  if ((right-left)/k >= 2) then
412  call sample_second_max_arg(a, x, left, right, k, pivot)
413  else
414  middle = (left + right) / 2
415  call median_of_three_arg(a, x, left, middle, right, pivot)
416  end if
417  call partition_desc_arg(a, x, left, right, pivot)
418  if (k < pivot) then
419  call sort_quickselect_desc_arg(a, x, left, pivot-1, k)
420  else if (k > pivot) then
421  call sort_quickselect_desc_arg(a, x, pivot+1, right, k)
422  end if
423  end if
424 
425  return
426  end subroutine sort_quickselect_desc_arg
427 
428  ! private
429 
430  subroutine swap(x, y)
431  implicit none
432  real(rp), intent(inout) :: x, y
433  real(rp) :: tmp
434  tmp = x
435  x = y
436  y = tmp
437 
438  return
439  end subroutine swap
440 
441  subroutine swap_i(x, y)
442  implicit none
443  integer, intent(inout) :: x, y
444  integer :: tmp
445  tmp = x
446  x = y
447  y = tmp
448 
449  return
450  end subroutine swap_i
451 
452  subroutine partition(A, left, right, pivot)
453  implicit none
454  real(rp), intent(inout) :: a(:)
455  integer, intent(in) :: left, right
456  integer, intent(inout) :: pivot
457  real(rp) :: a_pivot
458  integer :: idx, store_idx
459 
460  a_pivot = a(pivot)
461  call swap(a(pivot), a(right))
462 
463  store_idx = left
464  do idx = left, right-1
465  if (a(idx) < a_pivot) then
466  call swap(a(store_idx), a(idx))
467  store_idx = store_idx + 1
468  end if
469  end do
470 
471  call swap(a(right), a(store_idx))
472  pivot = store_idx
473 
474  return
475  end subroutine partition
476 
477  subroutine partition_arg(A, X, left, right, pivot)
478  implicit none
479  real(rp), intent(in) :: a(:)
480  integer, intent(inout) :: x(:)
481  integer, intent(in) :: left, right
482  integer, intent(inout) :: pivot
483  real(rp) :: a_pivot
484  integer :: idx, store_idx
485 
486  a_pivot = a(x(pivot))
487  call swap_i(x(pivot), x(right))
488 
489  store_idx = left
490  do idx = left, right-1
491  if (a(x(idx)) < a_pivot) then
492  call swap_i(x(store_idx), x(idx))
493  store_idx = store_idx + 1
494  end if
495  end do
496 
497  call swap_i(x(right), x(store_idx))
498  pivot = store_idx
499 
500  return
501  end subroutine partition_arg
502 
503  subroutine partition_desc(A, left, right, pivot)
504  implicit none
505  real(rp), intent(inout) :: a(:)
506  integer, intent(in) :: left, right
507  integer, intent(inout) :: pivot
508  real(rp) :: a_pivot
509  integer :: idx, store_idx
510 
511  a_pivot = a(pivot)
512  call swap(a(pivot), a(right))
513 
514  store_idx = left
515  do idx = left, right-1
516  if (a(idx) > a_pivot) then
517  call swap(a(store_idx), a(idx))
518  store_idx = store_idx + 1
519  end if
520  end do
521 
522  call swap(a(right), a(store_idx))
523  pivot = store_idx
524 
525  return
526  end subroutine partition_desc
527 
528  subroutine partition_desc_arg(A, X, left, right, pivot)
529  implicit none
530  real(rp), intent(in) :: a(:)
531  integer, intent(inout) :: x(:)
532  integer, intent(in) :: left, right
533  integer, intent(inout) :: pivot
534  real(rp) :: a_pivot
535  integer :: idx, store_idx
536 
537  a_pivot = a(x(pivot))
538  call swap_i(x(pivot), x(right))
539 
540  store_idx = left
541  do idx = left, right-1
542  if (a(x(idx)) > a_pivot) then
543  call swap_i(x(store_idx), x(idx))
544  store_idx = store_idx + 1
545  end if
546  end do
547 
548  call swap_i(x(right), x(store_idx))
549  pivot = store_idx
550 
551  return
552  end subroutine partition_desc_arg
553 
554  subroutine median_of_three(A, i1, i2, i3, i)
555  implicit none
556  real(rp), intent(in) :: a(:)
557  integer, intent(in) :: i1, i2, i3
558  integer, intent(out) :: i
559 
560  if (a(i1) < a(i2)) then
561  if (a(i2) < a(i3)) then
562  i = i2
563  else if (a(i1) < a(i3)) then
564  i = i3
565  else
566  i = i1
567  end if
568  else
569  if (a(i1) < a(i3)) then
570  i = i1
571  else if (a(i2) < a(i3)) then
572  i = i3
573  else
574  i = i2
575  end if
576  end if
577 
578  return
579  end subroutine median_of_three
580 
581  subroutine median_of_three_arg(A, X, i1, i2, i3, i)
582  implicit none
583  real(rp), intent(in) :: a(:)
584  integer, intent(in) :: x(:)
585  integer, intent(in) :: i1, i2, i3
586  integer, intent(out) :: i
587 
588  if (a(x(i1)) < a(x(i2))) then
589  if (a(x(i2)) < a(x(i3))) then
590  i = i2
591  else if (a(x(i1)) < a(x(i3))) then
592  i = i3
593  else
594  i = i1
595  end if
596  else
597  if (a(x(i1)) < a(x(i3))) then
598  i = i1
599  else if (a(x(i2)) < a(x(i3))) then
600  i = i3
601  else
602  i = i2
603  end if
604  end if
605 
606  return
607  end subroutine median_of_three_arg
608 
609  subroutine sample_second_min(A, left, right, K, i_2)
610  implicit none
611  real(rp), intent(in) :: a(:)
612  integer, intent(in) :: left, right
613  integer, intent(in) :: k
614  integer, intent(out) :: i_2
615  real(rp) :: a_min, a_min_2
616  integer :: i, j
617 
618  i = left
619  i_2 = left
620  a_min = huge(a)
621  a_min_2 = huge(a)
622  do j = left, right, k
623  if (a(j) < a_min) then
624  a_min_2 = a_min
625  a_min = a(j)
626  i_2 = i
627  i = j
628  else if (a(j) < a_min_2) then
629  a_min_2 = a(j)
630  i_2 = j
631  end if
632  end do
633 
634  return
635  end subroutine sample_second_min
636 
637  subroutine sample_second_min_arg(A, X, left, right, K, i_2)
638  implicit none
639  real(rp), intent(in) :: a(:)
640  integer, intent(in) :: x(:)
641  integer, intent(in) :: left, right
642  integer, intent(in) :: k
643  integer, intent(out) :: i_2
644  real(rp) :: a_min, a_min_2
645  integer :: i, j
646 
647  i = left
648  i_2 = left
649  a_min = huge(a)
650  a_min_2 = huge(a)
651  do j = left, right, k
652  if (a(x(j)) < a_min) then
653  a_min_2 = a_min
654  a_min = a(x(j))
655  i_2 = i
656  i = j
657  else if (a(x(j)) < a_min_2) then
658  a_min_2 = a(x(j))
659  i_2 = j
660  end if
661  end do
662 
663  return
664  end subroutine sample_second_min_arg
665 
666  subroutine sample_second_max(A, left, right, K, i_2)
667  implicit none
668  real(rp), intent(in) :: a(:)
669  integer, intent(in) :: left, right
670  integer, intent(in) :: k
671  integer, intent(out) :: i_2
672  real(rp) :: a_max, a_max_2
673  integer :: i, j
674 
675  i = left
676  i_2 = left
677  a_max = 0.0d0 - huge(a)
678  a_max_2 = 0.0d0 - huge(a)
679  do j = left, right, k
680  if (a(j) > a_max) then
681  a_max_2 = a_max
682  a_max = a(j)
683  i_2 = i
684  i = j
685  else if (a(j) > a_max_2) then
686  a_max_2 = a(j)
687  i_2 = j
688  end if
689  end do
690 
691  return
692  end subroutine sample_second_max
693 
694  subroutine sample_second_max_arg(A, X, left, right, K, i_2)
695  implicit none
696  real(rp), intent(in) :: a(:)
697  integer, intent(in) :: x(:)
698  integer, intent(in) :: left, right
699  integer, intent(in) :: k
700  integer, intent(out) :: i_2
701  real(rp) :: a_max, a_max_2
702  integer :: i, j
703 
704  i = left
705  i_2 = left
706  a_max = 0.0d0 - huge(a)
707  a_max_2 = 0.0d0 - huge(a)
708  do j = left, right, k
709  if (a(x(j)) > a_max) then
710  a_max_2 = a_max
711  a_max = a(x(j))
712  i_2 = i
713  i = j
714  else if (a(x(j)) > a_max_2) then
715  a_max_2 = a(x(j))
716  i_2 = j
717  end if
718  end do
719 
720  return
721  end subroutine sample_second_max_arg
722 
723 end module scale_sort
scale_sort::sort_quickselect_desc_arg
recursive subroutine, public sort_quickselect_desc_arg(A, X, left, right, K)
Definition: scale_sort.F90:403
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_sort
module SORT
Definition: scale_sort.F90:11
scale_sort::sort_quickselect_desc
recursive subroutine, public sort_quickselect_desc(A, left, right, K)
Definition: scale_sort.F90:378
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_io
module STDIO
Definition: scale_io.F90:10
scale_sort::sort_quickselect_arg
recursive subroutine, public sort_quickselect_arg(A, X, left, right, K)
Definition: scale_sort.F90:352
scale_sort::sort_quicksort
recursive subroutine, public sort_quicksort(n, array)
Definition: scale_sort.F90:219
scale_sort::sort_quickselect
recursive subroutine, public sort_quickselect(A, left, right, K)
Definition: scale_sort.F90:327
scale_sort::sort_heapsort
subroutine, public sort_heapsort(n, array)
Definition: scale_sort.F90:284
scale_sort::sort_uniq_int_sorted
subroutine, public sort_uniq_int_sorted(n, ary, c)
Definition: scale_sort.F90:200
Generated by   doxygen 1.8.17