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