Fortran で関数・サブルーチンを引数とする方法

# Fortran で副プログラムを引数とする

Fortran で関数やサブルーチンという副プログラムを引数とすることは、FORTRAN II の頃から出来るのですが（1961年の改訂マニュアルに記述あり）、あまり知られていないようなので、台形積分のサブルーチンに被積分関数を引数として渡す例を示しておきます。

特に Fortran2008 になって、内部副プログラム (internal procedure) を引数として利用できるようになったので、便利度が上がっています。さらに Fortran2018 では内部副プログラムの変数スコープの制御が import 文によって可能になるので、今後更に安全に利用できるようになると思われます。

副プログラムを引数として渡す方法は、FORTRAN II と FORTRAN66/77 と Fortran90 以降で異なるので、新しい方からプログラム例を示します。

## Fortran 2008

Fortran90 以降では、引数の型に相当するものを interface として記述します。ここでは、Fortran2008 の例を示します。

Fortran90/95/2003 では内部副プログラムに制限があって、これを他の副プログラムの引数として利用することが出来ません。その理由は、内部副プログラムは Fortran90 の時に文関数を代替するものとして導入されたもので、変数スコープを親プログラムと共有するなどインライン展開できるような軽い計算に利用されることを想定していたためではないかと推測されます。

Section 10.4
https://wg5-fortran.org/N1801-N1850/N1828.pdf

### プログラム

ここでは、interface 部分をモジュールの定義部に取り出して、関数定義部では　procedure 文を用いてスッキリさせています。

:::note info
追記 コメント欄で @septcolor さんに教えてもらったのですが、interface 内では module の頭で宣言した implicit none が適用されないということで、interface に implicit none を付ける方が安心・安全な気がします。なお import では module 内の定義や use は拾えますが implicit 宣言は拾えないようです。 
:::

```fortran
    module quadrature_m
        implicit none
        
        interface
            real function f_t(x)
                implicit none
                real, intent(in) :: x
            end function f_t
        end interface 
        
    contains
        subroutine trapez(f, n, x0, x1, s)
            procedure(f_t) :: f 
            integer, intent(in) :: n
            real, intent(in ) :: x0, x1
            real, intent(out) :: s
            real :: y(0:n), h
            integer :: i
            h = (x1 - x0) / n
            do i = 0, n
                y(i) = f(x0 + i * h)
            end do
            s = h * ( y(0) / 2 + sum(y(1:n-1)) + y(n) / 2 )
        end  subroutine 
    end module quadrature_m
    
    program quadature
        use quadrature_m
        implicit none
        real :: x0, x1, s0, s1, s
        x0 = 0.0
        x1 = 4 * atan(1.0) ! pi
        call trapez(func, 32, x0, x1, s0)
        call trapez(func, 64, x0, x1, s1)
        s = (4 * s1 - s0) / 3 ! Richardson extrapolation
        
        print *, '\int_0^\pi sin(x) dx'
        print *, 'x0, x1, S', x0, x1, s
        print *, 's0, s1', s0, s1
    contains 
        real function func(x)
            real, intent(in) :: x
            func = sin(x) 
        end function func
    end program quadature
```

interface を、引数の定義の所に直接書くこともできます。
        
```fortran
        subroutine trapez(f, n, x0, x1, s)
            interface
                real function f(x)
                    real, intent(in) :: x
                end function f
            end interface 
            integer, intent(in) :: n
            real, intent(in ) :: x0, x1
            real, intent(out) :: s
            real :: y(0:n), h
            integer :: i
            h = (x1 - x0) / n
            do i = 0, n
                y(i) = f(x0 + i * h)
            end do
            s = h * ( y(0) / 2 + sum(y(1:n-1)) + y(n) / 2 )
        end  subroutine 
```

### 実行結果

$\int_0^\pi \sin(x) dx = 2$ を計算しています。台形積分を刻み幅を半分にして二回行い、これを Richardson 補外して Simpson 則相当の結果を得ています。 

```
 \int_0^\pi sin(x) dx
 x0, x1, S  0.0000000E+00   3.141593       2.000000
 s0, s1   1.998394       1.999598
```

# 付録

以下に FORTRAN IV と FORTRAN II での例を IBM7094 Emulator での実行例とともに与えておきます。

### emulator

http://simh.trailing-edge.com/

## FORTRAN IV

FORTRAN66/77 では external 文によって引数として与える副プログラム名を呼び出し側で指名します。Classical FORTRAN は独立コンパイル形式になっているので、リンク時の引数の数や型のチェックが行われません。したがって interface もありません。


### JOB CARD

