Skip to content

Commit 8b3d48c

Browse files
Add integer array sorting
1 parent 48fb6a3 commit 8b3d48c

File tree

4 files changed

+361
-0
lines changed

4 files changed

+361
-0
lines changed

src/linalg.f90

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3458,6 +3458,8 @@ module linalg
34583458
module procedure :: sort_cmplx_array_ind
34593459
module procedure :: sort_eigen_cmplx
34603460
module procedure :: sort_eigen_dbl
3461+
module procedure :: sort_int32_array
3462+
module procedure :: sort_int32_array_ind
34613463
end interface
34623464

34633465
!> @brief Computes the LQ factorization of an M-by-N matrix.
@@ -5307,6 +5309,17 @@ module subroutine sort_eigen_dbl(vals, vecs, ascend, err)
53075309
class(errors), intent(inout), optional, target :: err
53085310
end subroutine
53095311

5312+
module subroutine sort_int32_array(x, ascend)
5313+
integer(int32), intent(inout), dimension(:) :: x
5314+
logical, intent(in), optional :: ascend
5315+
end subroutine
5316+
5317+
module subroutine sort_int32_array_ind(x, ind, ascend, err)
5318+
integer(int32), intent(inout), dimension(:) :: x
5319+
integer(int32), intent(inout), dimension(:) :: ind
5320+
logical, intent(in), optional :: ascend
5321+
class(errors), intent(inout), optional, target :: err
5322+
end subroutine
53105323
end interface
53115324

53125325
! ******************************************************************************

src/linalg_sorting.f90

Lines changed: 299 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -275,6 +275,73 @@ module subroutine sort_eigen_dbl(vals, vecs, ascend, err)
275275
! Formatting
276276
100 format(A, I0, A, I0, A, I0, A, I0, A)
277277
end subroutine
278+
279+
! ------------------------------------------------------------------------------
280+
module subroutine sort_int32_array(x, ascend)
281+
! Arguments
282+
integer(int32), intent(inout), dimension(:) :: x
283+
logical, intent(in), optional :: ascend
284+
285+
! Local Variables
286+
logical :: dir
287+
288+
! Initialization
289+
if (present(ascend)) then
290+
dir = ascend
291+
else
292+
dir = .true.
293+
end if
294+
295+
! Process
296+
call qsort_int32(dir, x)
297+
end subroutine
298+
299+
! ------------------------------------------------------------------------------
300+
module subroutine sort_int32_array_ind(x, ind, ascend, err)
301+
! Arguments
302+
integer(int32), intent(inout), dimension(:) :: x
303+
integer(int32), intent(inout), dimension(:) :: ind
304+
logical, intent(in), optional :: ascend
305+
class(errors), intent(inout), optional, target :: err
306+
307+
! Local Variables
308+
class(errors), pointer :: errmgr
309+
type(errors), target :: deferr
310+
character(len = :), allocatable :: errmsg
311+
integer(int32) :: n
312+
logical :: dir
313+
314+
! Initialization
315+
n = size(x)
316+
if (present(err)) then
317+
errmgr => err
318+
else
319+
errmgr => deferr
320+
end if
321+
if (present(ascend)) then
322+
dir = ascend
323+
else
324+
dir = .true. ! Ascend == true
325+
end if
326+
327+
! Input Check
328+
if (size(ind) /= n) then
329+
allocate(character(len = 256) :: errmsg)
330+
write(errmsg, 100) &
331+
"Expected the tracking array to be of size ", n, &
332+
", but found an array of size ", size(ind), "."
333+
call errmgr%report_error("sort_int32_array_ind", trim(errmsg), &
334+
LA_ARRAY_SIZE_ERROR)
335+
return
336+
end if
337+
if (n <= 1) return
338+
339+
! Process
340+
call qsort_int32_ind(dir, x, ind)
341+
342+
! Formatting
343+
100 format(A, I0, A, I0, A)
344+
end subroutine
278345

