Fixing and adding unit test coverage to merge_sort indices permutation results.

radixalgorithm/mergesort.f90

@@ -15,7 +15,12 @@
     real, intent(inout) :: a(n)
     integer, intent(out) :: indices(n)
     real, allocatable :: b(:)
+    integer, allocatable :: ai(:)
     allocate (b(n))
+    allocate (ai(n))
+    do i=1, n
+      ai(i) = i
+    end do
     L1 = 1
     m  = 1
     do while (m < n)    ! Determine L1 so that 2**(L1-1) < n < 2**L1
@@ -36,32 +41,34 @@
         do while(i < imax  .and.  j < jmax)
           if(a(i) < a(j)) then
             b(k) = a(i)
-            indices(k) = i
+            indices(k) = ai(i)
             i = i + 1
           else
             b(k) = a(j)
-            indices(k) = j
+            indices(k) = ai(j)
             j = j + 1
           end if
           k = k + 1
         end do
         do while(i < imax)
           b(k) = a(i)
-          indices(k) = i
+          indices(k) = ai(i)
           i = i + 1
           k = k + 1
         end do
         do while(j < jmax)
           b(k) = a(j)
-          indices(k) = j
+          indices(k) = ai(j)
           j = j + 1
           k = k + 1
         end do
       end do
       m = m + m
       a = b
+      ai = indices
     end do
     if(allocated(b)) deallocate(b)
+    if(allocated(ai)) deallocate(ai)
     return
     end subroutine merge_sort
 

radixalgorithm/tests/tstMergeSort.f90

-    Program run_merge
-    implicit none
-    interface
-      subroutine merge_sort(n, a, indices)
-      integer :: n
-      real, intent(inout) :: a(n)
-      integer, intent(out) :: indices(n)
-      end subroutine
-    end interface
-    integer :: n, i
-    character (len=10) :: d1, d2, d3
-    integer :: time(8)
-    real :: dtime1, dtime2
-    real, allocatable :: a(:)
-    integer, allocatable :: indices(:)
-    call DATE_AND_TIME(d1, d2, d3, time)
-    dtime1 = time(6)*60.+time(7)+time(8)/1000.
-    n = 10000
-    allocate (a(n))
-    allocate (indices(n))
-    ! avoid unitialised value warning
-    indices = 0
-    do i=1, time(7)
-      call RANDOM_NUMBER(dtime2)
-    end do
-    call RANDOM_NUMBER(a)
-    do i=1, n
-      a(i) = a(i) * n / 10.
-    end do
-!    print "(10f12.3)", (a(i,1), i=1, 10), (a(i,2), i=1, 10)
-!    a = a * N / 10.
-    call DATE_AND_TIME(d1, d2, d3, time)
-    dtime2 = time(6)*60.+time(7)+time(8)/1000.
-    dtime1 = dtime2 - dtime1
-    print *, "setup time (s) =", dtime1
-    print *, "First/Last 10 elements before sorting"
-    print "(10f12.3/10i10)", (a(i), i=1, 10), (indices(i), i=1, 10)
-    print "(10f12.3/10i10)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
-    call merge_sort(n, a, indices)
-    call DATE_AND_TIME(d1, d2, d3, time)
-    print *, "First/Last 10 elements after sorting"
-    print "(10f12.3/10i10)", (a(i), i=1, 10), (indices(i), i=1, 10)
-    print "(10f12.3/10i10)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
-    dtime2 = time(6)*60.+time(7)+time(8)/1000.-dtime2
-    print *, "sort time (s) =", dtime2
-    do i=1, (n-1)
-      ! check order of entire array
-      if(a(i) > a(i+1)) then
-         print *, "Element (",i,") is greater than element (",(i+1),")"
-         stop 99
-      endif
-    end do
-    stop
-    end program run_merge
+    Program run_merge
+    implicit none
+    interface
+      subroutine merge_sort(n, a, indices)
+      integer :: n
+      real, intent(inout) :: a(n)
+      integer, intent(out) :: indices(n)
+      end subroutine
+    end interface
+    integer :: n, i
+    character (len=10) :: d1, d2, d3
+    integer :: time(8)
+    real :: dtime1, dtime2
+    real, allocatable :: a(:)
+    real, allocatable :: blessed(:)
+    integer, allocatable :: indices(:)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    dtime1 = time(6)*60.+time(7)+time(8)/1000.
+    n = 1000
+    allocate (a(n))
+    allocate (blessed(n))
+    allocate (indices(n))
+    ! avoid unitialised value warning
+    indices = 0
+    do i=1, time(7)
+      call RANDOM_NUMBER(dtime2)
+    end do
+    call RANDOM_NUMBER(a)
+    do i=1, n
+      a(i) = a(i) * n / 10.
+    end do
+    ! save this original version of a as blessed version
+    blessed = a
+!    print "(10f12.3)", (a(i,1), i=1, 10), (a(i,2), i=1, 10)
+!    a = a * N / 10.
+    call DATE_AND_TIME(d1, d2, d3, time)
+    dtime2 = time(6)*60.+time(7)+time(8)/1000.
+    dtime1 = dtime2 - dtime1
+    print *, "setup time (s) =", dtime1
+    print *, "First/Last 10 elements before sorting"
+    print "(10f12.3/10i12)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10i12)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    call merge_sort(n, a, indices)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    print *, "First/Last 10 elements after sorting"
+    print "(10f12.3/10i12)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10i12)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    dtime2 = time(6)*60.+time(7)+time(8)/1000.-dtime2
+    print *, "sort time (s) =", dtime2
+    do i=1, (n-1)
+      ! check order of entire array
+      if(a(i) > a(i+1)) then
+         print *, "Element (",i,") is greater than element (",(i+1),")"
+         stop 99
+      endif
+    end do
+    ! check that the lookup works now
+
+    !print *, "unsorted:(10f3.2)", (blessed(i), i=1, n)
+    !print *, "Indices:(10i5)", (indices(i), i=1, n)
+    !print *, "sorted:(10f3.2)", (a(i), i=1, n)
+    do i=1, n
+      if(blessed(indices(i)) /= a(i)) then
+         print *, "Index (",i,") doesn't not match unsorted(",blessed(indices(i)),") with sorted(",a(i),")"
+         stop 100
+      endif
+    end do
+    stop
+    end program run_merge