```fortran
$JOB           FUNCTION ARGUMENT
$EXECUTE       IBJOB
$IBJOB         GO,LOGIC,MAP,FIOCS
$IBFTC MAIN    
C MAIN 
C
      EXTERNAL SIN, FUN
      PI = 4.0 * ATAN(1.0)
C QUADRATURE SIN(X) (0..PI)
      CALL TRPZD(SIN, 32, 0.0, PI, S0)
      CALL TRPZD(SIN, 64, 0.0, PI, S1)
C RICHARDSON EXTRAPOLATION (SIMPSON''S RULE)
      S2 = (4.0 * S1 - S0) / 3.0
      WRITE(6,100) S0, S1, S2
  100 FORMAT(1H ,11HINTEGRATION, 3F12.6)


C QUADRATURE X**4(1-X)**4/(1+X**2)=22/7-PI (0..1)
      CALL TRPZD(FUN, 32, 0.0, 1.0, S0)
      CALL TRPZD(FUN, 64, 0.0, 1.0, S1)
C RICHARDSON EXTRAPOLATION (SIMPSON''S RULE)
      S2 = (4.0 * S1 - S0) / 3.0
      A = 22.0 / 7.0 - PI
      WRITE(6,100) S0, S1, S2, A
      STOP
      END
$IBFTC FUNC    
C
      FUNCTION FUN(X)
      FUN = X**4*(1.0-X)**4/(1.0+X*X)
      RETURN
      END
C

$IBFTC TRAP    
C TRAPEZOID RULE (FOR LIBRARY FUNCTION)
C
      SUBROUTINE TRPZD(FUNCY, N, X0, X1, S)
      H = (X1 - X0) / FLOAT(N)
      S = 0.5 * (FUNCY(X0) + FUNCY(X1))
      N1 = N - 1
      DO 10 I = 1, N1
          X = X0 + FLOAT(I) * H
          S = S + FUNCY(X)
   10 CONTINUE   
      S = S * H 
      RETURN
      END  
```

syntax coloring が崩れるので Simpson's rule のところに余分なアポストロフィを足しておきました。

### 実行結果　抜粋
$\int_0^\pi \sin x dx = 2$ の他に、$\pi$ の近似値の比較 $22/7-\pi$ の計算も行っています。$$\int_0^1{x^4(1-x)^4\over1+x^2} dx = {22\over7}-\pi>0$$


```
INTEGRATION    1.998393    1.999598    2.000000
INTEGRATION    0.001264    0.001264    0.001264
INTEGRATION    0.001264
        511 LINES OUTPUT.
           $IBSYS
           $STOP
```

## FORTRAN II

FORTRAN II では、F66/77 における external 文に相当するものが、第一カラムに F の文字を書くことになっています。つまりコンパイラ指示行のような形式になっています。これは FORTRAN II の倍精度や複素数などでも用いられる形式です。

FORTRAN II の reference manual の初版には、副プログラムの引数としての利用の記述はありませんが、1961年1月発行の "reference manual 709/7090 FORTRAN Programming System" C28-6054-2 には記述があります。また裏表紙のところに改訂にあたって "704 and 709 FORTRAN: Using Function and Subroutine Names as Arguments, fonn J28 - 6082" の内容を取り込んだとあるので、これ以前にパッチのような形で、副プログラムの引数化の拡張がなされていたと思われます。

![102663112.05.01.acc.pdf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44522/79e71f74-d6ce-5111-0315-49529caf12a7.png)

https://www.computerhistory.org/collections/catalog/102663112
  　
### JOB CARD

FORTRAN II では組み込み関数はライブラリ関数として与えられていて、その名前は呼び出し時には型に合わせて末尾や先頭に F(Floating point）や X(fiXed point) の文字を付加せねばならないので、利用者定義関数と命名規則が微妙に異なっているようです。これに対応してふたつの SUBROUTINE を用意しました。

```FORTRAN
$JOB           TEST
$EXECUTE       FORTRAN
*      ID      FCARD
*      LIST8
*      XEQ
C
      FUNCTION FUN(X)
      FUN = X**4*(1.0-X)**4/(1.0+X*X)
      RETURN
      END
C
C TRAPEZOIDAL RULE (FOR LIBRARY FUNCTION)
C
      SUBROUTINE TRPZD1(FUNC, N, X0, X1, S)
      H = (X1 - X0) / FLOATF(N)
      S = 0.5 * (FUNCF(X0) + FUNCF(X1))
      N1 = N - 1
      DO 10 I = 1, N1
          X = X0 + FLOATF(I) * H
          S = S + FUNCF(X)
   10 CONTINUE   
      S = S * H 
      RETURN
      END  
C
C TRAPEZOIDAL RULE (FOR USER FUNCTION)
C
      SUBROUTINE TRPZD2(FUNC, N, X0, X1, S)
      H = (X1 - X0) / FLOATF(N)
      S = 0.5 * (FUNC(X0) + FUNC(X1))
      N1 = N - 1
      DO 10 I = 1, N1
          X = X0 + FLOATF(I) * H
          S = S + FUNC(X)
   10 CONTINUE   
      S = S * H 
      RETURN
      END  
C
C MAIN 
C
F     FUN, SIN
      PI = 4.0 * ATANF(1.0)
C QUADRATURE SIN(X) (0..PI)
      CALL TRPZD1(SIN, 32, 0.0, PI, S0)
      CALL TRPZD1(SIN, 64, 0.0, PI, S1)
C RICHARDSON EXTRAPOLATION (SIMPSON''S RULE)
      S2 = (4.0 * S1 - S0) / 3.0
      PRINT 100, S0, S1, S2
  100 FORMAT(1H ,11HTEST F CARD, 3F12.6)


C QUADRATURE X**4(1-X)**4/(1+X**2)=22/7-PI (0..1)
      CALL TRPZD2(FUN,  8, 0.0, 1.0, S0)
      CALL TRPZD2(FUN, 16, 0.0, 1.0, S1)
C RICHARDSON EXTRAPOLATION (SIMPSON''S RULE)
      S2 = (4.0 * S1 - S0) / 3.0
      A = 22.0 / 7.0 - PI
      PRINT 100, S0, S1, S2, A
      STOP
      END
```



### 実行結果

```
 TEST F CARD    1.998393    1.999598    2.000000
 TEST F CARD    0.001265    0.001264    0.001264
 TEST F CARD    0.001265
```

![fcardpng.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44522/07cfb48e-9a91-c0a1-0a09-e194b8847f1e.png)