module sort
!---------------------------------------------
! SORT subroutines
!   Bubble/quick sort for integer/double real
!---------------------------------------------
! HISTORY
!  2020/10/31 [v1.0  ] <new> author = Koji KOBAYASHI
!  2020/11/ 4 [v1.0.1] <comment> usage
!---------------------------------------------
! REQUIREMENTS
!  swaps.f90 (v1.0)
!---------------------------------------------
! USAGE
!  use sort
!  integer :: ivec(5)
!  real(kind(0.0d0)) :: dvec(5)
!    ivec = (/ 5,1,2,4,3 /)
!    call sort__bubble(ivec)  !-- Bubble sort
!    print '(5i3)', ivec
!
!    dvec = (/ 5d0,1d0,2d0,4d0,3d0 /)
!    call sort__quick(dvec)  !-- Quick sort
!    print '(5f4.1)', dvec
!---------------------------------------------
! ALGORITHMS
!  Bubble sort : O(n^2)
!    suitable for a small 'n'
!  Quick sort : O(n log n) - O(n^2)
!    suitable for a large 'n' (unstable sort)
!---------------------------------------------
 use swaps
 implicit none
 private

 INTERFACE sort__bubble
    module procedure bubbleSort_i,bubbleSort_d
 end INTERFACE
 INTERFACE sort__quick
    module procedure quickSort_i,quickSort_d
 end INTERFACE

 public :: sort__bubble, sort__quick

CONTAINS

!---------------------------------------------
! Bubble sort: integer
!---------------------------------------------
pure subroutine bubbleSort_i(vec)
 integer, intent(inout) :: vec(:)
 integer :: i,j,n
   n = size(vec)
   do i=1,n-1
      do j=2,n-i+1
         if( vec(j-1) > vec(j) ) then
            call swap(vec(j-1),vec(j))
         end if
      end do
   end do
end subroutine
!---------------------------------------------
! Bubble sort: double
!---------------------------------------------
pure subroutine bubbleSort_d(vec)
 real(kind(0d0)), intent(inout) :: vec(:)
 integer :: i,j,n
   n = size(vec)
   do i=1,n-1
      do j=2,n-i+1
         if( vec(j-1) > vec(j) ) then
            call swap(vec(j-1),vec(j))
         end if
      end do
   end do
end subroutine
!---------------------------------------------



!---------------------------------------------
! quick sort: integer 
!---------------------------------------------
recursive pure subroutine quickSort_i(a,left,right)
 integer, intent(inout) :: a(:)
 integer, intent(in), optional :: left, right
 integer :: pivot, i, j, imin, jmax
   if(present(left)) then;  imin = left
   else;                    imin = lbound(a,1); end if
   if(present(right)) then; jmax = right
   else;                    jmax = ubound(a,1); end if
   if (imin < jmax) then
      i = imin; j = jmax
      pivot = this__median3_i(a(i), a(i+(j-i)/2), a(j))
      do
         do while( a(i) < pivot )
            i = i + 1
         end do
         do while( a(j) > pivot )
            j = j - 1
         end do
         if ( i >= j )  EXIT
         call swap(a(i),a(j))
         i = i + 1
         j = j - 1
      end do
      call quickSort_i(a(:),imin,i-1) !-- sort smaller part
      call quickSort_i(a(:),j+1,jmax) !-- sort larger part
   end if
end subroutine
pure function this__median3_i(a,b,c) result(med)
 integer, intent(in) :: a,b,c
 integer :: med
   if(a < b) then
      if(c < a) then
         med = a; RETURN
      else if(c < b) then
         med = b; RETURN
      end if
   else if(a < c) then
      med = a; RETURN
   else if(c < b) then
      med = b; RETURN
   end if
   med = c
end function
!---------------------------------------------
! quick sort: double 
!---------------------------------------------
recursive pure subroutine quickSort_d(a,left,right)
 real(kind(0d0)), intent(inout) :: a(:)
 integer, intent(in), optional :: left, right
 real(kind(0d0)) :: pivot
 integer :: i, j, imin, jmax
   if(present(left)) then;  imin = left
   else;                    imin = lbound(a,1); end if
   if(present(right)) then; jmax = right
   else;                    jmax = ubound(a,1); end if
   if (imin < jmax) then
      i = imin; j = jmax
      pivot = this__median3_d(a(i), a(i+(j-i)/2), a(j))
      do
         do while( a(i) < pivot )
            i = i + 1
         end do
         do while( a(j) > pivot )
            j = j - 1
         end do
         if ( i >= j )  EXIT
         call swap(a(i),a(j))
         i = i + 1
         j = j - 1
      end do
      call quickSort_d(a(:),imin,i-1) !-- sort smaller part
      call quickSort_d(a(:),j+1,jmax) !-- sort larger part
   end if
end subroutine
pure function this__median3_d(a,b,c) result(med)
 real(kind(0d0)), intent(in) :: a,b,c
 real(kind(0d0)) :: med
   if(a < b) then
      if(c < a) then
         med = a; RETURN
      else if(c < b) then
         med = b; RETURN
      end if
   else if(a < c) then
      med = a; RETURN
   else if(c < b) then
      med = b; RETURN
   end if
   med = c
end function
!---------------------------------------------

end module