#### Topic: Merge Sort Routine

Hi Jeff,

While in lock-down and to keep my mind active, I've written a sorting routine for real arrays. I'm posting the sort routine here.  Although it's far from being the fastest sort routine, it's simple. I've sorted up to 65,000 random numbers with the routine with no problem. However, speed wise, it's more practical in the range of 3 to 4000 numbers (e.g., Random numbers =1034, t = 0.04 sec.; Random numbers =4108, t = 0.4 sec.).  Like many routines, sometimes we learn that similar logic has been applied before and are not aware of it. The program is self explanatory and can be compiled as is and run with SimplyFortran.

I would welcome any constructive suggestions from forum members on how to improve the routine.

Thanks,
Frank

!-------------------------------------------------------------------
!----------------- Start of MergSort Routine -----------------------
!-------------------------------------------------------------------
MODULE nrtype
IMPLICIT NONE

!--- Listing valid kind parameters types
INTEGER, PARAMETER :: I1B = KIND(2)
INTEGER, PARAMETER :: I2B = KIND(4)
INTEGER, PARAMETER :: I4B = KIND(9)

INTEGER, PARAMETER :: SP = KIND(1.0E0)
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 200)
INTEGER, PARAMETER :: QP = SELECTED_REAL_KIND(2*PRECISION(1.0_DP))

INTEGER, PARAMETER :: LGT = KIND(.TRUE.)
END MODULE nrtype

!-------------------------------------------------------------------
Module MSort
USE nrtype
Implicit None

PRIVATE

PUBLIC :: MergSort

INTERFACE MergSort
MODULE PROCEDURE MergSort
END INTERFACE

!-------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------

!---------------------------------------------------------------
FUNCTION MergSort( array ) RESULT( sorted )
!---------------------------------------------------------------
!---    MergSort version 1.0
!---
!---    Developed by F.W. Perrella, Ph.D.  November 27, 2020
!---            Silver Spring, Maryland U.S.A.
!---------------------------------------------------------------
! NAME:
!       MergSort
!
! PURPOSE:
!       Program to sort a numeric array in ascending order.
!
! MODULE:
!       nrType: Definition for kind variable types.
!---------------------------------------------------------------
!
!    CALLING SEQUENCE:
!       Result = MergSort( array )
!
!    INPUT ARGUMENTS:
!       array:      Real array of numbers to sort.
!       TYPE:       REAL( KIND = Double )
!
!    FUNCTION RESULT:
!       TYPE:    Real( double )
!
!       RESULT: A real array of values sorted in ascending order.
!---------------------------------------------------------------
!
!   BASIC MERGSORT ROUTINE:
!     Function MergSort(array) Result(sorted)
!       c = 123456789.0
!       unsorted = array
!       sorted = 0.0
!       Do i = 1, n
!         j = MINLOC(unsorted,1)
!         vecMin = unsorted(j)
!         unsorted = MERGE(unsorted,c,MASK=(unsorted /= vecMin))
!         sorted(i) = vecMin
!       End Do
!     End Function
!---------------------------------------------------------------
use nrtype
implicit none
real(DP), dimension(:), intent(in) :: array
real(DP), dimension(1:size(array)) :: sorted

!--- Local orking arrays
real(DP), allocatable, dimension(:) :: unsorted
real(DP), allocatable, dimension(:) :: c

!--- Local variables
real(DP)                            :: vecMin
integer(I4B)                        :: i, j, n

!--- Number of values to sort
n = SIZE( array )

!--- Array too small to be sorted
if ( n == 1 ) then
sorted = array
return
endif

!--- Allocate working arrays
allocate( unsorted(n), c(n) )

!--- Large number array place holder
c = 123456789.0

!--- Begin sorting array
unsorted = array
sorted = 0.0

do i = 1, n
!--- Find minimum value in array
j = MINLOC( unsorted, 1 )
vecMin = unsorted(j)

!--- Merge array eliminating minimum value
unsorted = MERGE( unsorted, c, MASK=(unsorted /= vecMin) )

