fortran66のブログ

fortran について書きます。

<del datetime="2010-02-24T23:43:13+09:00">Shell sort実行時間がO(N^2)で謎だった件</del>Shell sort かわいいよ Shell sort!の巻


動的に記憶領域を確保できなかったFORTRAN77時代に重宝していたShell Sortのプログラムを書いてみました。ギャップのとり方として、Knuthの数列[tex:\{g_{n+1}=3*g_n+1|2*g_ha102549)。しかし、実行時間をデータ数の関数としてみると大体O(N^2)に比例していて、Shell sortの計算量の理論値O(N^1.2〜1.5)と食い違っています。

少し調べたところ、Bubble sortやKnuth数列以外のギャップでの場合と比較すると、データスワップ回数はBubble sortではO(N^2)、Shell sortはO(N^1.3)前後で理論値通りになっていました。またスワップしない場合も含めた総ループ回数は、どの場合もO(N^2)に比例していて、ファクターだけがそれぞれの場合で異なっており、それが大体実行時間の比になっておりました。

これが本当なら、「現今のIntel Fortran実行環境下ではデータのスワップは無視できる時間的寄与しかしないため、shell sortの計算時間はデータスワップ回数(O(N^1.2〜1.5))ではなく、総ループ回数(O(N^2))に比例する。」ということがいえると思います。

よく確かめていないので、この結論は一切保障しませんw 追試してください。

FORTRAN77時代には、動的に記憶領域をとらなくてよいShell sortが好きで、よく使っておりました。何か釈然としない?www

余りに釈然としないので77時代のソースと比較したところELSE節が抜けていることを見出しました。我ながらひどすぎw Shell sort かわいいよ Shell sort! 少しでも疑った僕を許しておくれww

■実行結果

Knuthのgap間隔に対するShell sortの計算量の理論値O(N^1.2)と大体整合します。10億要素までやってみました。

■原始プログラム

MODULE m_sort
   IMPLICIT NONE
 CONTAINS

   PURE FUNCTION ssort(x) ! Shell Sort
     REAL, ALLOCATABLE :: ssort(:)
     REAL, INTENT(IN) :: x(:)  
     REAL :: tmp
     INTEGER :: i, j, igap
     
     ssort = x
     igap = (3**INT( LOG(SIZE(x) + 1.0) / LOG(3.0) ) - 1) / 2 ! Knuth
     DO 
      IF (igap < 1) EXIT     
      DO i = igap, SIZE(x)
       DO j = i - igap, 1, -igap
        IF ( ssort(j) > ssort(j + igap) ) THEN
          tmp             = ssort(j)
          ssort(j)        = ssort(j + igap)
          ssort(j + igap) = tmp
        ELSE
          EXIT
        END IF
       END DO
      END DO
      igap = (igap - 1) / 3
     END DO
    
     RETURN
   END FUNCTION ssort

END MODULE m_sort
!========================================================
PROGRAM sort
 USE m_sort
 IMPLICIT NONE
 REAL, ALLOCATABLE :: x(:), y(:)
 REAL :: t0, t1
 INTEGER :: i

 CALL RANDOM_SEED()

 DO i = 0, 9
   ALLOCATE( x(10**i) )
   CALL RANDOM_NUMBER(x)
   CALL CPU_TIME(t0)
   y = ssort(x)
   CALL CPU_TIME(t1)
   PRINT *, i, 10**i, t1 - t0
   DEALLOCATE( x, y )
 END DO

 STOP
END PROGRAM sort