1!
2!     SLATEC: public domain
3!
4*deck isort
5      subroutine isortiddc(ix,dy1,dy2,cy,n,kflag)
6!
7!     modified to make the same interchanges in two
8!     double (dy1 and dy2) and a char*20 aray (cy)
9!
10C***BEGIN PROLOGUE  ISORT
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  N6A2A
17C***TYPE      INTEGER (SSORT-S, DSORT-D, ISORT-I)
18C***KEYWORDS  SINGLETON QUICKSORT, SORT, SORTING
19C***AUTHOR  Jones, R. E., (SNLA)
20C           Kahaner, D. K., (NBS)
21C           Wisniewski, J. A., (SNLA)
22C***DESCRIPTION
23C
24C   ISORT sorts array IX and optionally makes the same interchanges in
25C   array IY.  The array IX may be sorted in increasing order or
26C   decreasing order.  A slightly modified quicksort algorithm is used.
27C
28C   Description of Parameters
29C      IX - integer array of values to be sorted
30C      IY - integer array to be (optionally) carried along
31C      N  - number of values in integer array IX to be sorted
32C      KFLAG - control parameter
33C            =  2  means sort IX in increasing order and carry IY along.
34C            =  1  means sort IX in increasing order (ignoring IY)
35C            = -1  means sort IX in decreasing order (ignoring IY)
36C            = -2  means sort IX 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***ROUTINES CALLED  XERMSG
42C***REVISION HISTORY  (YYMMDD)
43C   761118  DATE WRITTEN
44C   810801  Modified by David K. Kahaner.
45C   890531  Changed all specific intrinsics to generic.  (WRB)
46C   890831  Modified array declarations.  (WRB)
47C   891009  Removed unreferenced statement labels.  (WRB)
48C   891009  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 IX,IY. (M. McClain)
52C   920501  Reformatted the REFERENCES section.  (DWL, WRB)
53C   920519  Clarified error messages.  (DWL)
54C   920801  Declarations section rebuilt and code restructured to use
55C           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
56!   100411  changed the dimension of IL and IU from 21 to 31.
57!
58!     field IL and IU have the dimension 31. This is log2 of the largest
59!     array size to be sorted. If arrays larger than 2**31 in length have
60!     to be sorted, this dimension has to be modified accordingly
61!
62C***END PROLOGUE  ISORT
63!
64      implicit none
65C     .. Scalar Arguments ..
66      integer kflag, n,iside,istat
67C     .. Array Arguments ..
68      integer ix(2,*)
69      real*8 dy1(2,*),dy2(2,*)
70      character*20 cy(*)
71C     .. Local Scalars ..
72      real r
73      integer i, ij, j, k, kk, l, m, nn, t, tt
74      real*8 tty11,tty12,ty11,ty12,tty21,tty22,ty21,ty22,ttx2,tx2
75      character*20 uuy,uy
76C     .. Local Arrays ..
77      integer il(31), iu(31)
78C     .. External Subroutines ..
79!      EXTERNAL XERMSG
80C     .. Intrinsic Functions ..
81      intrinsic abs, int
82C***FIRST EXECUTABLE STATEMENT  ISORT
83!
84      do i=1,n
85         read(cy(i)(2:2),'(i1)',iostat=istat) iside
86         if(istat.gt.0) iside=0
87         ix(1,i)=10*ix(1,i)+iside
88      enddo
89!
90      nn = n
91      if (nn .lt. 1) then
92!         CALL XERMSG ('SLATEC', 'ISORT',
93!     +      'The number of values to be sorted is not positive.', 1, 1)
94         return
95      endif
96C
97      kk = abs(kflag)
98      if (kk.ne.1 .and. kk.ne.2) then
99!         CALL XERMSG ('SLATEC', 'ISORT',
100!     +      'The sort control parameter, K, is not 2, 1, -1, or -2.', 2,
101!     +      1)
102         return
103      endif
104C
105C     Alter array IX to get decreasing order if needed
106C
107      if (kflag .le. -1) then
108         do 10 i=1,nn
109            ix(1,i) = -ix(1,i)
110   10    continue
111      endif
112C
113      if (kk .eq. 2) go to 100
114C
115C     Sort IX only
116C
117      m = 1
118      i = 1
119      j = nn
120      r = 0.375e0
121C
122   20 if (i .eq. j) go to 60
123      if (r .le. 0.5898437e0) then
124         r = r+3.90625e-2
125      else
126         r = r-0.21875e0
127      endif
128C
129   30 k = i
130C
131C     Select a central element of the array and save it in location T
132C
133      ij = i + int((j-i)*r)
134      t = ix(1,ij)
135C
136C     If first element of array is greater than T, interchange with T
137C
138      if (ix(1,i) .gt. t) then
139         ix(1,ij) = ix(1,i)
140         ix(1,i) = t
141         t = ix(1,ij)
142      endif
143      l = j
144C
145C     If last element of array is less than than T, interchange with T
146C
147      if (ix(1,j) .lt. t) then
148         ix(1,ij) = ix(1,j)
149         ix(1,j) = t
150         t = ix(1,ij)
151C
152C        If first element of array is greater than T, interchange with T
153C
154         if (ix(1,i) .gt. t) then
155            ix(1,ij) = ix(1,i)
156            ix(1,i) = t
157            t = ix(1,ij)
158         endif
159      endif
160C
161C     Find an element in the second half of the array which is smaller
162C     than T
163C
164   40 l = l-1
165      if (ix(1,l) .gt. t) go to 40
166C
167C     Find an element in the first half of the array which is greater
168C     than T
169C
170   50 k = k+1
171      if (ix(1,k) .lt. t) go to 50
172C
173C     Interchange these elements
174C
175      if (k .le. l) then
176         tt = ix(1,l)
177         ix(1,l) = ix(1,k)
178         ix(1,k) = tt
179         go to 40
180      endif
181C
182C     Save upper and lower subscripts of the array yet to be sorted
183C
184      if (l-i .gt. j-k) then
185         il(m) = i
186         iu(m) = l
187         i = k
188         m = m+1
189      else
190         il(m) = k
191         iu(m) = j
192         j = l
193         m = m+1
194      endif
195      go to 70
196C
197C     Begin again on another portion of the unsorted array
198C
199   60 m = m-1
200      if (m .eq. 0) go to 190
201      i = il(m)
202      j = iu(m)
203C
204   70 if (j-i .ge. 1) go to 30
205      if (i .eq. 1) go to 20
206      i = i-1
207C
208   80 i = i+1
209      if (i .eq. j) go to 60
210      t = ix(1,i+1)
211      if (ix(1,i) .le. t) go to 80
212      k = i
213C
214   90 ix(1,k+1) = ix(1,k)
215      k = k-1
216      if (t .lt. ix(1,k)) go to 90
217      ix(1,k+1) = t
218      go to 80
219C
220C     Sort IX and carry IY along
221C
222  100 m = 1
223      i = 1
224      j = nn
225      r = 0.375e0
226C
227  110 if (i .eq. j) go to 150
228      if (r .le. 0.5898437e0) then
229         r = r+3.90625e-2
230      else
231         r = r-0.21875e0
232      endif
233C
234  120 k = i
235C
236C     Select a central element of the array and save it in location T
237C
238      ij = i + int((j-i)*r)
239      t = ix(1,ij)
240      ty11 = dy1(1,ij)
241      ty21 = dy1(2,ij)
242      ty12 = dy2(1,ij)
243      ty22 = dy2(2,ij)
244      tx2 = ix(2,ij)
245      uy = cy(ij)
246C
247C     If first element of array is greater than T, interchange with T
248C
249      if (ix(1,i) .gt. t) then
250         ix(1,ij) = ix(1,i)
251         ix(1,i) = t
252         t = ix(1,ij)
253         dy1(1,ij) = dy1(1,i)
254         dy1(2,ij) = dy1(2,i)
255         dy2(1,ij) = dy2(1,i)
256         dy2(2,ij) = dy2(2,i)
257         ix(2,ij) = ix(2,i)
258         cy(ij) = cy(i)
259         dy1(1,i) = ty11
260         dy1(2,i) = ty21
261         dy2(1,i) = ty12
262         dy2(2,i) = ty22
263         ix(2,i) = tx2
264         cy(i) = uy
265         ty11 = dy1(1,ij)
266         ty21 = dy1(2,ij)
267         ty12 = dy2(1,ij)
268         ty22 = dy2(2,ij)
269         tx2 = ix(2,ij)
270         uy = cy(ij)
271      endif
272      l = j
273C
274C     If last element of array is less than T, interchange with T
275C
276      if (ix(1,j) .lt. t) then
277         ix(1,ij) = ix(1,j)
278         ix(1,j) = t
279         t = ix(1,ij)
280         dy1(1,ij) = dy1(1,j)
281         dy1(2,ij) = dy1(2,j)
282         dy2(1,ij) = dy2(1,j)
283         dy2(2,ij) = dy2(2,j)
284         ix(2,ij) = ix(2,j)
285         cy(ij) = cy(j)
286         dy1(1,j) = ty11
287         dy1(2,j) = ty21
288         dy2(1,j) = ty12
289         dy2(2,j) = ty22
290         ix(2,j) = tx2
291         cy(j) = uy
292         ty11 = dy1(1,ij)
293         ty21 = dy1(2,ij)
294         ty12 = dy2(1,ij)
295         ty22 = dy2(2,ij)
296         tx2 = ix(2,ij)
297         uy = cy(ij)
298C
299C        If first element of array is greater than T, interchange with T
300C
301         if (ix(1,i) .gt. t) then
302            ix(1,ij) = ix(1,i)
303            ix(1,i) = t
304            t = ix(1,ij)
305            dy1(1,ij) = dy1(1,i)
306            dy1(2,ij) = dy1(2,i)
307            dy2(1,ij) = dy2(1,i)
308            dy2(2,ij) = dy2(2,i)
309            ix(2,ij) = ix(2,i)
310            cy(ij) = cy(i)
311            dy1(1,i) = ty11
312            dy1(2,i) = ty21
313            dy2(1,i) = ty12
314            dy2(2,i) = ty22
315            ix(2,i) = tx2
316            cy(i) = uy
317            ty11 = dy1(1,ij)
318            ty21 = dy1(2,ij)
319            ty12 = dy2(1,ij)
320            ty22 = dy2(2,ij)
321            tx2 = ix(2,ij)
322            uy = cy(ij)
323         endif
324      endif
325C
326C     Find an element in the second half of the array which is smaller
327C     than T
328C
329  130 l = l-1
330      if (ix(1,l) .gt. t) go to 130
331C
332C     Find an element in the first half of the array which is greater
333C     than T
334C
335  140 k = k+1
336      if (ix(1,k) .lt. t) go to 140
337C
338C     Interchange these elements
339C
340      if (k .le. l) then
341         tt = ix(1,l)
342         ix(1,l) = ix(1,k)
343         ix(1,k) = tt
344         tty11 = dy1(1,l)
345         tty21 = dy1(2,l)
346         tty12 = dy2(1,l)
347         tty22 = dy2(2,l)
348         ttx2 = ix(2,l)
349         uuy = cy(l)
350         dy1(1,l) = dy1(1,k)
351         dy1(2,l) = dy1(2,k)
352         dy2(1,l) = dy2(1,k)
353         dy2(2,l) = dy2(2,k)
354         ix(2,l) = ix(2,k)
355         cy(l) = cy(k)
356         dy1(1,k) = tty11
357         dy1(2,k) = tty21
358         dy2(1,k) = tty12
359         dy2(2,k) = tty22
360         ix(2,k) = ttx2
361         cy(k) = uuy
362         go to 130
363      endif
364C
365C     Save upper and lower subscripts of the array yet to be sorted
366C
367      if (l-i .gt. j-k) then
368         il(m) = i
369         iu(m) = l
370         i = k
371         m = m+1
372      else
373         il(m) = k
374         iu(m) = j
375         j = l
376         m = m+1
377      endif
378      go to 160
379C
380C     Begin again on another portion of the unsorted array
381C
382  150 m = m-1
383      if (m .eq. 0) go to 190
384      i = il(m)
385      j = iu(m)
386C
387  160 if (j-i .ge. 1) go to 120
388      if (i .eq. 1) go to 110
389      i = i-1
390C
391  170 i = i+1
392      if (i .eq. j) go to 150
393      t = ix(1,i+1)
394      ty11 = dy1(1,i+1)
395      ty21 = dy1(2,i+1)
396      ty12 = dy2(1,i+1)
397      ty22 = dy2(2,i+1)
398      tx2 = ix(2,i+1)
399      uy = cy(i+1)
400      if (ix(1,i) .le. t) go to 170
401      k = i
402C
403  180 ix(1,k+1) = ix(1,k)
404      dy1(1,k+1) = dy1(1,k)
405      dy1(2,k+1) = dy1(2,k)
406      dy2(1,k+1) = dy2(1,k)
407      dy2(2,k+1) = dy2(2,k)
408      ix(2,k+1) = ix(2,k)
409      cy(k+1) = cy(k)
410      k = k-1
411      if (t .lt. ix(1,k)) go to 180
412      ix(1,k+1) = t
413      dy1(1,k+1) = ty11
414      dy1(2,k+1) = ty21
415      dy2(1,k+1) = ty12
416      dy2(2,k+1) = ty22
417      ix(2,k+1) = tx2
418      cy(k+1) = uy
419      go to 170
420C
421C     Clean up
422C
423  190 if (kflag .le. -1) then
424         do 200 i=1,nn
425            ix(1,i) = -ix(1,i)
426  200    continue
427      endif
428!
429      do i=1,nn
430         read(cy(i)(2:2),'(i1)',iostat=istat) iside
431         if(istat.gt.0) iside=0
432         ix(1,i)=(ix(1,i)-iside)/10
433      enddo
434!
435      return
436      end
437