1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
!
! QUICKSORT - recursive version of the QuickSort algorithm
!
! by Wirth,N: Algorithm + Data Structure = Programs, Prentice-Hall, 1975
!
! Copyright © 1997-2011, 2017-8 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with Munipack. If not, see <http://www.gnu.org/licenses/>.
module quicksort
! precision of real numbers
integer, parameter, private :: rp = selected_real_kind(15)
private :: trid, tridi
contains
subroutine qsort(a,idx)
real(rp), dimension(:), intent(in out) :: a
integer, dimension(:), optional, intent (in out) :: idx
integer :: n
n = size(a)
if( n < 2 ) return
if( present(idx) ) then
if( n /= size(idx) ) stop 'qsort: n /= size(idx)'
call tridi(1,n,a,idx)
else
call trid(1,n,a)
end if
end subroutine qsort
recursive subroutine trid(l,r,a)
integer,intent(in) :: l,r
real(rp), dimension(:),intent(in out) :: a
integer :: i,j
real(rp) :: x,w ! internal buffers
i = l
j = r
x = a((l+r)/2)
do
do while ( a(i) < x )
i = i + 1
end do
do while ( x < a(j) )
j = j - 1
end do
if( i <= j )then
w = a(i); a(i) = a(j); a(j) = w
i = i + 1; j = j - 1
endif
if( i > j ) exit
end do
if( l < j ) call trid(l,j,a)
if( i < r ) call trid(i,r,a)
end subroutine trid
!--------------------------------------------------------------------------
recursive subroutine tridi(l,r,a,idx)
integer,intent(in) :: l, r
integer, dimension(:), intent(in out) :: idx
real(rp), dimension(:), intent(in out) :: a
integer :: i,j,m
real(rp) :: x,w ! internal buffers
i = l
j = r
x = a((l+r)/2)
do
do while ( a(i) < x )
i = i + 1
end do
do while ( x < a(j) )
j = j - 1
end do
if( i <= j )then
w = a(i); a(i) = a(j); a(j) = w
m = idx(i); idx(i) = idx(j); idx(j) = m
i = i + 1; j = j - 1
endif
if( i > j ) exit
end do
if( l < j ) call tridi(l,j,a,idx)
if( i < r ) call tridi(i,r,a,idx)
end subroutine tridi
end module quicksort
|
