1!
2!     SLATEC: public domain
3!
4*deck dsort
5      subroutine dsort (dx, iy, n, kflag)
6c
7c    slight change: XERMSG was removed; error messages are
8c                   led to the screen;
9c
10C***BEGIN PROLOGUE  DSORT
11C***PURPOSE  Sort an array and optionally make the same interchanges in
12C            an auxiliary array.  The array may be sorted in increasing
13C            or decreasing order.  A slightly modified QUICKSORT
14C            algorithm is used.
15C***LIBRARY   SLATEC
16C***CATEGORY  N6A2B
17C***TYPE      DOUBLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
18C***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
19C***AUTHOR  Jones, R. E., (SNLA)
20C           Wisniewski, J. A., (SNLA)
21C***ROUTINES CALLED  XERMSG
22C***DESCRIPTION
23C
24C   DSORT sorts array DX and optionally makes the same interchanges in
25C   array IY.  The array DX may be sorted in increasing order or
26C   decreasing order.  A slightly modified quicksort algorithm is used.
27C
28C   Description of Parameters
29C      DX - array of values to be sorted   (usually abscissas)
30C      IY - array to be (optionally) carried along
31C      N  - number of values in array DX to be sorted
32C      KFLAG - control parameter
33C            =  2  means sort DX in increasing order and carry IY along.
34C            =  1  means sort DX in increasing order (ignoring IY)
35C            = -1  means sort DX in decreasing order (ignoring IY)
36C            = -2  means sort DX in decreasing order and carry IY along.
37C
38C***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
39C                 for sorting with minimal storage, Communications of
40C                 the ACM, 12, 3 (1969), pp. 185-187.
41C***REVISION HISTORY  (YYMMDD)
42C   761101  DATE WRITTEN
43C   761118  Modified to use the Singleton quicksort algorithm.  (JAW)
44C   890531  Changed all specific intrinsics to generic.  (WRB)
45C   890831  Modified array declarations.  (WRB)
46C   891009  Removed unreferenced statement labels.  (WRB)
47C   891024  Changed category.  (WRB)
48C   891024  REVISION DATE from Version 3.2
49C   891214  Prologue converted to Version 4.0 format.  (BAB)
50C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
51C   901012  Declared all variables; changed X,Y to DX,IY; changed
52C           code to parallel SSORT. (M. McClain)
53C   920501  Reformatted the REFERENCES section.  (DWL, WRB)
54C   920519  Clarified error messages.  (DWL)
55C   920801  Declarations section rebuilt and code restructured to use
56C           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
57!   100411  changed the dimension of IL and IU from 21 to 31.
58!   150514  inserted intent statements
59!
60!     field IL and IU have the dimension 31. This is log2 of the largest
61!     array size to be sorted. If arrays larger than 2**31 in length have
62!     to be sorted, this dimension has to be modified accordingly
63!
64C***END PROLOGUE  DSORT
65      implicit none
66C     .. Scalar Arguments ..
67      integer kflag, n,iy(*),ty,tty
68C     .. Array Arguments ..
69      double precision dx(*)
70C     .. Local Scalars ..
71      double precision r, t, tt
72      integer i, ij, j, k, kk, l, m, nn
73C     .. Local Arrays ..
74      integer il(31), iu(31)
75C     .. External Subroutines ..
76c      EXTERNAL XERMSG
77C     .. Intrinsic Functions ..
78      intrinsic abs, int
79!
80!
81C***FIRST EXECUTABLE STATEMENT  DSORT
82      nn = n
83      if (nn .lt. 1) then
84         write(*,*) '*error in dsort: the number of values to be'
85         write(*,*) '       sorted is not positive: ',nn
86         call exit(201)
87      endif
88C
89      kk = abs(kflag)
90      if (kk.ne.1 .and. kk.ne.2) then
91         write(*,*) '*error in dsort: the sort control parameter is'
92         write(*,*) '       not 2, 1, -1 or -2'
93         call exit(201)
94      endif
95C
96C     Alter array DX to get decreasing order if needed
97C
98      if (kflag .le. -1) then
99         do 10 i=1,nn
100            dx(i) = -dx(i)
101   10    continue
102      endif
103C
104      if (kk .eq. 2) go to 100
105C
106C     Sort DX only
107C
108      m = 1
109      i = 1
110      j = nn
111      r = 0.375d0
112C
113   20 if (i .eq. j) go to 60
114      if (r .le. 0.5898437d0) then
115         r = r+3.90625d-2
116      else
117         r = r-0.21875d0
118      endif
119C
120   30 k = i
121C
122C     Select a central element of the array and save it in location T
123C
124      ij = i + int((j-i)*r)
125      t = dx(ij)
126C
127C     If first element of array is greater than T, interchange with T
128C
129      if (dx(i) .gt. t) then
130         dx(ij) = dx(i)
131         dx(i) = t
132         t = dx(ij)
133      endif
134      l = j
135C
136C     If last element of array is less than than T, interchange with T
137C
138      if (dx(j) .lt. t) then
139         dx(ij) = dx(j)
140         dx(j) = t
141         t = dx(ij)
142C
143C        If first element of array is greater than T, interchange with T
144C
145         if (dx(i) .gt. t) then
146            dx(ij) = dx(i)
147            dx(i) = t
148            t = dx(ij)
149         endif
150      endif
151C
152C     Find an element in the second half of the array which is smaller
153C     than T
154C
155   40 l = l-1
156      if (dx(l) .gt. t) go to 40
157C
158C     Find an element in the first half of the array which is greater
159C     than T
160C
161   50 k = k+1
162      if (dx(k) .lt. t) go to 50
163C
164C     Interchange these elements
165C
166      if (k .le. l) then
167         tt = dx(l)
168         dx(l) = dx(k)
169         dx(k) = tt
170         go to 40
171      endif
172C
173C     Save upper and lower subscripts of the array yet to be sorted
174C
175      if (l-i .gt. j-k) then
176         il(m) = i
177         iu(m) = l
178         i = k
179         m = m+1
180      else
181         il(m) = k
182         iu(m) = j
183         j = l
184         m = m+1
185      endif
186      go to 70
187C
188C     Begin again on another portion of the unsorted array
189C
190   60 m = m-1
191      if (m .eq. 0) go to 190
192      i = il(m)
193      j = iu(m)
194C
195   70 if (j-i .ge. 1) go to 30
196      if (i .eq. 1) go to 20
197      i = i-1
198C
199   80 i = i+1
200      if (i .eq. j) go to 60
201      t = dx(i+1)
202      if (dx(i) .le. t) go to 80
203      k = i
204C
205   90 dx(k+1) = dx(k)
206      k = k-1
207      if (t .lt. dx(k)) go to 90
208      dx(k+1) = t
209      go to 80
210C
211C     Sort DX and carry IY along
212C
213  100 m = 1
214      i = 1
215      j = nn
216      r = 0.375d0
217C
218  110 if (i .eq. j) go to 150
219      if (r .le. 0.5898437d0) then
220         r = r+3.90625d-2
221      else
222         r = r-0.21875d0
223      endif
224C
225  120 k = i
226C
227C     Select a central element of the array and save it in location T
228C
229      ij = i + int((j-i)*r)
230      t = dx(ij)
231      ty = iy(ij)
232C
233C     If first element of array is greater than T, interchange with T
234C
235      if (dx(i) .gt. t) then
236         dx(ij) = dx(i)
237         dx(i) = t
238         t = dx(ij)
239         iy(ij) = iy(i)
240         iy(i) = ty
241         ty = iy(ij)
242      endif
243      l = j
244C
245C     If last element of array is less than T, interchange with T
246C
247      if (dx(j) .lt. t) then
248         dx(ij) = dx(j)
249         dx(j) = t
250         t = dx(ij)
251         iy(ij) = iy(j)
252         iy(j) = ty
253         ty = iy(ij)
254C
255C        If first element of array is greater than T, interchange with T
256C
257         if (dx(i) .gt. t) then
258            dx(ij) = dx(i)
259            dx(i) = t
260            t = dx(ij)
261            iy(ij) = iy(i)
262            iy(i) = ty
263            ty = iy(ij)
264         endif
265      endif
266C
267C     Find an element in the second half of the array which is smaller
268C     than T
269C
270  130 l = l-1
271      if (dx(l) .gt. t) go to 130
272C
273C     Find an element in the first half of the array which is greater
274C     than T
275C
276  140 k = k+1
277      if (dx(k) .lt. t) go to 140
278C
279C     Interchange these elements
280C
281      if (k .le. l) then
282         tt = dx(l)
283         dx(l) = dx(k)
284         dx(k) = tt
285         tty = iy(l)
286         iy(l) = iy(k)
287         iy(k) = tty
288         go to 130
289      endif
290C
291C     Save upper and lower subscripts of the array yet to be sorted
292C
293      if (l-i .gt. j-k) then
294         il(m) = i
295         iu(m) = l
296         i = k
297         m = m+1
298      else
299         il(m) = k
300         iu(m) = j
301         j = l
302         m = m+1
303      endif
304      go to 160
305C
306C     Begin again on another portion of the unsorted array
307C
308  150 m = m-1
309      if (m .eq. 0) go to 190
310      i = il(m)
311      j = iu(m)
312C
313  160 if (j-i .ge. 1) go to 120
314      if (i .eq. 1) go to 110
315      i = i-1
316C
317  170 i = i+1
318      if (i .eq. j) go to 150
319      t = dx(i+1)
320      ty = iy(i+1)
321      if (dx(i) .le. t) go to 170
322      k = i
323C
324  180 dx(k+1) = dx(k)
325      iy(k+1) = iy(k)
326      k = k-1
327      if (t .lt. dx(k)) go to 180
328      dx(k+1) = t
329      iy(k+1) = ty
330      go to 170
331C
332C     Clean up
333C
334  190 if (kflag .le. -1) then
335         do 200 i=1,nn
336            dx(i) = -dx(i)
337  200    continue
338      endif
339      return
340      end
341