Adding Fortran mergesort.f90 to radixalgorithm.

radixalgorithm/CMakeLists.txt

@@ -11,6 +11,11 @@ marchingsquares.hh
 marchingsquares.i.hh
 )
 
+IF(${PROJECT_NAME}_ENABLE_Fortran)
+   SET(SOURCE ${SOURCE}
+       mergesort.f90)
+ENDIF()
+
 TRIBITS_ADD_LIBRARY(radixalgorithmlib
    SOURCES ${SOURCE}
    NOINSTALLHEADERS ${HEADERS}

radixalgorithm/mergesort.f90

+    subroutine merge_sort(n, a, indices)
+! This routine takes an input array 'a' of dimension (n)
+! and an input array 'indices' of dimension (n)
+! and conducts a merge sort on array 'a' while saving the original order
+! to array 'indices'
+
+! The first pass sorts each pair of consecutive elements [a(i) vs. a(i+1)].
+! The second pass sorts adjacent sequences of length 2 [a(i), a(i+1) vs. a(i+2), a(i+3)]; the sequences are already ordered, due to step 1.
+! The third pass sorts adjacent sequences of length 4 [a(i),...,a(i+3) vs. a(i+4),...,a(i+7)]; the sequences are already ordered, from step 2.
+! This continues for L1 passes, where 2**(L1-1) < n < 2**L1.  Note that n does not need to be a power of 2.
+
+    implicit none
+    integer :: i, imax, i0, j, jmax, k, L, L1, m, n
+    real    :: a(n)
+    real    :: indices(n)
+    real, allocatable :: b(:)
+    allocate (b(n))
+    L1 = 1
+    m  = 1
+    do while (m < n)    ! Determine L1 so that 2**(L1-1) < n < 2**L1
+      m = m + m
+      L1 = L1 + 1
+    end do
+    L1 = L1 - 1
+    m = 1
+    do L=1, L1
+      k=1
+      do i0=1, n-m+1, m+m
+        i=i0
+        j=i+m
+        imax=j
+        jmax=j+m
+        if(imax > n) imax = n + 1
+        if(jmax > n) jmax = n + 1
+        do while(i < imax  .and.  j < jmax)
+          if(a(i) < a(j)) then
+            b(k) = a(i)
+            indices(k) = i
+            i = i + 1
+          else
+            b(k) = a(j)
+            indices(k) = j
+            j = j + 1
+          end if
+          k = k + 1
+        end do
+        do while(i < imax)
+          b(k) = a(i)
+          indices(k) = i
+          i = i + 1
+          k = k + 1
+        end do
+        do while(j < jmax)
+          b(k) = a(j)
+          indices(k) = j
+          j = j + 1
+          k = k + 1
+        end do
+      end do
+      m = m + m
+      a = b
+    end do
+    if(allocated(b)) deallocate(b)
+    return
+    end subroutine merge_sort
+

radixalgorithm/tests/CMakeLists.txt

@@ -3,3 +3,9 @@ INCLUDE(GoogleTest)
 ADD_GOOGLE_TEST(tstOrdering.cc NP 1)
 ADD_GOOGLE_TEST(tstMarchingSquares.cc NP 1)
 ADD_GOOGLE_TEST(tstChaikins.cc NP 1)
+IF(${PROJECT_NAME}_ENABLE_Fortran)
+TRIBITS_ADD_EXECUTABLE_AND_TEST(tstMergeSortFortran
+    SOURCES tstMergeSort.f90
+    LINKER_LANGUAGE Fortran
+    )
+ENDIF()

radixalgorithm/tests/tstMergeSort.f90

+    Program run_merge
+    implicit none
+    integer :: n, i
+    character (len=10) :: d1, d2, d3
+    integer :: time(8)
+    real :: dtime1, dtime2
+    real, allocatable :: a(:)
+    real, allocatable :: indices(:)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    dtime1 = time(6)*60.+time(7)+time(8)/1000.
+    n = 100000000
+    allocate (a(n))
+    allocate (indices(n))
+    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/10f12.0)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10f12.0)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    call merge_sort(n, a, indices)
+    print "(10f12.3/10f12.0)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10f12.0)", (a(i), i=n-9, n), (indices(i), i=n-9, n)
+    call DATE_AND_TIME(d1, d2, d3, time)
+    print *, "First/Last 10 elements after sorting"
+    print "(10f12.3/10f12.0)", (a(i), i=1, 10), (indices(i), i=1, 10)
+    print "(10f12.3/10f12.0)", (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