!--- Save the minimum values in array
sorted(i) = vecMin
end do

!--- Deallocate working arrays
Deallocate( unsorted, c )

return
End Function MergSort

End Module MSort

!-------------------------------------------------------------------
PROGRAM Test_Sort
use nrtype
implicit none
character(1) :: a

call Test_MergSort()

WRITE(*,'(/1x,"Press return to exit",\$)')
READ(*,'(a1)') a
End Program Test_Sort

!-------------------------------------------------------------------
Subroutine Test_MergSort()
use nrtype
use MSort
implicit none

real(SP) :: t1, t2, tsum
real(DP),allocatable :: array1(:)
integer(I4B) :: j, k, m, n
integer(I4B) :: left, right, isize
logical(LGT) :: bflag, btot

m = 10
n = 14
btot = .TRUE.

WRITE(*,'(/1x,"Testing MergSort")')
WRITE(*,'(1x,a)') REPEAT("-", 35)

do j = m, n
!--- 1,034 to 16,398 random numbers
isize = 2**j
isize = isize + j
allocate( array1(isize) )

tsum = 0.0
do k = 1, m
call random_number( array1 )

left = 1
right = isize

!--- Start timer
t1 = DeltaTime( .TRUE. )

!--- Sort in Ascending Order
array1 = MergSort( array1 )

!--- End timer
t2 = DeltaTime( .FALSE. )

tsum = tsum + (t2 - t1)
end do

WRITE(*,'(1x,"Random numbers =",1x,I0)') isize
WRITE(*,'(1x,"Average elapsed time =",1x,F6.4,1x,"sec")') tsum/m

!--- Check that the array was sorted correctly.
call check_sort( array1, left, right, bflag )

btot = (btot .eqv. bflag)

deallocate(array1)
end do

WRITE(*, &
'(1x,"All tests were sorted correctly:",1x,L1,"RUE")') btot
WRITE(*,'(1x,a)') REPEAT("-", 35)

return
!End Subroutine Test_MergSort

!-------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------

!-------------------------------------------------------------------
Subroutine check_sort( array, left, right, bflag )
use nrtype
implicit none
integer(I4B), intent(in)  :: left, right
real(DP), intent(inout)   :: array(left:right)
logical(LGT), intent(out) :: bflag
integer(I4B) :: j

bflag = .true.
do j = 1, size(array)-1
if( array(j+1) < array(j) ) then
bflag = .false.
EXIT
end if
end do

if (bflag) then
print *, "Correctly sorted array"
else
print *, "Incorrectly sorted array"
end if
PRINT *

return
End Subroutine check_sort

!-------------------------------------------------------------------
FUNCTION DeltaTime(Start) RESULT(tm)
!-------------------------------------
!--- Routine to calculate run time
!    print *, 'Starttime:',  t0
!    print *, 'Finishtime:', t1
!    print *, 'Elapsed time=', t1 - t0
!-------------------------------------
IMPLICIT NONE
LOGICAL, INTENT(IN) :: Start
LOGICAL             :: Begin
REAL, SAVE          :: t0
REAL                :: t1
REAL                :: tm
REAL, SAVE          :: tt

!--- Start timer
!--- Start = .TRUE.
Begin = Start
IF ( Begin ) THEN
call cpu_time ( t0 )
t1 = t0
Begin = .FALSE.

!--- Stop timer
ELSE !--- Start = .FALSE.
call cpu_time ( t1 )
END IF

!--- Calculate time in seconds
tt = tt + (t1 - t0)

!--- Return total time in seconds
tm = tt
END FUNCTION DeltaTime

!-------------------------------------------------------------------
End Subroutine Test_MergSort
!-------------------------------------------------------------------
!---------------- End of MergSort Routine --------------------------
!-------------------------------------------------------------------

#### Re: Merge Sort Routine

MergSort Update:

The following revised sort routine is x2 faster than the previously posted routine.
I realize that it's simple, but it works fairly well.

Frank

