
      SUBROUTINE ISORTS(A,IDX,N)  ! Shell sort (integer)
*     integer keys A(*), unchanged at return 
*     N indices IDX(1) ... IDX(N), changed at return such that
*        A(IDX(L)) <= A(IDX(L+1)) 
      INTEGER I,J,M,IDXT,N,IDX(*)
      INTEGER A(*)                ! array of keys 
      M=4
 10   IF(3*M+1.LT.N) THEN    ! determine initial increment
         M=3*M+1
         GOTO 10
      END IF
 20   M=M/3                  ! next increment
      DO I=1+M,N             ! ... used in insert sort
       IDXT=IDX(I) 
       J=I
 30    IF(A(IDX(J-M)).GT.A(IDXT)) THEN
          IDX(J)=IDX(J-M)
          J=J-M
          IF(J.GT.M) GOTO 30 
       END IF
       IDX(J)=IDXT
      END DO 
      IF(M.GT.1) GOTO 20
      END

      SUBROUTINE ISORTQ(A,IDX,N)    ! quick sort (integer)
*     integer keys A(*), unchanged at return 
*     N indices IDX(1) ... IDX(N), changed at return such that
*        A(IDX(L)) <= A(IDX(L+1)) 
      INTEGER NLEV          ! stack size
      PARAMETER (NLEV=2*32) ! ... for N = 2**32 = 4.3 10**9
      INTEGER I,J,L,R,IDXT,N,IDX(N),LEV,LR(NLEV)
      INTEGER A(*),AT       ! array of keys and pivot key
      LEV=0
      L=1
      R=N
 10   IF(R-L.EQ.1) THEN     ! sort two elements
         IF(A(IDX(L)).GT.A(IDX(R))) THEN
            IDXT   = IDX(L) ! exchange
            IDX(L) = IDX(R)
            IDX(R) = IDXT  
         END IF 
         R=L 
      END IF
      IF(R.EQ.L) THEN  
         IF(LEV.LE.0) RETURN   
         LEV=LEV-2 
         L=LR(LEV+1)
         R=LR(LEV+2)
      ELSE
         AT=A(IDX((L+R)/2))
         I=L-1       ! find limits [J,I] with [L,R]
         J=R+1
 20      I=I+1
         IF(A(IDX(I)).LT.AT) GOTO 20 
 30      J=J-1
         IF(A(IDX(J)).GT.AT) GOTO 30 
         IF(I.LE.J) THEN
            IDXT   = IDX(I)
            IDX(I) = IDX(J)
            IDX(J) = IDXT
            GO TO 20
         END IF
         IF(LEV+2.GT.NLEV) STOP ' SORTQS: stack overflow '
         IF(R-I.LT.J-L) THEN 
            LR(LEV+1)=L
            LR(LEV+2)=J
            L=I
         ELSE 
            LR(LEV+1)=I
            LR(LEV+2)=R
            R=J
         END IF 
         LEV=LEV+2
      END IF
      GO TO 10
      END

      SUBROUTINE FSORTH(A,IDX,N) ! heap sort (real)
*     real keys A(*), unchanged at return 
*     N indices IDX(1) ... IDX(N), changed at return such that
*        A(IDX(L)) <= A(IDX(L+1)) 
      INTEGER I,J,L,R,IDXT,N,IDX(N) 
      REAL A(*),AT    ! array of keys and pivot key value
      L=N/2+1
      R=N
   10 IF(L.GT.1) THEN
         L=L-1
         IDXT=IDX(L)
         AT  =A(IDXT)
      ELSE
         IDXT=IDX(R)
         AT  =A(IDXT)
         IDX(R)=IDX(1)
         R=R-1
         IF(R.EQ.1) THEN
            IDX(1)=IDXT
            RETURN
         END IF
      END IF
      I=L
      J=L+L
   20 IF(J.LE.R) THEN
         IF(J.LT.R) THEN
            IF(A(IDX(J)).LT.A(IDX(J+1))) J=J+1
         END IF
         IF(AT.LT.A(IDX(J))) THEN
            IDX(I)=IDX(J)
            I=J
            J=J+J
         ELSE
            J=R+1
         END IF
         GOTO 20
      END IF
      IDX(I)=IDXT
      GOTO 10
      END

      FUNCTION LEFT(X,TAB,N) ! binary search in ordered table TAB(N)
