FORTRAN77での実引数の扱い

FORTRAN77 での再帰を使ったプログラムで気になるは、引数が参照渡しになっていることです。局所変数は、みな AUTOMATIC 変数として呼び出しごとに確保されるからよいわけですが、引数が参照渡しである場合、呼び出しの深さに関わらず実引数の内容が書き換えられることになってしまいます。したがって、引数は INTENT(IN) の属性を持っていないと不都合がでることになります。

FORTRAN77 で再帰を使った例

階乗

      PROGRAM TEST
      DO 10 I = 1, 12
        PRINT *, I, IFACT(I, IFACT)
10    CONTINUE
      STOP
      END

      FUNCTION IFACT(N, MYSELF)
      EXTERNAL MYSELF
      IF (N .LE. 1) THEN
        IFACT = 1
      ELSE
        IFACT = N * MYSELF(N - 1, MYSELF)
      END IF
      RETURN
      END

ところで、FORTRAN77 以降では、実引数に数式を書くことができるわけですが、その場合に仮引数が INTENT(IN OUT) 属性で返り値があるとすると、返り値を戻す変数が存在しないことになって困難が生じます。ゆえに、数式を渡していい仮引数は INTENT(IN) 属性を持たなければならないとされています。同様なことは関数を実引数とする時も問題になります。

ところが FORTRAN77 では、INTENT 指定ができず、リンク時のチェックも弱いので、この手の問題が思いがけず顔を出します。

有名な FORTRAN 昔話に、数値べたうちの副プログラム呼び出しの常数実引数が、副プログラム内部書き換えられてしまって、それが主プログラム中の別の場所の同じ数字のべたうちと共有されていて、以降の主プログラムの結果がみなおかしくなってしまった事件があります。

例

以下に規約を破って、計算式を INTENT(IN OUT) 属性を持つ引数に入れて、サブルーチン内で書き換えを行った結果を示します。サブルーチンでは、引数から１引かれます。

変数 i を実引数としてそのまま入れた場合には、i は１引かれて戻ってきます。ところが、実引数として i*1 の数式を与えると、あたかも call by value 型の引数のように、変数 i の値は保護されています。これは多分 i*1 の結果がおかれた番地が書き変わっているものと思われます。

このように色々予想外の結果が生じるので、INTENT 属性を明示しない副プログラムでは注意が必要になります。

実行結果とソースプログラム

      program test
      i = 1
      print *, 'main: before dec', i
      call dec(i * 1)
      print *, 'main: after  dec', i
      
      print *
      
      print *, 'main: before dec', i
      call dec(i)
      print *, 'main: after  dec', i
      stop
      end
      
      subroutine dec(n)
      print *, 'dec : before n-1', n
      n = n - 1
      print *, 'dec : after  n-1', n
      end

77 での再帰 QuickSort を昔の処理系で

MS-FORTRAN Ver.5.1

オプション指定に Automatic 変数が見当たりませんが、とりあえず動きます。症状からするとコンパイラの default が 77 規格通りに Automatic になっている模様。マニュアルを見れば調べられますが、段ボールの奥底に眠っているので略。

ソース・プログラム

      PROGRAM MAIN
      PARAMETER(NMAX = 30)
      REAL X(NMAX), WK(NMAX)
      EXTERNAL QSORT
      N = NMAX
      DO 10 I = 1, N
	CALL RANDOM(X(I))
   10 CONTINUE
      CALL QSORT(N, X, WK, QSORT)
      PRINT *, (X(I), I = 1, N)
      STOP
      END

C  FORTRAN77 RECURSIVE SUBROUTINE
C  based on the idea by Andrew J. Miller http://www.esm.psu.edu/~ajm138/fortranexamples.html
      SUBROUTINE QSORT(N, X, WK, DUMSUB)
      REAL X(N), WK(N)
      EXTERNAL DUMSUB
      IF (N .LE. 1) RETURN
      K = 1
      J = N
      PIVOT = X(1)
      DO 10 I = 2, N
	IF (X(I) .LT. PIVOT) THEN
	  WK(K) = X(I)
	  K = K + 1
	ELSE
	  WK(J) = X(I)
	  J = J - 1
	END IF
   10 CONTINUE
      DO 20 I = 1, K - 1
	X(I) = WK(I)
   20 CONTINUE
      X(K) = PIVOT
      DO 30 I = K + 1, N
	X(I) = WK(I)
   30 CONTINUE
      CALL DUMSUB(K - 1, X,	   WK,	      DUMSUB)
      CALL DUMSUB(N - K, X(K + 1), WK(K + 1), DUMSUB)
      RETURN
      END

実行結果

Microsoft FORTRAN Powerstation 1.0a

Automatic 変数化オプションのタブを開いている。

ソース・プログラム

      PROGRAM MAIN
      PARAMETER(NMAX = 10**2)
      REAL X(NMAX), WK(NMAX)
      EXTERNAL QSORT
      N = NMAX        
      DO 10 I = 1, N
        CALL RANDOM(X(I))
   10 CONTINUE     
      CALL QSORT(N, X, WK, QSORT)  
      PRINT *, (X(I), I = 1, N) 
      STOP
      END

C  FORTRAN77 RECURSIVE SUBROUTINE
C  based on the idea by Andrew J. Miller http://www.esm.psu.edu/~ajm138/fortranexamples.html
      SUBROUTINE QSORT(N, X, WK, DUMSUB)
      REAL X(N), WK(N)
      EXTERNAL DUMSUB
      IF (N .LE. 1) RETURN
      K = 1
      J = N
      PIVOT = X(1)
      DO 10 I = 2, N
        IF (X(I) .LT. PIVOT) THEN
          WK(K) = X(I)
          K = K + 1
        ELSE
          WK(J) = X(I)
          J = J - 1
        END IF
