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