*     (multiple table values allowed, e.g. for B-splines)
*     returned function value L with 
*        TAB(L) <= X < TAB(L+1), i.e. outlow: L=0; outhigh: L=N 
*     exception if X =TAB(N): TAB(L) <  X <= TAB(L+1) 
      INTEGER L,R,N,LR
      REAL TAB(N),X
      L=0
      R=N+1
   10 LR=(L+R)/2    ! middle index
      IF(LR.NE.L) THEN
         IF(X.GT.TAB(LR)) THEN
            L=LR     ! new left limit
         ELSE
            R=LR     ! new right limit
         END IF
         GOTO 10
      END IF
   20 IF(LR.LT.N) THEN ! check multiple keys
         IF(X.EQ.TAB(LR+1)) THEN
            LR=LR+1
            GOTO 20
         END IF
      ELSE 
         LR=N       ! check exception X=TAB(N)
   30    IF(X.LE.TAB(LR)) THEN ! multiple keys
            LR=LR-1
            GOTO 30
         END IF
      END IF
      LEFT=LR
      END

      SUBROUTINE EQUICL(NC,N,LISTJK,M) ! equivalence classes
*     Elements are numbered 1 ... N and result is in array NC(N).
*     Input is LISTJK(2,M) with M pairs of equivalent elements,  
*     for example: LISTJK(1,7)=1  LISTJK(2,7)=3 means:
*     elements 1 and 3 are equivalent. 
*     Result is in array NC(N): NC(I) is the equivalence class number
*     for element I, and all equivalent elements will have the same
*     equivalence class number.
      INTEGER NC(N),LISTJK(2,M),J,K,L,M,N
      DO K=1,N
       NC(K)=K              ! initialize classes
      END DO 
      DO L=1,M              ! loop on pairs 
       J=LISTJK(1,L)        
 10    IF(NC(J).NE.J) THEN 
          J=NC(J)
          GOTO 10
       END IF
       K=LISTJK(2,L)        
 20    IF(NC(K).NE.K) THEN
          K=NC(K)
          GOTO 20
       END IF
       IF(J.NE.K) NC(J)=K
      END DO 
      DO J=1,N              
 30    IF(NC(J).NE.NC(NC(J))) THEN
          NC(J)=NC(NC(J))
          GOTO 30
       END IF
      END DO 
      END

      SUBROUTINE MSTREE(NC,N,LISTJK,DISTJK,IDX,M) ! minimum spanning tree
*     Elements are numbered 1 ... N and result is in array NC(N).
*     Input is LISTJK(2,M) with M pairs of equivalent elements, 
*     and DISTJK(M), the non-negative distances between the elements.
*     IDX(M) is an auxiliary array used for sorting the distances.
*     Result is in array NC(N): NC(I) is the equivalence class number
*     for element I, and all equivalent elements will have the same
*     equivalence class number.
      INTEGER NC(N),IDX(M),LISTJK(2,M) ! pairs of items
      REAL DISTJK(M) ! non-negative distances between pairs
      INTEGER I,J,K,M,N
*     ...
      DO I=1,M
       IDX(I)=I                 ! prepare index vector
      END DO 
      CALL FSORTH(DISTJK,IDX,M) ! real version of heapsort
      DO J=1,N
       NC(J)=J                  ! initialize classes
      END DO
      DO I=1,M      ! loop in order of increasing distance
       J=LISTJK(1,IDX(I))       
 10    IF(NC(J).NE.J) THEN
          J=NC(J)
          GOTO 10
       END IF
       K=LISTJK(2,IDX(I))       
 20    IF(NC(K).NE.K) THEN
          K=NC(K)
          GOTO 20
       END IF
       IF(J.NE.K) THEN
          NC(J)=K               ! J and K now related
       ELSE 
          DISTJK(I)=-DISTJK(I)-1.0  ! flag as removed
       END IF
      END DO
      DO J=1,N                  
 30      IF(NC(J).NE.NC(NC(J))) THEN
          NC(J)=NC(NC(J))
         GOTO 30
       END IF
      END DO
      END