REVISED BASIC MERGSORT ROUTINE:
--------------------------------
Function MergSort(array) Result(sorted)
n = SIZE( array )
unsorted = array
sorted = 0.0
c = 123456789.0
vecMin = c
i = 0
Do While( ANY(unsorted < c) )
j = MINLOC( unsorted, 1 )
vecMin = unsorted(j)
unsorted(j) = c
i = i + 1
sorted(i) = vecMin
End Do
End Function

Listing of the Revised MergSort Module below:

Module MSort
USE nrtype
Implicit None

PRIVATE

PUBLIC :: MergSort

INTERFACE MergSort
MODULE PROCEDURE MergSort
END INTERFACE

CONTAINS

function MergSort( array ) result( sorted )
!---------------------------------------------------------------
!---    MergSort version 1.1
!---
!---    Developed by F.W. Perrella, Ph.D.  November 28, 2020
!---            Silver Spring, Maryland U.S.A.
!---------------------------------------------------------------
! NAME:
!       MergSort
!
! PURPOSE:
!       Program to sort a numeric array in ascending order.
!
! MODULE:
!       nrType: Definition for kind variable types.
!---------------------------------------------------------------
!
!    CALLING SEQUENCE:
!       Result = MergSort( array )
!
!    INPUT ARGUMENTS:
!       array:      Real array of numbers to sort.
!       TYPE:       REAL( KIND = Double )
!
!    FUNCTION RESULT:
!       TYPE:    Real( double )
!
!       RESULT: A real array of values sorted in ascending order.
!---------------------------------------------------------------
!---
!--- Execution Time: 0.003 sec. to sort 2000 random numbers
!---
use nrtype
implicit none
real(dp), dimension(:), intent(in)  :: array
real(dp), dimension(1:size(array))  :: sorted

!--- Local orking arrays
real(dp), allocatable, dimension(:) :: unsorted

!--- Local variables
integer(I4B)                        :: i, j, n
real(dp)                            :: vecMin
real(dp), parameter                 :: c = 123456789.0

!--- Number of values to sort
n = SIZE( array )

!--- Array too small to be sorted
if ( n == 1 ) then
sorted = array
return
endif

!--- Begin sorting array
unsorted = array
sorted = 0.0
vecMin = c
i = 0

DO WHILE( ANY(unsorted < c) )
!--- Find minimum value in array
j = MINLOC( unsorted, 1 )

!--- Minimum value
vecMin = unsorted(j)

!--- Replace array minimum with a very large number.
unsorted(j) = c

!--- Save the minimum values in array
i = i + 1
sorted(i) = vecMin
END DO

!--- Deallocate working arrays
Deallocate( unsorted )

return
end function MergSort

End Module MSort

#### Re: Merge Sort Routine

Here is a revised Floating Point sort (FPSORT) routine that sorts x2.5 faster than the previously posted sort routine and sorts in ascending or descending order. It utilizes MAXLOC and reduces the size (PACK) of the searched working array as the sort progresses.

!--- FPSort ROUTINE:
Function FPSort(array,bflag) Result(sorted)
bFlag = .true.
n0 = SIZE( array )
n = n0
n1 = n / 40
unsorted = array
sorted = 0.0
c = -HUGE(1.0)
vecMax = c
i = 0

DO WHILE( ANY(unsorted(:n) > c) )
j = MAXLOC( unsorted(:n), 1 )
vecMax = unsorted(j)
unsorted(j) = c
i = i + 1
sorted(i) = vecMax
IF ( (n > 400).AND.(MODULO(i, n1) == 0) ) &
unsorted = PACK(unsorted(:n), unsorted(:n) /= c)
n = SIZE( unsorted )
End Do

If (bFlag) Then
!--- Reverse order Ascending Order Sort
sorted = sorted(n0:1:-1)
End If
End Function FPSort

The complete listing is shown below.

Happy Holidays,
Frank

Module FPSort_module
Implicit None

PRIVATE

PUBLIC :: FPSort

INTERFACE FPSort
MODULE PROCEDURE FPSort
END INTERFACE

CONTAINS

