C   SUBROUTINE DSORT(A,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 ELEMENT A(N+1)
C   MUST CONTAIN A VALUE GREATER THAN OR EQUAL THE OTHER VALUES OF A.
C   CALLING OF THE PROGRAM:
C
C      CALL DSORT(A,N,NP,LL,P)
C
C   WHERE
C      A=  THE ARRAY TO BE SORTED (INPUT/OUTPUT); LENGTH N+1
C      N=  NUMBER OF ELEMENTS TO BE SORTED (INPUT)
C      NP= A SWITCH (INPUT)
C      LL= AN AUXILIARY ARRAY (INPUT)
C      P= AN AUXILIARY ARRAY (INPUT).
C
C      THE METHOD
C   THE SUBROUTINE HAS TWO OPERATING MODES,
C   RULED BY THE VALUES OF THE PARAMERERS N AND NP.
C   MODE 1: IF N IS LESS THAN 500 OR NP IS LESS THAN N, QUICK
C      SORT IS USED.
C   MODE 2: OTHERWISE, THE METHOD IS DISTRIBUTIVE SORTING IN WHICH
C      THE ELEMENTS OF A ARE GROUPED INTO BUCKETS AND THE
C      BUCKETS ARE SORTED BY QUICK SORT OR INSERTION SORT.
C   IF THE USER WANTS THAT THE PROGRAM WORKS IN MODE 1, HE
C   CAN CALL THE SUBROUTINE WITH THE VALUE NP=0 AND THEN
C   THE CALLING PROGRAM MUST DEFINE THE ARRAYS LL(1) AND P(1)
C   IF THE USER WANTS THAT THE PROGRAM OPERATES IN MODE 2, HE
C   MUST CALL IT WITH NP GREATER THAN EQUAL TO N AND THEN HAVE
C   IN THE CALLING PROGRAM A DEFINITION OF THE ARRAYS LL AND P
C   AT LEAST OF THE LENGTH LL(NP/5) AND P(NP), RESPECTIVELY.
C   A(N+1) MUST CONTAIN A VALUE GREATER THAN ANY OF THE
C   ELEMENTS A(1),..,A(N).
C
C   FOR FURTHER INFORMATION OF THE METHOD, SEE
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 AURHORS:
C      OLLI NEVALAINEN
C      DEPARTMENT OF COMPUTER SCIENCE,
C      UNIVERSITY OF TURKU,
C      SF-20500 TURKU
C      FINLAND
C
C   DATE: 30.11.1981
C
C
	SUBROUTINE DSORT(A,N,NP,LL,P)
	DIMENSION A(1),LL(1),P(1)
	DIMENSION IA(50),IB(50)
	INTEGER R
	DATA MR/14/,M/10/
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 THE A-ARRAY
	DO 20 I=1,N
	J=P(I)*AA+BB; A(LL(J))=P(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	V=A((L+R)/2); A((L+R)/2)=A(L+1); A(L+1)=V
	IF(A(L+1).LE.A(R)) GO TO 2
	V=A(L+1); A(L+1)=A(R); A(R)=V
2	IF(A(L).LE.A(R)) GO TO 3
	V=A(L); A(L)=A(R); A(R)=V
3	IF(A(L+1).LE.A(L)) GO TO 4
	V=A(L+1); A(L+1)=A(L); A(L)=V
4	I=L+1; J=R; 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
	B=A(I); A(I)=A(J); A(J)=B
	GO TO 5
7	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=IA(IND); R=IB(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	   IA(IND)=I; IB(IND)=R; IND=IND+1
	   R=J-1; GO TO 10
18	   IA(IND)=L; IB(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
	V=A(I); J=I+1
50	A(J-1)=A(J); J=J+1
	IF(A(J).LT.V) GO TO 50
	A(J-1)=V
40	CONTINUE
	END
C
C   PAAOHJELMA
	DIMENSION X(20010),L(4000)
C	DIMENSION Y(20010)
	DIMENSION P(20010)
	DSEED=123.0D0
	N=5000
	DO 3 KIE=1,100
C	CALL GGEXN(DSEED,1.,N,X)
C	CALL GGNML(DSEED,N,X)
	DO 1 I=1,N
	X(I)=RAN(Z)
C	Y(I)=X(I)
1	CONTINUE
	X(N+1)=100000000.
	NP=N
C        X(20)=1000000
	NN=N
	CALL DSORT(X,NN,NP,L,P)
C	CALL VSRTA(X,NN)
C	TYPE*,KIE
C	DO 4 I=1,N
C	DO 5 J=1,N
C5	IF(Y(I).EQ.X(J)) GO TO 4
C	TYPE *,Y(I),I
C4	CONTINUE
3	CONTINUE
	END