C   SUBROUTINE DSOPE(A,IR,N,NP,LL,P)
C   
C   THIS SUBROUTINE MAY BE USED FOR SORTING INTERNALLY AN ARRAY
C   A OF N REAL NUMBERS TO AN INCREASING ORDER. THE SUBROUTINE
C   ALSO DETERMINES A PERMUTATION ARRAY IR, WHICH GIVES FOR
C   THE FINAL ORDERING THE INDEXES OF THE ORIGINAL ORDERING.
C   CALLING OF THE PROGRAM:
C
C      CALL DSOPE(A,IR,N,NP,LL,P)
C
C   WHERE
C
C      A=  THE ARRAY TO BE SORTED (INPUT/OUTPUT); LENGTH N+1
C      IR= THE PERMUTATION ARRAY (INPUT/OUTPUT). THE USER MUST
C          INITIALIZE IR(I)=I FOR I= 1 TO N. AT THE
C          TERMINATION OF THE SUBROUTINE IR CONTAINS FOR THE
C          NEW ORDERING THE ORIGINAL INDEXES; K=IR(I) MEANS
C          THAT A(I) OF THE FINAL ORDERING WAS THE ELEMENT K
C          IN THE ORIGINAL ORDERING
C      NP= A SWITCH (INPUT)
C      LL= AN AUXILIARY ARRAY (INPUT), MUST BE DEFINED IN THE
C          CALLING PROGRAM TO BE AT LEAST OF LENGTH NP/5
C      P=  AN AUXILIARY ARRAY (INPUT), MUST BE DEFINED IN THE
C          CALLING PROGRAM TO BE AT LEAST OF THE LENGTH NP.
C
C      THE METHOD
C   THE SUBROUTINE HAS TWO MODES OF THE OPERATION,
C   RULED BY THE VALUES OF THE PARAMETERS N AND NP:
C    MODE 1: IF N IS LESS THAN 500 OR NP IS LESS THAN N
C   QUICK SORT IS USED.
C    MODE 2: OTHERWISE, THE METHOD IS DISTRIBUTIVE SORTING.
C   THE USER CAN SELECT THE MODE 1 BY CALLING THE PROGRAM WITH
C   NP=0. THEN THE CALLING PROGRAM MUST CONTAIN A DEFINITION
C   OF THE ARRAYS LL(1) AND P(1).
C   THE USER CAN SELECT THE MODE 2 TO BE USED BY CALLING THE
C   PROGRAM WITH THE PARAMETER VALUE NP GREATER THAN OR EQUAL
C   TO N. IN THAT CASE THE CALLING PROGRAM MUST CONTAIN A
C   DEFINITION OF THE ARRAYS LL(NP/5) AN P(NP). A(N+1) MUST
C   CONTAIN A VALUE WHICH IS GREATER THAN ANY OF VALUES OF
C   A(1),..A(N).

C
C   FOR FURTHER INFORMATION OF THE METHOD, SEE:
C
C      JARMO ERNVALL AND OLLI NEVALAINEN:
C      "PERFORMANCE TESTS WITH DISTRIBUTIVE SORTING PROGRAMS",
C      DEPT. OF COMP. SCI., UNIV. OF TURKU, FINLAND, 1981
C
C   ADDRESS OF THE AUTHORS:
C      OLLI NEVALAINEN
C      DEPERTMENT OF COMPUTER SCIENCE
C      UNIVERSITY OF TURKU
C      SF-20500 TURKU
C      FINLAND
C
C   DATE: 30.11.1981
C
C
 	SUBROUTINE DSOPE(A,IR,N,NP,LL,P)
	DIMENSION A(1),LL(1),P(1),IR(1)
	DIMENSION ISA(50),ISB(50)
	INTEGER R
	DATA M/10/,MR/14/
C   CHOOSE OF THE SORTING TECHNIQUE
	IF(N.LE.1) RETURN
	IF(N.LT.500) GO TO 88 
	IF(NP.GE.N) GO TO 100
C   MODE 1: ONLY ONE BUCKET
88	MM=1; LL(1)=0; GO TO 55
C   MODE 2: N/5 BUCKETS
C   DISTRIBUTIVE SORTING
100	MM=N/5
	DO 44 J=1,MM
44	LL(J)=0
C   THE MAXIMAL AND MINIMAL ELEMENTS
	XMIN=A(N); XMAX=XMIN
	DO 111 I=1,N-1,2
	IF(A(I).LT.A(I+1)) 222,333
222	IF(A(I).LT.XMIN) XMIN=A(I)
	IF(A(I+1).GT.XMAX) XMAX=A(I+1)
	GO TO 111