Function FPSort( array,bFlag ) Result( sorted )
!---------------------------------------------------------------
!---    FPSort version 2.0
!---    (x4 Faster than a Bubble sort routine)
!---
!---    Developed by F.W. Perrella, Ph.D.  December 7, 2020
!---    Silver Spring, Maryland U.S.A.
!---------------------------------------------------------------
! NAME:
!       FPSort
!
! PURPOSE:
!       Sorts a numeric array in ascending or descending order.
!
! MODULE:
!       nrType: Definition for kind variable types.
!---------------------------------------------------------------
!
!    CALLING SEQUENCE:
!       Result = FPSort( array,bFlag )
!
!    INPUT ARGUMENTS:
!       array: Real array of numbers to sort
!       TYPE:  REAL( KIND = Double )
!
!       bFlag: .TRUE.= Ascending sort, .FALSE.= Descending sort
!       TYPE:  LOGICAL
!
!    FUNCTION RESULT:
!       TYPE:    Real( Double )
!
!       RESULT: Real array of values sorted (Ascending,Descending)
!---------------------------------------------------------------
use nrtype
implicit none
logical(LGT), intent(in)                  :: bFlag
real(DP), dimension(:), intent(in)   :: array
real(DP), dimension(1:size(array)) :: sorted

!--- Local working array
real(DP), allocatable, dimension(:) :: unsorted

!--- Local variables
integer(I4B)                        :: i, j
integer(I4B)                        :: n, n0, n1
real(DP)                              :: vecMax
real(DP), parameter             :: c = -HUGE(1.0)

!--- Number of values to sort
n = SIZE( array )

!--- Allocate working array
IF (ALLOCATED( unsorted )) DEALLOCATE( unsorted )
ALLOCATE( unsorted(1:n) )

!--- Save n to n0
n0 = n

!--- Array too small to be sorted
IF ( n0 == 1 ) THEN
sorted = array
GO TO 900
!return
END IF

!--- Initialize variable
sorted = array
vecMax = c
n = n0
n1 = n
i = 0

!--- Working array size reduction factor
n1 = n0 / 40

!--- Working array
unsorted(:n) = sorted(:n)

!--- Begin sorting array
!--- Sort 1 to n element array
DO WHILE( ANY(unsorted(:n) > c) )
!--- Find maximum value in the array
j = MAXLOC( unsorted(:n), 1 )

!--- Maximum value
vecMax = unsorted(j)

!--- Replace array maximum with a huge negative number
unsorted(j) = c

!--- Save the maximum values in the array
i = i + 1
sorted(i) = vecMax

!--- Reduce the size of the working array
IF ( (n > 400).AND.(MODULO(i, n1) == 0) ) &
unsorted = PACK(unsorted(:n), unsorted(:n) /= c)

!--- Revised working array size
n = SIZE( unsorted )
END DO

!--- Ascending (.true.) or Descending (.false.) Sort Order
IF (bFlag) THEN
!--- Ascending Order Sort
!--- Reverse the order of elements in descending order array
sorted = sorted(n0:1:-1)
ELSE !--- Descending Order Sort
CONTINUE
END IF

900     CONTINUE

!--- Deallocate working array
Deallocate( unsorted )

RETURN
End Function FPSort

End Module FPSort_module

#### Re: Merge Sort Routine

Revised and updated FPSort that's faster than the previous post using pointers and reducing the size of the working array during the search.  FPSort can sort in both ascending and descending order.

Happy Holidays,
Dr.Frank

Fortran Listing:

MODULE FPSort_module
IMPLICIT NONE

PRIVATE

PUBLIC :: FPSort

INTERFACE FPSort
MODULE PROCEDURE FPSort
END INTERFACE

CONTAINS