279346
! ******************************************************************************
280347
! PRIVATE HELPER ROUTINES
@@ -655,5 +722,237 @@ subroutine cmplx_partition_ind(ascend, x, ind, marker)
655722
end if
656723
end subroutine
657724

725+
! ------------------------------------------------------------------------------
726+
!> @brief A recursive quick sort algorithm.
727+
!!
728+
!! @param[in] ascend Set to true to sort in ascending order; else, false
729+
!! to sort in descending order.
730+
!! @param[in,out] x On input, the array to sort. On output, the sorted
731+
!! array.
732+
!!
733+
!! @par Notes
734+
!! This implementation is a slight modification of the code presented at
735+
!! http://www.fortran.com/qsort_c.f95.
736+
recursive subroutine qsort_int32(ascend, x)
737+
! Arguments
738+
logical, intent(in) :: ascend
739+
integer(int32), intent(inout), dimension(:) :: x
740+
741+
! Local Variables
742+
integer(int32) :: iq
743+
744+
! Process
745+
if (size(x) > 1) then
746+
call int32_partition(ascend, x, iq)
747+
call qsort_int32(ascend, x(:iq-1))
748+
call qsort_int32(ascend, x(iq:))
749+
end if
750+
end subroutine
751+
752+
! ------------------------------------------------------------------------------
753+
!> @brief A routine to perform the partioning necessary for the quick sort
754+
!! algorithm.
755+
!!
756+
!! @param[in] ascend Set to true to sort in ascending order; else, false
757+
!! to sort in descending order.
758+
!! @param[in,out] x On input, the array to sort. On output, the sorted
759+
!! array.
760+
!! @param[out] marker The partioning marker.
761+
!!
762+
!! @par Notes
763+
!! This implementation is a slight modification of the code presented at
764+
!! http://www.fortran.com/qsort_c.f95
765+
subroutine int32_partition(ascend, x, marker)
766+
! Arguments
767+
logical, intent(in) :: ascend
768+
integer(int32), intent(inout), dimension(:) :: x
769+
integer(int32), intent(out) :: marker
770+
771+
! Local Variables
772+
integer(int32) :: i, j, temp, pivot
773+
774+
! Process
775+
pivot = x(1)
776+
i = 0
777+
j = size(x) + 1
778+
if (ascend) then
779+
! Ascending Sort
780+
do
781+
j = j - 1
782+
do
783+
if (x(j) <= pivot) exit
784+
j = j - 1
785+
end do
786+
i = i + 1
787+
do
788+
if (x(i) >= pivot) exit
789+
i = i + 1
790+
end do
791+
if (i < j) then
792+
! Exchage X(I) and X(J)
793+
temp = x(i)
794+
x(i) = x(j)
795+
x(j) = temp
796+
else if (i == j) then
797+
marker = i + 1
798+
return
799+
else
800+
marker = i
801+
return
802+
end if
803+
end do
804+
else
805+
! Descending Sort
806+
do
807+
j = j - 1
808+
do
809+
if (x(j) >= pivot) exit
810+
j = j - 1
811+
end do
812+
i = i + 1
813+
do
814+
if (x(i) <= pivot) exit
815+
i = i + 1
816+
end do
817+
if (i < j) then
818+
! Exchage X(I) and X(J)
819+
temp = x(i)
820+
x(i) = x(j)
821+
x(j) = temp
822+
else if (i == j) then
823+
marker = i + 1
824+
return
825+
else
826+
marker = i
827+
return
828+
end if
829+
end do
830+
end if
831+
end subroutine
832+
833+
! ------------------------------------------------------------------------------
834+
!> @brief A recursive quick sort algorithm.
835+
!!
836+
!! @param[in] ascend Set to true to sort in ascending order; else, false
837+
!! to sort in descending order.
838+
!! @param[in,out] x On input, the array to sort. On output, the sorted
839+
!! array.
840+
!! @param[in,out] ind On input, a tracking array of the same length as @p x.
841+
!! On output, the same array, but shuffled to match the sorting order of
842+
!! @p x.
843+
!!
844+
!! @par Notes
845+
!! This implementation is a slight modification of the code presented at
846+
!! http://www.fortran.com/qsort_c.f95.
847+
recursive subroutine qsort_int32_ind(ascend, x, ind)
848+
! Arguments
849+
logical, intent(in) :: ascend
850+
integer(int32), intent(inout), dimension(:) :: x
851+
integer(int32), intent(inout), dimension(:) :: ind
852+
853+
! Local Variables
854+
integer(int32) :: iq
855+
856+
! Process
857+
if (size(x) > 1) then
858+
call int32_partition_ind(ascend, x, ind, iq)
859+
call qsort_int32_ind(ascend, x(:iq-1), ind(:iq-1))
860+
call qsort_int32_ind(ascend, x(iq:), ind(iq:))
861+
end if
862+
end subroutine
863+
864+
! ------------------------------------------------------------------------------
865+
!> @brief A routine to perform the partioning necessary for the quick sort
866+
!! algorithm.
867+
!!
868+
!! @param[in] ascend Set to true to sort in ascending order; else, false
869+
!! to sort in descending order.
870+
!! @param[in,out] x On input, the array to sort. On output, the sorted
871+
!! array.
872+
!! @param[in,out] ind On input, a tracking array of the same length as @p x.
873+
!! On output, the same array, but shuffled to match the sorting order of
874+
!! @p x.
875+
!! @param[out] marker The partioning marker.
876+
!!
877+
!! @par Notes
878+
!! This implementation is a slight modification of the code presented at
879+
!! http://www.fortran.com/qsort_c.f95
880+
subroutine int32_partition_ind(ascend, x, ind, marker)
881+
! Arguments
882+
logical, intent(in) :: ascend
883+
integer(int32), intent(inout), dimension(:) :: x
884+
integer(int32), intent(inout), dimension(:) :: ind
885+
integer(int32), intent(out) :: marker
886+
887+
! Local Variables
888+
integer(int32) :: i, j, itemp, temp, pivot
889+
890+
! Process
891+
pivot = x(1)
892+
i = 0
893+
j = size(x) + 1
894+
if (ascend) then
895+
! Ascending Sort
896+
do
897+
j = j - 1
898+
do
899+
if (x(j) <= pivot) exit
900+
j = j - 1
901+
end do
902+
i = i + 1
903+
do
904+
if (x(i) >= pivot) exit
905+
i = i + 1
906+
end do
907+
if (i < j) then
908+
! Exchage X(I) and X(J)
909+
temp = x(i)
910+
x(i) = x(j)
911+
x(j) = temp
912+
913+
itemp = ind(i)
914+
ind(i) = ind(j)
915+
ind(j) = itemp
916+
else if (i == j) then
917+
marker = i + 1
918+
return
919+
else
920+
marker = i
921+
return
922+
end if
923+
end do
924+
else
925+
! Descending Sort
926+
do
927+
j = j - 1
928+
do
929+
if (x(j) >= pivot) exit
930+
j = j - 1
931+
end do
932+
i = i + 1
933+
do
934+
if (x(i) <= pivot) exit
935+
i = i + 1
936+
end do
937+
if (i < j) then
938+
! Exchage X(I) and X(J)
939+
temp = x(i)
940+
x(i) = x(j)
941+
x(j) = temp
942+
943+
itemp = ind(i)
944+
ind(i) = ind(j)
945+
ind(j) = itemp
946+
else if (i == j) then
947+
marker = i + 1
948+
return
949+
else
950+
marker = i
951+
return
952+
end if
953+
end do
954+
end if
955+
end subroutine
956+
658957
! ------------------------------------------------------------------------------
659958
end submodule

tests/linalg_test.f90

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -354,6 +354,10 @@ program main
354354
rst = test_pgmres_1()
355355
if (.not.rst) flag = 106
356356

357+
! More Tests
358+
rst = test_int32_ascend_sort()
359+
if (.not.rst) flag = 107
360+
357361
! End
358362
if (flag /= 0) stop flag
359363
end program

0 commit comments

Comments
 (0)