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 文を用いてスッキリさせています。
追記 コメント欄で @septcolor さんに教えてもらったのですが、interface 内では module の頭で宣言した implicit none が適用されないということで、interface に implicit none を付ける方が安心・安全な気がします。なお import では module 内の定義や use は拾えますが implicit 宣言は拾えないようです。
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 を、引数の定義の所に直接書くこともできます。
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
FORTRAN IV
FORTRAN66/77 では external 文によって引数として与える副プログラム名を呼び出し側で指名します。Classical FORTRAN は独立コンパイル形式になっているので、リンク時の引数の数や型のチェックが行われません。したがって interface もありません。
JOB CARD
$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" の内容を取り込んだとあるので、これ以前にパッチのような形で、副プログラムの引数化の拡張がなされていたと思われます。
https://www.computerhistory.org/collections/catalog/102663112
JOB CARD
FORTRAN II では組み込み関数はライブラリ関数として与えられていて、その名前は呼び出し時には型に合わせて末尾や先頭に F(Floating point）や X(fiXed point) の文字を付加せねばならないので、利用者定義関数と命名規則が微妙に異なっているようです。これに対応してふたつの SUBROUTINE を用意しました。
$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