FUNCTION FPSort( Array, bFlag ) RESULT( sorted )
!---------------------------------------------------------------
!---    FPSort version 2.0
!---    (Four-times faster than a bubble sort of random numbers)
!---
!---    Developed by F.W. Perrella, Ph.D.  December 13, 2020
!---    Silver Spring, Maryland U.S.A.
!---------------------------------------------------------------
! NAME:
!       FPSort
!
! PURPOSE:
!       Sorts a numeric array in ascending or descending order.
!
! MODULE:
!       nrType: Definition for kind variable types.
!       INTEGER, PARAMETER :: I4B = KIND(9)
!       INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 200)
!
!---------------------------------------------------------------
!
!    CALLING SEQUENCE:
!      Result = FPSort( Array,bFlag )
!
!    INPUT ARGUMENTS:
!      Array: Real array of numbers to sort
!      TYPE:  REAL( KIND = Double )
!
!      bFlag: .TRUE.= Ascending sort, .FALSE.= Descending sort
!      TYPE:  LOGICAL
!
!    FUNCTION RESULT:
!      TYPE:    Real( Double )
!      RESULT: Real array of values sorted (Ascending,Descending)
!---------------------------------------------------------------
!    FPSort:
!
!    Random numbers =521
!    average time=   0.0000
!
!    Random numbers =1034
!    average time=   0.0000
!
!    Random numbers =2059
!    average time=   0.001563
!
!    Random numbers =4108
!    average time=   0.004688
!
!    Random numbers =8205
!    average time=   0.02656    (x3.9 faster than Bubble Sort)
!
!    Random numbers =16398
!    average time=  0.1016      (x4.1 faster than Bubble Sort)
!
!    Random numbers =32783
!    average time=  0.3906      (x4.4 faster than Bubble Sort)
!
!    Random numbers =65552
!    average time=   1.4547     (x4.1 faster than Bubble Sort)
!
!---------------------------------------------------------------
USE nrtype
IMPLICIT NONE
LOGICAL(LGT), INTENT(IN)            :: bFlag
REAL(DP), DIMENSION(:), INTENT(IN)  :: Array
REAL(DP), DIMENSION(1:SIZE(Array))  :: sorted

!--- Local working array
TYPE :: Data_Type
REAL(DP) :: rval
END TYPE Data_Type
TYPE (Data_Type), DIMENSION(:), ALLOCATABLE, TARGET :: Array2

!--- Pointer to local working array
TYPE :: BOX
TYPE (Data_Type), POINTER :: Ptr
END TYPE BOX
!--- A()%Ptr is a pointer to Data_Type.
TYPE (BOX), DIMENSION(:), ALLOCATABLE :: A

!--- Local variables
INTEGER(I4B)                        :: i, j
INTEGER(I4B)                        :: n, n0

!--- Number of values to sort
n = SIZE( Array )

!--- Save n in n0
n0 = n

!--- Allocate working arrays
ALLOCATE( Array2(1:n), A(1:n) )

!--- Array too small to be sorted
IF ( n == 1 ) THEN
sorted(:) = Array(:)
GO TO 900
!return
END IF

!--- Save data in working array
Array2(:)%rval = Array(:)

!--- Pointer to working array
FORALL (i=1:n) A(i)%Ptr => Array2(i)

!--- Initialize variables
sorted = 0.0
i = 0
j = 0

!--- Begin sorting 1 to n element array
DO WHILE( n > 0 )
!--- Find maximum value in the array
j = MAXLOC( Array2(:n)%rval, 1 )

!--- Save the maximum values in the sorted array
i = i + 1
sorted(i) = A(j)%Ptr%rval

!--- Assign last array element to current maximum location
Array2(j)%rval = A(n)%Ptr%rval

!--- Reduce the size of the working array
!--- eliminating the last element.
n = n - 1
END DO

!--- Ascending (.true.) or Descending (.false.) Sort Order
IF (bFlag) THEN
!--- Ascending Order Sort
!--- Reverse the order of elements in descending order array
sorted = sorted(n0:1:-1)
ELSE !--- Descending Order Sort
CONTINUE
END IF

900     CONTINUE

!--- Deallocate working array
DEALLOCATE( Array2 )

RETURN
END FUNCTION FPSort

END MODULE FPSort_module

#### Re: Merge Sort Routine

Frank,

I keep meaning to spend some time trying this out!  Sorry I've been quiet!

Jeff Armstrong
Approximatrix, LLC