333	IF(A(I+1).LT.XMIN) XMIN=A(I+1)
	IF(A(I).GT.XMAX) XMAX=A(I)
111	CONTINUE
	IF(XMAX.EQ.XMIN) RETURN
	AA=(MM-1)/(XMAX-XMIN); BB=1.1-AA*XMIN
C  THE HISTOGRAM
	DO 54 I=1,N
	P(I)=A(I); J=A(I)*AA+BB
54	LL(J)=LL(J)+1
C   THE INDEX OF THE BUCKETS
	DO 24 J=2,MM
24	LL(J)=LL(J-1)+LL(J)
C   THE REARRANGEMENT OF A
	DO 20 I=1,N
	J=P(I)*AA+BB; A(LL(J))=P(I); IR(LL(J))=I
20	LL(J)=LL(J)-1
C   SORTING OF THE BUCKETS
55	R=N
	DO 34 JJ=MM,1,-1
	IF(R-LL(JJ).GT.MR) 35,34
C   QUICK SORT
35	L=LL(JJ)+1; IND=1
10	I1=(L+R)/2; I2=L+1
	IV=IR(I1); IR(I1)=IR(I2); IR(I2)=IV
	V=A(I1); A(I1)=A(I2); A(I2)=V
	IF(A(I2).LE.A(R)) GO TO 2
	IV=IR(I2); IR(I2)=IR(R); IR(R)=IV
	V=A(I2); A(I2)=A(R); A(R)=V
2	IF(A(L).LE.A(R)) GO TO 3
	IV=IR(L); IR(L)=IR(R); IR(R)=IV
	V=A(L); A(L)=A(R); A(R)=V
3	I2=L+1
	IF(A(I2).LE.A(L)) GO TO 4
	IV=IR(I2); IR(I2)=IR(L); IR(L)=IV
	V=A(I2); A(I2)=A(L); A(L)=V
4	I=L+1; J=R; IV=IR(L); V=A(L)
5	I=I+1
	IF(A(I).LT.V) GO TO 5
6	J=J-1
	IF(A(J).GT.V) GO TO 6
	IF(J.LT.I) GO TO 7
	IB=IR(I); IR(I)=IR(J); IR(J)=IB
	B=A(I); A(I)=A(J); A(J)=B
	GO TO 5
7	IV=IR(L); IR(L)=IR(J); IR(J)=IV
	V=A(L); A(L)=A(J); A(J)=V
	NA=J-L; NB=R-I+1
	IF((NA.LE.M).AND.(NB.LE.M)) 11,12
11	 IF(IND.EQ.1) GO TO 34
	 IND=IND-1; L=ISA(IND); R=ISB(IND); GO TO 10
12	 IF((NA.LE.M).OR.(NB.LE.M)) 13,14
13	  IF(NA.LT.NB) 15,16
15	   L=I; GO TO 10
16	   R=J-1; GO TO 10
14	  IF(NA.LT.NB) 17,18
17	   ISA(IND)=I; ISB(IND)=R; IND=IND+1
	   R=J-1; GO TO 10
18	   ISA(IND)=L; ISB(IND)=J-1; IND=IND+1
	   L=I; GO TO 10
34	R=LL(JJ)
C   INSERTION SORT FOR THE WHOLE A
	DO 40 I=N-1,1,-1
	IF(A(I).LE.A(I+1)) GO TO 40
	IV=IR(I); V=A(I); J=I+1
50	IR(J-1)=IR(J); A(J-1)=A(J); J=J+1
	IF(A(J).LT.V) GO TO 50
	IR(J-1)=IV; A(J-1)=V
40	CONTINUE
	END
C
C   PAAOHJELMA
	DIMENSION X(20010),L(4000)
	DIMENSION IR(20010)
	DIMENSION Y(20010)
	DIMENSION P(20010)
	DSEED=123.0D0
	N=600
	DO 3 KIE=1,1
C	CALL GGEXN(DSEED,1.,N,X)
C	CALL GGNML(DSEED,N,X)
	DO 1 I=1,N
	IR(I)=I
	X(I)=RAN(Z)
	Y(I)=X(I)
1	CONTINUE
	X(N+1)=100000000.
	NP=N
	IR(N+1)=N+1
C        X(20)=1000000
	NN=N
	CALL DSOPE(X,IR,NN,NP,L,P)
C	CALL VSRTA(X,NN)
C	TYPE*,KIE
C	DO 4 I=1,N
C4	IF(Y(IR(I)).NE.X(I)) TYPE*,I,X(I),IR(I),Y(IR(I))
3	CONTINUE
	END
