再帰を使わないquicksortその２

前回のものに比べて、必要なワーク配列を減らしました。ここでは最悪の場合に必要になるサイズ、すなわち並べ替えるデータと同じ大きさの整数配列二つを用意しています。しかし、実験的に調べてみると。ランダムなデータの場合、並べ替える数が十分に増えると２０％以下しか実際は使われません。
（最悪の場合とは、ソートの途中で[2,1,4,3,6,5,....N,N-1]のような並びが現れる状況だと思います。）
クイックソートの分割が十分に進んだ段階で、挿入ソートなり別の方法を使うようにすれば、必要なワーク配列を2N/nにできると思われます。ここでNは並べ替えるデータの総数、nはソート方法を切り替えるデータ数の閾値。

再帰版に比べて１０倍くらい速いようです。90で書いていますが77にもすぐ直せる命令しか使っていません。

ソース・プログラム

module m_quick
  implicit none
contains
  subroutine qsort(x)
    real, intent(in out) :: x(:)
    integer, allocatable :: indx(:), iwk(:) 
    integer :: n
    allocate( indx(size(x)) )
    allocate( iwk (size(x)) ) 
    n = size(x)
    call qsort2(n, x, indx, iwk)
    return
  end subroutine qsort
  
  subroutine qsort2(n, x, indx, iwk)
    integer, intent(in) :: n
    real, intent(in out) :: x(n)
    integer, intent(out) :: indx(n), iwk(n)
    integer :: i, j, kk0, kk1, k0, k1, ix0, ix1, iindx, nindx
    real :: tmp
    if (n < 2) return
    indx(1) = 1
    indx(2) = n !size(x)  
    nindx = 2
    do
      if (nindx == 0) exit
      iindx = 0
      do j = 1, nindx - 1, 2
        kk0 = indx(j) 
        kk1 = indx(j + 1) 
        k0 = kk0 + 1
        k1 = kk1
        do
          if (k0 > k1) exit
          do ix0 = k0, k1              ! L
            if (x(ix0) > x(kk0)) exit 
          end do
          do ix1 = k1, k0, -1          ! R
            if (x(ix1) < x(kk0)) exit 
          end do
          if (ix0 >= ix1) exit
          tmp = x(ix0)                 !swap
          x(ix0) = x(ix1)
          x(ix1) = tmp
          k0 = ix0 + 1
          k1 = ix1 - 1
        end do
        tmp = x(kk0)                   ! CSHIFT(x(kk0:ix1), 1)
        do i = kk0, ix1 - 1
          x(i) = x(i + 1)
        end do
        x(ix1) = tmp
        if (kk0 < ix1 - 1) then         ! L
          iindx = iindx + 1
          iwk(iindx) = kk0 
          iindx = iindx + 1
          iwk(iindx) = ix1 - 1
        end if
        if (kk1 > ix0) then             ! R
          iindx = iindx + 1
          iwk(iindx) = ix0
          iindx = iindx + 1
          iwk(iindx) = kk1
        end if 
      end do
      nindx = iindx
      indx(1:nindx) = iwk(1:nindx)
    end do
    return
  end subroutine qsort2
end module m_quick

module m_kahan
  implicit none
contains
  pure function sum_kahan(xx) result(s)
    real, intent(in) :: xx(:)
    real :: c, s, x, t
    integer :: i
    c = 0.0
    s = 0.0
    do i = 1, size(xx)
      x = xx(i) - c
      t = s + x
      c = (t - s) - x
      s = t
    end do
    return
  end function sum_kahan
end module m_kahan
 
program main
  use m_quick
  use m_kahan
  implicit none
  integer, parameter :: nmax = 10**7
  real :: x(nmax), s0, s1
   REAL(8) :: t0, t1
  integer :: i
  call RANDOM_SEED()
  do i = 1, 1000
    call RANDOM_NUMBER(x)
    s0 = sum_kahan(x)
    call CPU_TIME(t0)
    call qsort(x)
    call CPU_TIME(t1)
    s1 = sum_kahan(x)
    print *, t1 - t0, ANY( x(:size(x) - 1) > x(2:) ), s0 - s1, s0
  end do
  stop
end program main