C **********************************************************************
C SHELK.F
C                                                             07.06.80 *        
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
C=*                                                                    *
C=* This program is free software; you can redistribute it and/or      *
C=* modify it under the terms of the GNU General Public License as     *
C=* published by the Free Software Foundation; either version 2 of the *
C=* License, or (at your option) any later version.                    *
C=*                                                                    *
C=* This program is distributed in the hope that it will be useful,    *
C=* but WITHOUT ANY WARRANTY; without even the implied warranty of     *
C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C
C  PURPOSE:  SORT VECTOR X() IN ASCENDING ORDER                                  *        
C            WARNING: ORIGINAL VECTOR(X) GETS CLOBBERED UPON RETURN         
C            BUT THE ORIGINAL LOCATIONS ARE SAVED IN VECTOR KX()                 *        
C            REFERENCES:
C            (1) J.BOOTHROYD  SHELLSORT ALGORITHM.201 
C                COMM.ACM  VOL.6  (1963),NO.8,PP.445   
C            (2) D.A.SHELL  A HIGH-SPEED SORTING PROCEDURE         
C                COMM.ACM  VOL.2(1959),PP.30-32                                       *        
C                                                                    
C*---------------------------------------------------------------------*        

      SUBROUTINE SHELK(N, X, KX )
                                            
C     I DO NOT KNOW IF SAVE IS NEEDED FEB 99 al
      SAVE

      DIMENSION  X(N) , KX(N) 
                                                  
      DO   J = 1,N                                                          
         KX(J) = J                                                                 
      ENDDO

      I     = 1                                                                 
   20 I     = I + I                                                             
      IF (I .LE. N) GO TO  20                                          
      M     = I - 1                                                             
   30 M     = M / 2                                                             
      IF (M .EQ. 0) GO TO  70                                          
      K     = N - M  
                                                           
      DO 60 J = 1,K                                                          
         JM = J + M                                                             
   40    JM = JM - M                                                            
         IF (JM .LE. 0) GO TO  60                                          
         L = JM + M                                                            
         IF (X(L) .GE. X(JM)) GO TO  60                                          
         PIV   = X(JM)                                                             
         X(JM) = X(L)                                                              
         X(L)  = PIV                                                               
         KPIV  = KX(JM)                                                            
         KX(JM)= KX(L)                                                             
         KX(L) = KPIV                                                              
         GO TO  40 
                                         
   60   CONTINUE                                                                
        GO TO  30 
                                         
   70   RETURN                                                                  
        END