この文書の URL は http://www.cc.kyoto-su.ac.jp/~mtkg/lecture/comp_B/2012/07.html です。
副プログラム（サブルーチン）
大きなプログラムは機能的にまとまったいくつかの「部品」に分割すると作りやすくなります。 このような部品のことを副プログラム (subprogram) といいます。 Fortran の副プログラムには「関数」と「手続き」の二種類があります。
関数
これまでにもいくつか組込み関数を利用してきましたが、関数は自分でも作成することができます。
簡単な例として、2つの整数を受け取り、和の値を返す関数 add を作ってみましょう。 以下のプログラムは関数 add の定義と、add を利用するメインプログラムの2つの部分から成っています。
!------------------------------------------------------------ ! 和を返す関数 !------------------------------------------------------------ integer function add(p, q) ! 関数の型 function 関数名(引数リスト) implicit none ! おまじない integer :: p, q ! 仮引数の型 add = p + q ! 関数名に値をセット return ! 呼び出し元に値を返して終了 end function add ! end function 関数名 !------------------------------------------------------------ ! メインプログラム !------------------------------------------------------------ program work0701 implicit none integer :: nx, ny integer :: add write(*,*) 'Input two integers:' read(*,*) nx, ny write(*,*) add(nx, ny) stop end program work0701
関数の定義
関数 add の構造を詳しくみてみましょう。
■ integer function add(p, q)
関数は function 文で定義します。 function 文は次のように書きます。
関数の型 function 関数名(a1, a2, ..., an)
「関数の型」には関数が返すデータ型を指定します。
関数 add は整数を返すので関数の型は integer となります。
a1, a2, ..., an は関数が受け取る「仮引数」 (dummy argument) と呼ばれます。
仮引数が複数個あるときはカンマ , で区切って並べます。
引数を持たない関数を定義することもできます。
関数 add では p, q が仮引数になります。
メインプログラム中では add(nx, ny) と呼び出されているので、実行時に
p, q は nx, ny の値を持つことになります。
メインプログラム側の引数のことを「実引数」 (actual argument) といいます。
■ integer :: p, q
仮引数の型を宣言します。
■ add = p + q
関数名に値をセットすることで、関数の戻り値が設定されます。 この代入は関数の中で何度行っても構いませんが、最後に代入された値が関数の値となります。
■ return
関数を終了し、関数の呼び出しもと（この場合はメインプログラム）に制御を戻します。 return 文が end 文の直前にあるときは省略することができます。
■ end function add
「 end function 関数名 」で関数定義の終わりを示します。
メインプログラムの注意点
関数 add をメインプログラムや他の関数から利用する場合は、関数 add が整数型の関数であることをコンパイラに教えてあげなくてはいけません。 宣言方法は変数宣言の場合とまったく同様です。
integer :: add
組込み関数の場合はコンパイラが関数の型を知っているので宣言する必要はありません。
練習問題 (1) — 2数の最大値
2つの整数を受け取り、大きい方の値を返す関数 max2 を作れ。
練習問題 (2) — 3数の最大値
3つの整数を受け取り、最大値を返す関数 max3 を作成せよ。
階乗
0 以上の整数 n を受け取り、階乗 n! の値を返す関数 fact を考えてみましょう。
整数型の範囲では狭いので、関数の返す値は real(8) 型とします。
!------------------------------------------------------------ ! 階乗 n! !------------------------------------------------------------ real(8) function fact(n) implicit none integer, intent(IN) :: n ! 仮引数 ! local integer :: i ! ローカル変数（関数の中だけで使う作業変数） fact = 1.0d+0 do i = n, 1, -1 fact = fact * i end do return end function fact !------------------------------------------------------------ ! メインプログラム !------------------------------------------------------------ program work0704 implicit none real(8) :: fact integer :: n do n = 0, 20, 2 write(*,'(I5,E20.12)') n, fact(n) end do stop end program work0704
仮引数の有効範囲
ある関数で宣言された変数は、基本的にその関数の中だけで有効です。 他の関数やメインプログラムに同じ名前の変数があっても、まったく別のものとして扱われます。
この例では n という変数がメインプログラムと関数 fact の両方に現れていますが、関数の仮引数は関数の中だけで有効で、メインプログラムの n とは何の関係もありません。 この性質があるため、関数の設計はメインプログラムや他の関数とまったく独立に行うことができます。
ただし、Fortran には「仮引数の値を変更すると呼び出し元の実引数の値も変更される」という決まりがあります。 不用意に仮引数の値を変更しないように注意が必要です。 （以下の「仮引数の値の変更」を参照してください。）
手続き
値を返さない副プログラムを「手続き」（サブルーチン）といいます。 手続きを使うといくつかの処理をひとまとめにすることができ、プログラムの見通しがよくなります。
例として、呼び出された回数を画面に出力する手続き count を作ってみましょう。
!------------------------------------------------------------ ! 呼び出された回数を表示する !------------------------------------------------------------ subroutine count() ! subroutine 手続き名(引数リスト) implicit none integer, save :: n = 1 ! ローカル変数の宣言 write(*,*) n n = n + 1 return ! 呼び出し元に戻る end subroutine count ! end subroutine 手続き名 !------------------------------------------------------------ ! メインプログラム !------------------------------------------------------------ program work0705 implicit none call count() call count() call count() stop end program work0705
このプログラムを実行すると次のようになります。
% ./a.out
1
2
3
手続きの定義
手続きは subroutine 文で定義します。 手続きの構造は関数とほとんど同じですが、値を返さないので subroutine の型というものはありません。
subroutine 手続き名(a1, a2, ..., an) implicit none : return end subroutine 手続き名
手続きの呼び出し
手続き副プログラムは call 文によって呼び出します。
call 手続き名(引数リスト)
引数がない場合は、引数リストを省略し、次のように呼び出すこともできます。
call 手続き名() call 手続き名
save 属性
副プログラム（関数と手続き）の中で宣言された変数や配列の内容は、値を設定する前は「未定義」です。 また、副プログラムの中でいったん値を設定しても、return 文または end 文によって副プログラムを終了した時点で再び未定義となり、次にまたその副プログラムが呼ばれても前回の値は残っていません。 （残っていることもありますが、それはたまたまです。）
ただし、型宣言のときに save 属性を付ければ、変数の値を副プログラムの実行後も保持しておくことができます。
サブルーチン count では変数 n に save 属性が付けられており、初期値 1 が設定されています。
したがって、count がはじめて呼ばれたとき、
write(*,*) n の実行時の n は 1 になり、
次の行の n = n + 1 によって値が 2 に更新されます。
2回目に呼ばれたときは write(*,*) n の実行時の n は 2 であり、
次の行の n = n + 1 によって値が 3 に更新されることになります。
仮引数の値の変更
Fortran には「副プログラムの中で仮引数の値を変更すると呼び出し元の実引数の値も変更される」という重要な規則があります（これを引数の「参照渡し」といいます）。 次のプログラムでは手続き twice の中で仮引数の値が2倍されます。 結果的に呼び出し元の実引数の値も2倍されることを確認して下さい。
!------------------------------------------------------------ ! 仮引数を2倍する !------------------------------------------------------------ subroutine twice(n) implicit none integer, intent(INOUT) :: n n = n * 2 return end subroutine twice !------------------------------------------------------------ ! メインプログラム !------------------------------------------------------------ program work0706 implicit none integer :: n = 10 write(*,*) 'Before: ', n call twice(n) write(*,*) 'After : ', n stop end program work0706
このプログラムを実行すると次のようになります。
% ./a.out Before: 10 After : 20
手続き twice の実行の前後で、メインプログラム中の変数 n の値が 10 から 20 に変化したことに注意してください。
注意
twice の引数に定数を与えるとエラーになります。
次の例では twice の中で 10 = 10 * 2 という無効な代入文が実行されることになります。
call twice(10) ! Bus error
引数の変更は慎重に行いましょう。
引数を通じて値を返す
「副プログラムの中で仮引数の値を変更すると呼び出し元の実引数の値も変更される」という特徴を利用すると、手続きと呼び出し元の間でデータのやり取りをすることができます。 関数は呼び出し元にひとつの値しか返せませんが、手続きの引数を利用すると複数の値を同時に返すことができます。 次のプログラムは与えられた配列 x の最大値と最小値を同時に返す手続き副プログラムです。
!------------------------------------------------------------ ! 配列 x(1:n) の最小値 xmin と最大値 xmax を返す !------------------------------------------------------------ subroutine minmax(n, x, xmin, xmax) implicit none integer, intent(IN) :: n ! 配列 x の寸法 real(8), intent(IN) :: x(n) ! 配列 x real(8), intent(OUT) :: xmin, xmax ! x の最小値と最大値 xmin = minval(x(1:n)) xmax = maxval(x(1:n)) return end subroutine minmax !------------------------------------------------------------ ! メインプログラム !------------------------------------------------------------ program work0707 implicit none integer, parameter :: nmax = 5 integer :: i real(8) :: vec(nmax) = (/ (i, i = 1, nmax) /) real(8) :: vmin, vmax write(*,'(5F12.7)') vec(1:nmax) call minmax(nmax, vec, vmin, vmax) write(*,'(2F12.7)') vmin, vmax stop end program work0707
配列の渡し方
配列を副プログラムの仮引数にする場合、配列の寸法も仮引数として副プログラムに渡すことによって、いろいろなサイズの配列に対応することができます。 これを「整合配列」といいます。
仮引数の intent 属性
副プログラムの中で仮引数の値を設定（変更）すると、呼び出し元の実引数の値も変更されます。 このような動作を「副作用がある」ということがあります。 意図しない値の変更（うっかりミス）を防ぐため、Fortran 90 以降では副プログラム中での仮引数の取り扱いを intent 属性によって指定することができます。
intent 属性には次の3種類があります。
intent(IN) ..... 仮引数が副プログラム中で参照されるだけであり、値は変更されない。
intent(OUT) .... 仮引数を参照する前に副プログラム中で必ず値を設定しなければならない。
仮引数によって呼び出し元から情報を受け取るのではなく、副プログラム
から呼び出し元に情報を受け渡す場合に用いられる。
intent(INOUT) .. 副プログラム中での仮引数の参照・値の設定に特に制限を設けない。
intent 属性の省略
intent 属性は省略することができます。 その場合は intent(INOUT) が指定されたものとみなされます。
自動配列
Fortran 90 以降では副プログラムのローカル変数（その副プログラム内でのみ参照可能な変数）として「自動配列」 (automatic array) を使うことができます。 自動配列の寸法は副プログラムに引数として与えられた整数などを使って指定します。
次の例は与えられた配列を交換する手続きです。
subroutine swap(n, a, b) implicit none integer, intent(IN) :: n real(8), intent(INOUT) :: a(1:n), b(1:n) ! local real(8) :: work(1:n) ! 自動配列の宣言 work = a a = b b = work return end subroutine swap
自動配列は副プログラムの実行開始時に記憶領域が確保され、終了時に解放されます。 したがって、自動配列に save 属性をつけることはできません。
割り付け配列 (allocatable array)
自動配列の代わりに「割り付け配列」 (allocatable array) を使うこともできます。 割り付け配列は変数宣言時に配列の階数 (rank) だけを指定しておき、必要な大きさのメモリを実行時に割り付けることができます。
subroutine swap(n, a, b) implicit none integer, intent(IN) :: n real(8), intent(INOUT) :: a(1:n), b(1:n) ! local real(8), allocatable :: work(:) ! 割り付け配列の宣言 allocate(work(1:n)) ! サイズ指定（メモリ確保） work = a a = b b = work deallocate(work) ! メモリ解放 return end subroutine swap
自動配列 or 割り付け配列？
割り付け配列には allocate, deallocate が必要になる点を除くと、自動配列と割り付け配列の機能はほぼ同じです。 副プログラムの中ではどちらを使うべきでしょうか？
自動配列はメモリの「スタック」と呼ばれる領域を使用する場合が多いのに対し、割り付け配列はメモリの「ヒープ」領域を使用する場合が多いとされています。 このため、割り付け配列よりも自動配列のほうがメモリ割り付けが高速であると期待されます。 しかしながら、使用可能なスタックのサイズには制限がある場合が多く、自動配列で使用できる配列サイズが制限されることがあります。
したがって、まず自動配列を試してみて、スタック不足が起こるようなら割り付け配列に切り替えるのがよさそうです。
ソースファイルの分割
機能的に独立した副プログラム（関数やサブルーチン）ごとにソースファイルを分けておくと、プログラムの再利用が簡単になります。
以下の例では、階乗を計算する関数 fact() だけを含むソースファイル fact.f90 と、それを利用する主プログラムのソースファイル main.f90 を別々に作成し、この2つのソースファイルから実行可能ファイル a.out を作っています。 fact.f90 や main.f90 といったファイル名は好きに付けることができますが、ファイルに含まれる副プログラムや主プログラムの名前と同じにしておくと分かりやすいでしょう。
fact.f90 の内容
real(8) function fact(n) implicit none integer, intent(IN) :: n integer :: i fact = 1 do i = 2, n fact = fact * i end do return end function fact
main.f90 の内容
program main implicit none real(8) :: fact write(*,*) '20! = ', fact(20) stop end program main
コンパイル方法
複数のソースファイルをコンパイルするにはいくつか方法がありますが、一番簡単なのは必要なソースファイルをすべてコンパイラに渡してしまうことです。 こうすればあとは勝手にコンパイラが処理してくれます。
% gfortran fact.f90 main.f90 % ./a.out 20! = 2.4329020081766400E+018
本日の課題
1. 1 から n までの和
整数 n を引数として受け取り、1 から n までの和を返す関数 sum1 を作成せよ。 和の計算には do ループを利用し、 n = 1, 10, 100 に対して検算を行うメインプログラムも作成すること。
2. メッセージの表示
画面に hello と表示する手続き型サブルーチン print_hello を作成せよ。
3. ルジャンドル多項式
n 次のルジャンドル多項式 Pn(x) を返す関数 lp(n, x) を作成せよ。
[-1, 1] の実数 x に対して Pn(x) は次の漸化式で定義される。
P0(x) = 1 P1(x) = x Pn(x) = [ (2*n-1) * x * Pn-1(x) - (n-1) * Pn-2(x) ] / n (n = 2, 3, ...)
解答例のプログラムを実行すると x, P1(x), P2(x), P3(x), P4(x) の値が表になって出力される。 適当なデータファイルに結果を保存し、Gnuplot などでグラフを描いてみよ。 P2(x) のグラフをプロットするには次のようにする。
% ./a.out > lp.dat % gnuplot gnuplot> plot 'lp.dat' using 1:3 w l
関数 lp() だけを含むソースファイル lp.f90 と、lp() を利用してデータを出力する主プログラム lp_main.f90 に分けておくと、lp() を再利用するとき便利である。
チャレンジ問題
1. Zeller の公式
Zeller の公式を用いると年月日から曜日を調べることができる。 y を西暦年、m を月、d を日とすると、
w = y + [y/4] - [y/100] + [y/400] + [2.6*m + 1.6] + d
で定義される w を 7 で割ったあまりが曜日となる（0 から 6 が日曜日から土曜日に対応する）。 ただし、1月と2月は前年の13月と14月とすること （例えば、2001年1月20日の曜日を計算する場合、2000年13月20日として上式に代入する）。 [x] はガウス記号で x を超えない整数を表す。
年月日を表す3つの整数 year, month, day を引数とし、曜日に対応する整数を返す関数
zeller(year, month, day) を定義せよ。
1月と2月の場合に、不用意に仮引数を変更しないように気をつけること。
動作検証用に今日の曜日を調べるメインプログラムも作成せよ。
integer function zeller(year, month, day) implicit none integer, intent(IN) :: year, month, day : return end function zeller
ヒント： 1月と2月の場合の修正を加えた年月日を y, m, d とすると、整数 w は以下で計算できる。
w = y + y/4 - y/100 + y/400 + (13*m+8)/5 + d
2. 13日の金曜日
西暦2000年から2999年の1000年間に13日の金曜日は何回あるか。
3. エルミート多項式
n 次のエルミート多項式 Hn(x) を返す関数を作成せよ。 任意の実数 x に対して Hn(x) は次の漸化式で定義される。
H0(x) = 1 H1(x) = 2 * x Hn(x) = 2 * x * Hn-1(x) - 2 * (n-1) * Hn-2(x) (n = 2, 3, 4, ...)
Gnuplot を用いて [-10, 10] の範囲で関数 Hn(x) * exp(-x*x/2) をプロットせよ。
