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