PROGRAM vecgen
!
!  This program reads in a vector from the file vector.dat, with the
!  first line being the number of elements in V, and all other entries
!  being the data for V, and then sorts it from greatest-to-least,
!  and writes the sorted vector to vecsorted.dat.
!
   IMPLICIT NONE
   INTEGER N
   DOUBLE PRECISION, ALLOCATABLE :: V(:), Z(:)
   EXTERNAL :: veccopysort

   OPEN(unit=10, file='vector.dat', status='old', action='read')
   READ (10,*) N
   PRINT *, 'N read in as:', N
   IF (N .lt. 1) THEN
      PRINT*, 'N must be greater than 0, but you used:', N
      STOP
   END IF
!
!  Allocate array, read in V, and close file
!
   ALLOCATE(V(N))               ! allocate array
   ALLOCATE(Z(N))               ! allocate array
   READ(10, *) V
   CLOSE(10)
!
!  Sort the array
!
   call veccopysort(N, V, Z)
!
!  Write sorted array out, and if it is small, print vector to screen
!
   OPEN(unit=10, file='vecsorted.dat', status='replace', action='write')
   WRITE(10,*) N
   WRITE(10,*) Z
   CLOSE(10)

   IF (N .LT. 8) PRINT *, Z
   DEALLOCATE(V)
   DEALLOCATE(Z)
   PRINT *, 'Sorted vector written to file ''vecsorted.dat'''

END PROGRAM vecgen

SUBROUTINE veccopysort(N, X, Y)
!
!  This subroutine copies the unsorted array to Y, and sorts Y
!
   INTEGER, INTENT(IN) :: N
   DOUBLE PRECISION, INTENT(IN) :: X(*)
   DOUBLE PRECISION, INTENT(INOUT) :: Y(*)
   EXTERNAL :: vecsort

   Y(1:N) = X(1:N)      ! copy X to Y
   CALL vecsort(N, Y)   ! sort Y

END SUBROUTINE veccopysort

SUBROUTINE vecsort(N, Y)
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: N
   DOUBLE PRECISION, INTENT(INOUT) :: Y(*)
   INTEGER :: i, j, imax
   DOUBLE PRECISION :: mymax
!
!  sort array from greatest-to-least using selection sort
!
   DO j=1, N-1         ! j points to position to find the largest # in
      mymax = Y(j)
      imax = j
      DO i = j+1, N    ! loop over all positions to right of sorted array
         IF (Y(i) .gt. mymax) THEN
            mymax = Y(i)
            imax = i
         END IF
      END DO
      IF (imax .NE. j) THEN    ! if greatest remeaning elt not already in pos j
         Y(imax) = Y(j)        ! then swap them
         Y(j) = mymax
      END IF
   END DO

END SUBROUTINE vecsort