再帰を使わない quick sort

再帰を使わない Quicksort を書いてみたのですが、ついでなのでFORTRAN77 に書き直してみました。再帰を使わない分、作業領域として入力配列とほぼ同じ長さの作業配列が３本使っています。

再帰を使っていないので大幅にスピードアップしました。また、同じ再帰を使わない版でも、F77版は同等のF90版のものより高速です。メモリーの受け渡しの問題だと思うのですが、どこで差がついているのか・・・・・調べていません。

ソース

      SUBROUTINE qsort77(n, x, w, indx, iwk)
      REAL x(n), w(n)
      INTEGER indx(n + 2), iwk(n + 2)
      DO 10 i = 1, n + 2
       indx(i) = 0 
       iwk (i) = 0
   10 CONTINUE
      DO 11 i = 1, n
        w(i) = x(i) 
   11 CONTINUE           
      indx(1) = 0
      indx(2) = n + 1
      nindx = 2
      iwk(1) = 0
      DO 100 kdummy = 1, 99999999  
        IF (nindx .eq. n + 2) goto 999
        iindx = 1
        DO 20 j = 1, nindx - 1
          m = indx(j + 1) - indx(j)
          IF (m .eq. 1) THEN
            iindx = iindx + 1
            iwk(iindx) = indx(j + 1)
          ELSE IF (m .eq. 2) THEN
            iindx = iindx + 1
            iwk(iindx) = indx(j) + 1
            iindx = iindx + 1
            iwk(iindx) = indx(j + 1)
          ELSE 
            k0 = indx(j) + 1
            k1 = indx(j + 1) - 1
            ix0 = k0 
            ix1 = k1
            DO 25 k = k0 + 1, k1
              IF (x(k) .lt. x(k0)) THEN 
                w(ix0) = x(k)
                ix0 = ix0 + 1
              ELSE
                w(ix1) = x(k)
                ix1 = ix1 - 1
              END IF  
   25       CONTINUE
            w(ix0) = x(k0)
            iindx = iindx + 1
            iwk(iindx) = ix0
            iindx = iindx + 1
            iwk(iindx) = indx(j + 1)
          END IF
   20   CONTINUE
        DO 30 i = 1, n
          x(i) = w(i)
   30   CONTINUE             
        nindx = iindx
        DO 40 i = 1, nindx
          indx(i) = iwk(i)
   40   CONTINUE  
  100 CONTINUE
  999 RETURN
      END 

      PROGRAM main
      use m_kahan
      PARAMETER( nmax = 5 * 10**7  )
      REAL :: x(nmax), wk(nmax)
      integer :: indx(nmax + 2), iwk(nmax + 2)
      REAL(8) :: t0, t1
      call RANDOM_SEED()
      call RANDOM_NUMBER(x)
      call CPU_TIME(t0)
      call qsort77(nmax, x, wk, indx, iwk)
      call CPU_TIME(t1)
      PRINT *, t1 - t0
      print *, any( x(:size(x) - 1) > x(2:) )
      STOP
      END