1      SUBROUTINE RANGEMENT(NEQ)
2      INCLUDE 'com_faces.f'
3C      REAL*8 SEC0,SEC1
4C
5C      CALL TEMPS(SEC0,IS)
6C      CALL MYSORT_PLUS(CENTR,CENTRMI,CENTRMA,IORDRE,NEQ)
7      CALL MYSORT(CENTR,IORDRE,NEQ)
8C      CALL TEMPS(SEC1,IS)
9C      PRINT*,'Sort    =',SEC1-SEC0
10      END
11C-----------------------------------------------------------------------
12      SUBROUTINE MYSORT(ARRAY,IARRAY,NEQ)
13CCBEGIN
14C-----------------------------------------------------------------------
15C
16C
17C     PURPOSE
18C     ------
19C
20C     CLASSE LES ELEMENTS DU TABLEAU DE REEL "ARRAY"
21C     ET CONSERVE LA TRACE DE LA PERMUTATION DANS LE TABLEAU
22C     D'ENTIERS "IARRAY".
23C
24CCDOC
25C     INPUT ARGUMENTS
26C     --------------
27C
28C     NEQ     = LENGTH OF IARRAY.
29C
30C     INPUT-OUTPUT ARGUMENTS
31C     ----------------------
32C
33C     IARRAY  = INTEGER ARRAY.
34C
35C     REFERENCE
36C     ---------
37C
38C     LORIN H.,SORTING AND SORTING SYSTEMS,
39C     ADDISON-WESLEY PUBLISHING COMPANY,1975.A-42.
40C     ACM ALGORITHM 347.
41C
42C     NOTES
43C     -----
44C
45C     INTERNAL WORK AREA IWR(2,K) PERMITS SORTING UP TO
46C     2**(K+1)-1 ELEMENTS. HERE K IS SET TO 32.
47C
48C-----------------------------------------------------------------------
49C
50CCEND
51C
52C
53      DIMENSION IARRAY(NEQ),IWR(2,32)
54      REAL*4 ARRAY(NEQ)
55C
56C     INITIALISATION
57C
58      DO 200 I=1,NEQ
59        IARRAY(I)=I
60  200 CONTINUE
61      M = 1
62      I = 1
63      J = NEQ
64    5 IF(I.GE.J)  GO TO 70
65C
66C--------BEGIN SPLIT OF SEGMENT--------------------------------------
67C
68C        FIND MEDIAN OF FIRST, MIDDLE AND LAST ELEMENTS OF SEGMENT, IMID
69C
70   10    K = I
71         IJ2 = (J+I)/2
72         IMID = IARRAY(IJ2)
73         AMID = ARRAY(IJ2)
74         IF(ARRAY(I).LE.AMID)  GOTO 20
75            IARRAY(IJ2) = IARRAY(I)
76            IARRAY(I) = IMID
77            IMID = IARRAY(IJ2)
78            ARRAY(IJ2) = ARRAY(I)
79            ARRAY(I) = AMID
80            AMID = ARRAY(IJ2)
81   20    L = J
82         IF(ARRAY(J).GE.AMID)  GO TO 40
83            IARRAY(IJ2) = IARRAY(J)
84            IARRAY(J) = IMID
85            IMID = IARRAY(IJ2)
86            ARRAY(IJ2) = ARRAY(J)
87            ARRAY(J) = AMID
88            AMID = ARRAY(IJ2)
89            IF(ARRAY(I).LE.AMID)  GO TO 40
90               IARRAY(IJ2) = IARRAY(I)
91               IARRAY(I) = IMID
92               IMID = IARRAY(IJ2)
93               ARRAY(IJ2) = ARRAY(I)
94               ARRAY(I) = AMID
95               AMID = ARRAY(IJ2)
96            GO TO 40
97C
98C        SPLIT SEGMENT AROUND IMID, SORTING ELEMENTS INTO THOSE
99C        GT. AND LT. IMID.
100C
101   30    IARRAY(L) = IARRAY(K)
102         IARRAY(K) = ITT
103         ARRAY(L) = ARRAY(K)
104         ARRAY(K) = ATT
105   40    L = L - 1
106         IF(ARRAY(L).GT.AMID)  GO TO 40
107         ITT = IARRAY(L)
108         ATT = ARRAY(L)
109   50    K = K + 1
110         IF(ARRAY(K).LT.AMID)  GO TO 50
111         IF(K.LE.L)  GO TO 30
112C
113C--------END OF SPLIT OF SEGMENT---------------------------------------
114C
115C     STORE LIMITS OF LARGER SEGMENT IN IWR
116C
117      IF(L-I.LE.J-K)  GO TO 60
118         IWR(2,M) = I
119         IWR(1,M) = L
120         I = K
121         M = M + 1
122         GO TO 80
123   60 IWR(2,M) = K
124      IWR(1,M) = J
125      J = L
126      M = M + 1
127      GO TO 80
128C
129C     WHEN SEGMENT HAS BEEN SORTED, GET LIMITS OF NEXT SEGMENT.
130C     WHEN ALL SEGMENTS HAVE BEEN SORTED, SORT IS COMPLETE.
131C
132   70 M = M - 1
133      IF(M.EQ.0)  RETURN
134      I = IWR(2,M)
135      J = IWR(1,M)
136C
137C     SORT SHORT SEQUENCES (LT.11 ELEMENTS) BY INTERCHANGE OF ADJACENT
138C     PAIRS (EXCEPT INITIAL SEGMENT)
139C
140   80 IF(J-I.GE.11)  GO TO 10
141      IF(I.EQ.1)  GO TO 5
142      I = I - 1
143   90 I = I + 1
144      IF(I.EQ.J)  GO TO 70
145      IMID = IARRAY(I+1)
146      AMID = ARRAY(I+1)
147      IF(ARRAY(I).LE.AMID)  GO TO 90
148      K = I
149  100 IARRAY(K+1) = IARRAY(K)
150      ARRAY(K+1) = ARRAY(K)
151      K = K - 1
152      IF(AMID.LT.ARRAY(K))  GO TO 100
153      IARRAY(K+1) = IMID
154      ARRAY(K+1) = AMID
155      GO TO 90
156C
157      END
158C-----------------------------------------------------------------------
159      SUBROUTINE MYSORT_PLUS(ARRAY,ARRAYMI,ARRAYMA,IARRAY,NEQ)
160CCBEGIN
161C-----------------------------------------------------------------------
162C
163C
164C     PURPOSE
165C     ------
166C
167C     CLASSE LES ELEMENTS DU TABLEAU DE REEL "ARRAY"
168C     ET CONSERVE LA TRACE DE LA PERMUTATION DANS LE TABLEAU
169C     D'ENTIERS "IARRAY".
170C
171CCDOC
172C     INPUT ARGUMENTS
173C     --------------
174C
175C     NEQ     = LENGTH OF IARRAY.
176C
177C     INPUT-OUTPUT ARGUMENTS
178C     ----------------------
179C
180C     IARRAY  = INTEGER ARRAY.
181C
182C     REFERENCE
183C     ---------
184C
185C     LORIN H.,SORTING AND SORTING SYSTEMS,
186C     ADDISON-WESLEY PUBLISHING COMPANY,1975.A-42.
187C     ACM ALGORITHM 347.
188C
189C     NOTES
190C     -----
191C
192C     INTERNAL WORK AREA IWR(2,K) PERMITS SORTING UP TO
193C     2**(K+1)-1 ELEMENTS. HERE K IS SET TO 32.
194C
195C-----------------------------------------------------------------------
196C
197CCEND
198C
199C
200      INTEGER IARRAY(*),IWR(2,32)
201      DIMENSION ARRAY(*),ARRAYMI(*),ARRAYMA(*)
202      logical pluspetit,plusgrand
203      include 'com_sort.f'
204C
205C INITIALISATION
206C
207      nboper = 0
208      nbfacile = 0
209      DO I=1,NEQ
210        IARRAY(I)=I
211      ENDDO
212      M = 1
213      I = 1
214      J = NEQ
215      do while (i.eq.1)
216        IF (I.lt.J) then
217C
218C--------BEGIN SPLIT OF SEGMENT--------------------------------------
219C
220C FIND MEDIAN OF FIRST, MIDDLE AND LAST ELEMENTS OF SEGMENT, IMID
221C
222 10       K = I
223          IJ2 = (J+I)/2
224          call affectecour(array,arraymi,arrayma,iarray,ij2)
225          if (plusgrand(array,arraymi,arrayma,i))
226     &         call echange(array,arraymi,arrayma,iarray,i,ij2)
227          L = J
228          IF (pluspetit(array,arraymi,arrayma,j)) then
229            call echange(array,arraymi,arrayma,iarray,j,ij2)
230            if (plusgrand(array,arraymi,arrayma,i)) then
231              call echange(array,arraymi,arrayma,iarray,i,ij2)
232              GOTO 40
233C
234C SPLIT SEGMENT AROUND IMID, SORTING ELEMENTS INTO THOSE
235C GT. AND LT. IMID.
236C
237Cfj              call affecte(array,arraymi,arrayma,iarray,l,k)
238Cfj              call affectetmpinv(array,arraymi,arrayma,iarray,k)
239            endif
240          endif
241 40       L = L-1
242          do while (plusgrand(array,arraymi,arrayma,l))
243            L = L-1
244          enddo
245          call affectetmp(array,arraymi,arrayma,iarray,l)
246          K = K+1
247          do while (pluspetit(array,arraymi,arrayma,k))
248            K = K+1
249          enddo
250          IF (K.LE.L) then
251            call affecte(array,arraymi,arrayma,iarray,l,k)
252            call affectetmpinv(array,arraymi,arrayma,iarray,k)
253            goto 40
254          endif
255C
256C--------END OF SPLIT OF SEGMENT---------------------------------------
257C
258C     STORE LIMITS OF LARGER SEGMENT IN IWR
259C
260          IF (L-I.gt.J-K) then
261            IWR(2,M) = I
262            IWR(1,M) = L
263            I = K
264            M = M + 1
265            IF (J-I.GE.11) GOTO 10
266          endif
267          IWR(2,M) = K
268          IWR(1,M) = J
269          J = L
270          M = M + 1
271          IF (J-I.GE.11) GOTO 10
272        endif
273C
274C WHEN SEGMENT HAS BEEN SORTED, GET LIMITS OF NEXT SEGMENT.
275C WHEN ALL SEGMENTS HAVE BEEN SORTED, SORT IS COMPLETE.
276C
277 70     M = M-1
278        IF (M.EQ.0) then
279          print*,'Operations',nboper,nbfacile
280     &         ,100.*real(nbfacile)/real(nboper)
281          RETURN
282        endif
283        I = IWR(2,M)
284        J = IWR(1,M)
285C
286C SORT SHORT SEQUENCES (LT.11 ELEMENTS) BY INTERCHANGE OF ADJACENT
287C PAIRS (EXCEPT INITIAL SEGMENT)
288C
289        IF (J-I.GE.11) GOTO 10
290      enddo
291c
292      I = I-1
293 90   I = I+1
294      IF (I.EQ.J) GOTO 70
295      call affectecour(array,arraymi,arrayma,iarray,i+1)
296      if ( .not.(pluspetit(array,arraymi,arrayma,i)) ) then
297        k = i
298        call affecte(array,arraymi,arrayma,iarray,k+1,k)
299        k = k-1
300        do while ( plusgrand(array,arraymi,arrayma,k) )
301          call affecte(array,arraymi,arrayma,iarray,k+1,k)
302          k = k-1
303        enddo
304        call affectecourinv(array,arraymi,arrayma,iarray,k+1)
305      endif
306      GOTO 90
307C
308      END
309C-----------------------------------------------------------------------
310      logical function pluspetit(array,arraymi,arrayma,i)
311      dimension array(*),arraymi(*),arrayma(*)
312      include 'com_sort.f'
313c
314      pluspetit = (array(i).lt.acour)
315      nboper = nboper+1
316      if (array(i).lt.acour.and.arrayma(i).lt.acourmi)
317     &     nbfacile=nbfacile+1
318      end
319C-----------------------------------------------------------------------
320      logical function plusgrand(array,arraymi,arrayma,i)
321      dimension array(*),arraymi(*),arrayma(*)
322      include 'com_sort.f'
323c
324      plusgrand = (array(i).gt.acour)
325      nboper = nboper+1
326      if (array(i).gt.acour.and.arraymi(i).gt.acourma)
327     &     nbfacile=nbfacile+1
328      end
329C-----------------------------------------------------------------------
330      subroutine echange(array,arraymi,arrayma,iarray,i,j)
331      dimension array(*),arraymi(*),arrayma(*)
332      integer   iarray(*)
333      include 'com_sort.f'
334c
335      iarray(j) = iarray(i)
336      iarray(i) = icour
337      icour     = iarray(j)
338      array(j)  = array(i)
339      array(i)  = acour
340      acour     = array(j)
341      arraymi(j) = arraymi(i)
342      arraymi(i) = acourmi
343      acourmi    = arraymi(j)
344      arrayma(j) = arrayma(i)
345      arrayma(i) = acourma
346      acourma    = arrayma(j)
347      end
348C-----------------------------------------------------------------------
349      subroutine affecte(array,arraymi,arrayma,iarray,i,j)
350      dimension array(*),arraymi(*),arrayma(*)
351      integer   iarray(*)
352      include 'com_sort.f'
353c
354      iarray(i) = iarray(j)
355      array(i)  = array(j)
356      arraymi(i)  = arraymi(j)
357      arrayma(i)  = arrayma(j)
358      end
359C-----------------------------------------------------------------------
360      subroutine affectecour(array,arraymi,arrayma,iarray,i)
361      dimension array(*),arraymi(*),arrayma(*)
362      integer   iarray(*)
363      include 'com_sort.f'
364c
365      icour = iarray(i)
366      acour = array(i)
367      acourmi = arraymi(i)
368      acourma = arrayma(i)
369      end
370C-----------------------------------------------------------------------
371      subroutine affectecourinv(array,arraymi,arrayma,iarray,i)
372      dimension array(*),arraymi(*),arrayma(*)
373      integer   iarray(*)
374      include 'com_sort.f'
375c
376      iarray(i) = icour
377      array(i)  = acour
378      arraymi(i)  = acourmi
379      arrayma(i)  = acourma
380      end
381C-----------------------------------------------------------------------
382      subroutine affectetmp(array,arraymi,arrayma,iarray,i)
383      dimension array(*),arraymi(*),arrayma(*)
384      integer   iarray(*)
385      include 'com_sort.f'
386c
387      itmp = iarray(i)
388      atmp = array(i)
389      atmpmi = arraymi(i)
390      atmpma = arrayma(i)
391      end
392C-----------------------------------------------------------------------
393      subroutine affectetmpinv(array,arraymi,arrayma,iarray,i)
394      dimension array(*),arraymi(*),arrayma(*)
395      integer   iarray(*)
396      include 'com_sort.f'
397c
398      iarray(i) = itmp
399      array(i)  = atmp
400      arraymi(i)  = atmpmi
401      arrayma(i)  = atmpma
402      end
403C-----------------------------------------------------------------------
404      SUBROUTINE CORRIGE(IORDRE,NBON,NSURF,NEIS,NSENS)
405      INTEGER IORDRE(*),NEIS(*),NSENS(*)
406ccc      real*8 S1,S2
407C
408      N = NBON+NSURF
409      NPROF = 1000
410ccc      CALL TEMPS(S1,IS)
411      DO I=2,N
412        IF (IORDRE(I).LE.NBON) THEN
413          IF (NSENS(IORDRE(I)).LT.0) THEN
414            NFIN = MAX(1,I-NPROF)
415            JSURF = 0
416            DO J=I-1,NFIN,-1
417              IF (IORDRE(J).LE.NBON) THEN
418                GOTO 1
419              ELSE
420                IF (NEIS(IORDRE(J)-NBON).EQ.-NSENS(IORDRE(I))) JSURF=J
421              ENDIF
422            ENDDO
423 1          JJ = JSURF
424            IF (JJ.NE.0) THEN
425cc                print*,IORDRE(I),' recule de',i,' a',jj,nsens(IORDRE(I))
426              IBID = IORDRE(I)
427              DO J=I,JJ+1,-1
428                IORDRE(J) = IORDRE(J-1)
429              ENDDO
430              IORDRE(JJ) = IBID
431            ENDIF
432          ENDIF
433        ENDIF
434      ENDDO
435      DO I=N-1,1,-1
436        IF (IORDRE(I).LE.NBON) THEN
437          IF (NSENS(IORDRE(I)).GT.0) THEN
438            NFIN = MIN(N,I+NPROF)
439            JSURF = 0
440            DO J=I+1,NFIN
441              IF (IORDRE(J).LE.NBON) THEN
442                GOTO 2
443              ELSE
444                IF (NEIS(IORDRE(J)-NBON).EQ.NSENS(IORDRE(I))) JSURF=J
445              ENDIF
446            ENDDO
447 2          JJ = JSURF
448            IF (JJ.NE.0) THEN
449cc                print*,IORDRE(I),' avance de',i,' a',jj,nsens(IORDRE(I))
450              IBID = IORDRE(I)
451              DO J=I,JJ-1
452                IORDRE(J) = IORDRE(J+1)
453              ENDDO
454              IORDRE(JJ) = IBID
455            ENDIF
456          ENDIF
457        ENDIF
458      ENDDO
459ccc      CALL TEMPS(S2,IS)
460ccc      print*,real(s2-s1)
461      END
462