10    CONTINUE
      DO 20 I = 1, K - 1
        X(I) = WK(I)
20    CONTINUE
      X(K) = PIVOT
      DO 30 I = K + 1, N 
        X(I) = WK(I)
30    CONTINUE
      CALL DUMSUB(K - 1, X,        WK,        DUMSUB)
      CALL DUMSUB(N - K, X(K + 1), WK(K + 1), DUMSUB)                
      RETURN
      END

Microsoft FORTRAN Powerstation 4.0

Automatic 変数化オプションのタブを開いている。

ソース・プログラム

      PROGRAM MAIN
      PARAMETER(NMAX = 10**3)
      REAL X(NMAX), WK(NMAX)
      EXTERNAL QSORT
C  Fortran90 : random_number     
      call random_number(X) 
      N = NMAX
      CALL QSORT(N, X, WK, QSORT)
      print *, x
      print *, any(x(1:n - 1) > x(2:))
      STOP
      END

C  FORTRAN77 RECURSIVE SUBROUTINE
C  based on the idea by Andrew J. Miller http://www.esm.psu.edu/~ajm138/fortranexamples.html
      SUBROUTINE QSORT(N, X, WK, DUMSUB)
      REAL X(N), WK(N)
      EXTERNAL DUMSUB
      IF (N .LE. 1) RETURN
      K = 1
      J = N
      PIVOT = X(1)
      DO 10 I = 2, N
        IF (X(I) .LT. PIVOT) THEN
          WK(K) = X(I)
          K = K + 1
        ELSE
          WK(J) = X(I)
          J = J - 1
        END IF
10    CONTINUE
      DO 20 I = 1, K - 1
        X(I) = WK(I)
20    CONTINUE
      X(K) = PIVOT
      DO 30 I = K + 1, N 
        X(I) = WK(I)
30    CONTINUE
      CALL DUMSUB(K - 1, X(1)    , WK(1)    , DUMSUB)
      CALL DUMSUB(N - K, X(K + 1), WK(K + 1), DUMSUB)
      RETURN
      END

FORTRAN77で再帰を使ってQUICK SORT

たまたま面白いサイトに出くわしました。本来再帰が許されていない FORTRAN77 で、シンタックスエラーを出さず、違法と合法のスレスレで再帰を実現する方法を発見した方がおられました。
http://www.esm.psu.edu/~ajm138/fortranexamples.html

自分自身を関数引数としてとるという、今まで聞いたことも見たこともない*1非常に興味深い方法で再帰を実現していて大変関心させられました。FORTRAN77 は規格の上では、副プログラムの変数は AUTOMATIC になっているので、引数の類をスタックに積むような事は自動でやってくれています。したがって、プログラムに何ら無理で不自然なところはありません。

ただ処理系によっては FORTRAN66 との互換性を保つために、デフォルトで副プログラムの変数を SAVE 扱いにしている場合もままあり、その場合はコンパイラオプションで本来の規格に忠実に AUTOMATIC 変数扱いにする必要があります。

非常に面白いアイデアで、私的には近年まれにみるヒット作で、FORTRAN77 もまだまだ奥深いなと思いました。

ついでなので、このアイデアに基づいて QuickSort を書いてみました。

実行結果

100個の実数を整列させています。DEBUG MODE ではサイズ０の配列を渡す時に配列はみだしエラーが起きますが、特に問題はありません。

ソース・プログラム

RANDOM_NUMBER() は Fortran90 のサブルーチンですが、面倒なので利用しました。それ以外はFORTRAN77 の範囲で記述したつもりです。（print *, x も 90 かｗ)

      PROGRAM MAIN
      PARAMETER(NMAX = 10**7)
      REAL X(NMAX), WK(NMAX)
      EXTERNAL QSORT
C  Fortran90 : random_number     
      call random_number(X) 
      N = NMAX
      CALL QSORT(N, X, WK, QSORT)
      print *, any(x(1:n - 1) > x(2:))
!      print *, x
      STOP
      END

C  FORTRAN77 RECURSIVE SUBROUTINE
C  based on the idea by Andrew J. Miller http://www.esm.psu.edu/~ajm138/fortranexamples.html
      SUBROUTINE QSORT(N, X, WK, DUMSUB)
      REAL X(N), WK(N)
      EXTERNAL DUMSUB
      IF (N .LE. 1) RETURN
      K = 1
      J = N
      PIVOT = X(1)
      DO 10 I = 2, N
        IF (X(I) .LT. PIVOT) THEN
          WK(K) = X(I)
          K = K + 1
        ELSE
          WK(J) = X(I)
          J = J - 1
        END IF
10    CONTINUE
      DO 20 I = 1, K - 1
        X(I) = WK(I)
20    CONTINUE
      X(K) = PIVOT
      DO 30 I = K + 1, N 
        X(I) = WK(I)
30    CONTINUE
      CALL DUMSUB(K - 1, X,        WK,        DUMSUB)
      CALL DUMSUB(N - K, X(K + 1), WK(K + 1), DUMSUB)                
      RETURN
      END

*1:でもどっかで見たような気がしたので、岩波FORTRAN辞典を見返してみたら、階乗の再帰的呼び出し例で、全く同じアイデアのプログラム例が載ってました。P.331

再帰を使わない 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