7. 関数とサブルーチン
これまでのプログラムは全ての処理が program から end program で囲まれた部分に記述されていたことと思う.これをメインプログラムと呼ぶ．これに対して，メインプログラム以外にもまとまった処理を1つのプログラム単位として記述しておくことが出来る．これをサブプログラムと呼んでいる．サブプログラムを用いることで同じ処理を何度も書かずに済むようになるため，プログラムを簡潔に記述することが出来る(従って無用な間違いも減らすことが出来る)．
なおFortranでのサブプログラムには関数( function )とサブルーチン( subroutine )の2種類が有る [1]．
関数は値を返すのに対してサブルーチンは値を返さないという違いが有るが，どちらも同じようなものである．
(実際にほとんどの言語でサブルーチンと関数の区別は存在しない．Fortranで言うところサブルーチンはC/C++では単に返値が void 型の関数でしかない.)
最初は少し取っ付きづらいかもしれないが，関数やサブルーチンを使いこなせるようになると格段にプログラムの開発が楽になるので積極的に利用しよう．
参考
sample1.f90 : 関数
sample2.f90 : サブルーチン
sample3.f90 : 変数のスコープ
sample4.f90 : intent, 配列渡し，save
sample5.f90 : optionalとキーワード引数
sample5.f90 : 再帰呼び出し
sample7.f90 : 内部手続き
sample8.f90 : 外部手続き
7.1. 概要
関数やサブルーチンなどのサブプログラムはひとまとまりの処理を実行する独立なプログラムと考えることが出来る．例えば，3次方程式 \(a x^3 + b x ^2 + c x + d = 0\) の数値解を与えられた初期値からなんらかの反復法を用いて求める処理を考えよう [2]．
ここで係数 \(a, b, c, d\) の値が何であっても行う処理は基本的に同じであろう．では2つの異なる係数の組に対して数値解を知りたい場合はどうすれば良いだろうか? 同じ処理なのだからコピー＆ペーストして係数の値だけ書き換えるというのも一つの案である．しかし，100個の異なる係数の組に対する解を知りたい場合にはどうすれば良いだろうか? 100回コピー＆ペーストするというのは得策とは言えない．このような場合には与えられた係数の組に対して処理を実行し，結果を返す 関数 が定義されていれば，それを100回呼び出すだけである． (\(\cos x\) のような数学関数を思い浮かべてもらえば良い.)
このように，何らかのまとまった処理を1つの独立したプログラム単位として定義することが出来れば圧倒的にプログラム作成が簡単になる．プログラム作成時に大事な事は人間は間違える動物であるという前提に立つことである．間違いを少なくするには問題を簡単にするしか方法は無い．サブプログラムはそれ単体で独立したプログラムであるから，仕様さえ確定してしまえばプログラマはそのサブプログラムの実装に集中すれば良く，それがどのように使われるかに気を使う必要は無いのである．また，サブプログラムが機能的に小さければ小さいほど(問題設定が単純になるので)実装が簡単で間違いが少ない [3]． 従って，規模の大きなプログラムを開発する際には，機能を出来る限り分割した小さな関数やサブルーチン群を実装し，それらを利用して最終的に所望の機能を実装する戦略の方が圧倒的にバグが入りにくいし，またあったとしても修正が容易になる．プログラムを作成する際には関数やサブルーチンにまとめられるような処理が無いか常に注意しておくべきである．
次図は関数とサブルーチンの概念図である．いずれもひとまとまりの処理を実行するものであり，何らかの入力を受け取り，出力を返すことが出来る．(もちろん入力や出力は無くても良い.) ひとたび正しく実装してしまえば，利用者(呼び出し側)はサブプログラムの内部の実装を気にする必要は無く，ブラックボックスとして用いることが出来る．(\(\cos x\) の値を内部でどのように計算するのか意識しながら用いることは無いのと同じである.) このため，利用者はその他の処理の実装に集中することが可能になるのである．
関数とサブルーチンの概念図．
サブプログラムはメインプログラムとは独立したプログラムなので，メインプログラムとは他の場所に定義して呼び出すことになる．以下の図はメインプログラムで全ての処理を実行する場合と，サブプログラムを利用する場合の処理の流れを概念的に示している．サブプログラムを呼び出すと，その定義に記述されている処理を実行し，終了すると呼び出し元に処理が返る仕組みになっている．もちろんサブプログラムが他のサブプログラムを呼び出すことも可能である．
サブプログラムを利用したプログラムの概念図．
7.2. 定義と呼び出し
まずはサブプログラムの定義の仕方と，その呼び出し方を学ぼう．
7.2.1. 定義場所
詳細は 内部手続きと外部手続き † に譲るが，基本的には関数やサブルーチンはメインプログラムの stop の後に contains を挿入し，そこから end program の間に定義すれば良い．
1program sample
2 implicit none
3
4 ! メインプログラムの処理
5
6 stop
7contains
8
9 ! ここに関数やサブルーチンを定義する
10
11end program sample
この方法で定義した関数やサブルーチンは特に準備せずにメインプログラムから呼出しが出来る．私見では，後に学ぶモジュールを使う場合(ソースファイルを複数に分割する場合)を除けば，このやり方が最も間違いを減らすことが出来る方法である．従って特に理由がない限りはこの方法で関数・サブルーチンを定義することを推奨する．
7.2.2. 関数
関数とは何らかの値を返すサブプログラムであり, 以下の様な形式で定義される．
1型 function 関数名(引数リスト)
2
3 変数や引数の型宣言
4
5 処理
6
7end function 関数名
例えば以下は倍精度実数の2乗を返す関数の定義である．
14 !
15 ! 関数の宣言(1) : 関数名と同じ名前の変数に返値を代入する形式
16 !
17 real(8) function square1(x) ! square1という名前で関数を宣言
18 implicit none ! 暗黙の型宣言の禁止
19 real(8) :: x ! 引数を宣言
20
21 square1 = x**2 ! 返値は関数名と同じ名前の変数に代入
22
23 return ! 呼び出し元に制御を戻す
24 endfunction square1
メインプログラムが program で始まり end program で終わるのと同様に，関数定義も function で始まり end function で終わる一つの独立したプログラム単位である．従って implicit none によって暗黙の型宣言を禁止し，使用する変数は明示的に宣言をする必要がある．メインプログラムとの大きな違いは，サブプログラムには引数(この場合は x)が入力として与えられることであり，引数についてもデータ型を宣言しなければならない．引数はカンマで区切って複数与えても良いし，それぞれが別のデータ型であっても構わない．
なおサブプログラムの定義時の引数のことを仮引数と呼び，呼び出すときに実際に与えられる引数のことを実引数と呼ぶ．一方，返値は 関数名と同じ名前の変数に値を代入する ことでそれが返値となる．返値のデータ型は1行目先頭の real(8) で指定される．return 文はサブプログラムからその呼び出し元へ制御を戻すことを意味しており，上のように他の処理が全て終わった後であれば省略しても構わない．処理の途中で return することも可能で，その場合にはそれ以降の処理は行われずに呼び出し元へと制御が戻される.例えば if 文で何らかの条件判定を行い，それ以降の処理を行う必要がないと判断した場合にはその時点で return によって関数から抜けることが出来る．
なお関数宣言時に function の前に返値のデータ型を指定するのでは無く，以下の様に result(変数名) のような形で指定することも出来る．
1function 関数名(引数リスト) result(変数名)
2
3 変数や引数の型宣言
4
5 処理
6
7end function 関数名
この場合は result で指定された変数に値を代入することで，それが関数の返値となる．従って，先ほどの square の定義は以下のように行うことも出来る．
26 !
27 ! 関数の宣言(2) : resultを使った形式
28 !
29 function square2(x) result(y) ! square2という名前で関数を宣言
30 implicit none ! 暗黙の型宣言の禁止
31 real(8) :: x ! 引数を宣言
32 real(8) :: y ! 返値となる変数(result)の宣言
33
34 y = x**2 ! 返値を代入
35
36 return ! 呼び出し元に制御を戻す
37 endfunction square2
関数の呼出は組込み関数と全く同じで，以下のように適宜引数を与えて呼び出せば良い．
1write(*,*) 'square1 => ', square1(2.0_8)
2write(*,*) 'square2 => ', square2(2.0_8)
ただし，引数の型は宣言したものと同一でなければならない．この場合は引数は real(8) で宣言されているので，
1write(*,*) square1(2.0), square2(2.0)
のような呼出しはコンパイルエラーとなることに注意して欲しい．
7.2.3. サブルーチン
サブルーチンは関数と良く似ているが，値を返さないという違いがある．定義は以下の様な形式となる．
1subroutine サブルーチン名(引数リスト)
2 変数や引数の型宣言
3
4 処理
5
6end subroutine サブルーチン名
例えば
10 !
11 ! サブルーチンの宣言
12 !
13 subroutine hello(name) ! helloという名前でサブルーチンを宣言
14 implicit none ! 暗黙の型宣言の禁止
15 character(len=*) :: name ! 引数を宣言 (任意長の文字列)
16
17 write(*, *) 'Hello ', name ! 内部の処理
18
19 return
20 endsubroutine hello
は引数で渡された文字列 name を Hello に続けて標準出力に表示するだけのサブルーチンである．関数と非常によく似た構造になっていることが分かるだろう．実際に，関数で必要だった返値の型指定が無いことを除くとほとんど同じである．
サブルーチンの呼び出しは関数と異なり call を用いて以下のような形でなければならない．
1call hello('Michel')
(逆に関数呼び出しに call を用いることはできない.)
7.3. 変数のスコープ
注意しなければならないのは，サブプログラムは1つの独立したプログラム単位であるので，その中で宣言する 変数は外部の変数とは完全に独立 であるという点である．例えばメインプログラムで宣言されている x という変数とサブプログラム中で宣言されている x という変数は完全に別のものとして扱われる．また当然サブプログラム中で使用している変数を外部から使用することは出来ない．例外(内部手続き を参照)はあるものの，基本的にはサブプログラムとメインプログラム及び他のサブプログラムは全く独立なものとして考えて良い [4]．サブプログラムと外部の情報のやり取りは基本的には引数と返値を通じて行う ものと理解して欲しい．
以下に挙げる例ではサブルーチン内部で変数を宣言し使用している．
9 ! sub1
10 subroutine sub1(exponent)
11 implicit none
12 real(8) :: exponent
13
14 !
15 ! 基本的に内部でのみ使う変数はここで宣言して使う
16 ! 特にループ変数(以下の例ではi)は必須
17 ! これらの変数は外部からは見えないので，メインプログラムや他のサブプログラム
18 ! で同名の変数が宣言されていてもそれらとは完全に独立であることに注意
19 !
20 ! (引数に値を代入して返す場合，メインプログラムの変数を参照する場合について
21 ! はレジュメのintent属性や内部手続きの節を参照のこと)
22 !
23 integer, parameter :: n = 10
24 integer :: i
25 real(8) :: x(n), sum
26
27 ! 代入
28 do i = 1, n
29 x(i) = i**exponent
30 enddo
31
32 ! 和を計算
33 sum = 0.0_8
34 do i = 1, n
35 sum = sum + x(i)
36 enddo
37
38 write(*, *) 'sum of array = ', sum
39
40 endsubroutine sub1
この変数は全てこのサブルーチン内部でのみ参照されるものであり，外部からはこれらの変数は参照出来ない．従って，別のサブルーチンが
42 ! sub2
43 subroutine sub2()
44 implicit none
45
46 ! このように定義されていても問題無い(sub1のnとは独立な変数である)
47 integer, parameter :: n = 100
48
49 write(*, *) 'n = ', n
50
51 endsubroutine sub2
のように定義されていても， sub1 中の n と sub2 中の n は全く別の変数である．
7.4. 引数の詳細
7.4.1. intent属性
関数は返値として値を返すことが出来るが，返値はあくまでも一つだけである．実用的には複数の値を結果として返して欲しい場合も多いが，そのような場合には引数に結果の値を代入して返すことが出来る [5]．このことからすぐ分かるように関数とサブルーチンには本質的な違いは無い．サブルーチンを使っても返したい値を引数に代入して返せば良いからだ．どちらを使うかは好みの問題であろう．(Fortranしか使わない人はあまり関数を使いたがらない傾向があるように思える．一方でC言語では関数の返値はエラーチェックに使うことが多いので，C言語から入った人は関数を好むかもしれない.)
さて，実際には引数で与えた変数の値を勝手に変更して欲しく無い場合もあるだろう．そのため，以下のようにサブプログラムの定義時に引数の入出力特性を指定することが出来る．
intent(in)入力用の変数に指定する．値は内部で参照されるのみで変更はされない． (変更しようとするとコンパイルエラーとなる．)
intent(out)出力用の変数に指定する． サブプログラム中で値が代入されることを意味する．
intent(inout)入出力のどちらにも用いる変数に指定する． 何も指定しない場合のデフォルト．
例えば以下のように引数に属性を指定することによって意図せず第1引数 a や第2引数 b の値が変更されてしまうバグを防ぐことが出来る．
29 !
30 ! <<< intent属性 >>>
31 !
32 ! * intent(in) => 入力用変数(値の変更不可)
33 ! * intent(out) => 出力用変数
34 ! * intent(inout) => 入出力
35 !
36 ! 以下は
37 !
38 ! c = a + b
39 !
40 ! のような処理を行うことを意図している．ユーザーはこの場合にaやbが変更されると
41 ! は予想しないであろう．誤ってサブルーチン内でaやbの値を変更するのを防ぐために
42 ! intent(in)を指定する．
43 !
44 subroutine add(a, b, c)
45 implicit none
46 real(8), intent(in) :: a, b ! 入力用変数(変更不可)
47 real(8), intent(out) :: c ! 出力用変数
48
49 ! 以下はコンパイルエラー
50 !a = 1.0_8
51
52 ! 出力用の変数に値を代入
53 c = a + b
54
55 endsubroutine add
C言語の経験者はC言語の関数の引数が値渡しなのに対してFortranの関数やサブルーチンでは参照渡しであることに注意して欲しい．C言語では明示的にポインタを(またはC++での参照を)渡さない限り呼び出し元の値が変更されることは無いが，Fortranではサブプログラム中で引数の値を変更すると呼び出し元の値まで変更されてしまうのである．
7.4.2. 配列渡し
配列も同様に関数やサブルーチンに引数として渡すことが可能である．以下の例の average1 では任意のサイズの配列を渡すことが出来る．(ただし次元は予め指定しておく必要がある.) 配列のサイズや形状が必要であれば size や shape などの組込み関数を使って求めることが出来る．一方で average2 では配列サイズを引数として明示的に渡している．配列の添字範囲を指定するなどの特別な事情が無い限りは average1 のような書き方の方がシンプルで良い．
57 !
58 ! <<< 形状引継ぎ配列の使い方 >>>
59 !
60 ! 引数の配列のサイズは自動的に呼出し時に与えた配列のサイズになる
61 ! サイズが必要な場合は組み込み関数sizeを用いて取得可能
62 !
63 function average1(x) result(ave)
64 implicit none
65 real(8), intent(in) :: x(:) ! サイズは自動的に決まる
66 real(8) :: ave
67
68 ave = sum(x) / size(x)
69
70 endfunction average1
71
72 !
73 ! <<< 配列サイズの引数渡し >>>
74 !
75 ! 配列のサイズを引数として明示的に受け取る
76 !
77 function average2(n, x) result(ave)
78 implicit none
79 integer, intent(in) :: n ! サイズを引数として受け取る
80 real(8), intent(in) :: x(n) ! サイズは引数として渡された整数
81 real(8) :: ave
82
83 ave = sum(x) / size(x)
84
85 endfunction average2
7.4.3. save属性
関数やサブルーチン内で save 属性付きで宣言された変数は前回の呼び出し時の値を記憶しておくことが出来る(C言語のstatic変数と同等である)．従って，例えば自分が呼び出された回数を保持することなどもできる．
save 属性付きの変数はプログラムの開始時に一度だけ宣言文で代入された値に初期化される．例えば以下のサブルーチン fibonacci ではプログラムの開始時に n = 1，f0 = 0，f1 = 0 と値が初期化されるが，呼び出しごとに値が変更され，プログラムが終了するまでその値を内部に保持し続ける．
87 !
88 ! <<< save属性 >>>
89 !
90 ! save属性付きの変数はプログラム実行中はその値を保持するので，複数回呼び出され
91 ! た場合には前回の呼出し時の値を記憶したままとなる
92 !
93 subroutine fibonacci()
94 implicit none
95 ! 以下の3つがsave属性付き (初回の呼出し時の値は宣言文で与える)
96 integer, save :: n = 1
97 integer, save :: f0 = 0
98 integer, save :: f1 = 0
99
100 integer :: f2
101
102 if(n == 1) then
103 write(*, *) 'Fibonacci number [', 0, '] = ', f0
104 f2 = 1
105 else
106 f2 = f0 + f1
107 endif
108
109 write(*, *) 'Fibonacci number [', n, '] = ', f2
110
111 ! 次回呼び出し用 (これらの値を記憶し続ける)
112 n = n + 1
113 f0 = f1
114 f1 = f2
115
116 endsubroutine fibonacci
なおFortranでは変数宣言時に同時に初期化を行うと，それを自動的に save 属性付きと扱うようである．(個人的にはこれは大変紛らわしい仕様だと思うのだが・・・)
すなわち
1integer :: n = 1
と宣言された変数には自動的に save 属性が付加されるため n = 1 に初期化されるのは一度だけである．一方で,
1integer :: n
2n = 1
では毎回 n = 1 に初期化される．混乱を防ぐために save 属性付きとしたい変数は明示的に save を指定し，それ以外の変数は宣言時の初期化は避けたほうが無難である．
7.4.4. optional属性とキーワード引数 †
引数の型宣言において optional 属性を指定した引数は，呼出し時に省略することが出来る． optional 属性付きの引数は，その引数が与えられたかどうかを検査する present という組込み関数と共に用いる．すなわち present(引数) は引数が与えられていれば真，そうでない場合には偽を返すので，if による条件分岐と組み合わせて用いれば良い．以下の例では引数 unit が与えられた場合にはその装置番号へ，与えられていない場合は標準出力へと出力を行う．
23 !
24 ! <<< optional属性 >>>
25 !
26 ! optional属性付きの引数は呼出し時に与えなくても良い．与えられなかった場合のデ
27 ! フォルトの振る舞いはユーザーの責任で実装しなければならない．
28 ! 以下では出力先を引数で与える装置番号にするか，デフォルトの標準出力にするかを
29 ! 選択することが出来る．
30 !
31 subroutine hello(name, unit)
32 implicit none
33 character(len=*), intent(in) :: name
34 integer, intent(in), optional :: unit ! optional属性付き引数
35
36 integer :: u
37
38 ! 組込み関数presentで引数が呼出し時に指定されたかどうかを調べることが出来る．
39 ! 返値が真なら指定有り，偽なら指定無し (偽の場合はデフォルトの動作を実装)
40 if(present(unit)) then
41 u = unit ! unitを指定
42 else
43 u = 6 ! デフォルトは標準出力
44 endif
45
46 write(u, *) 'Hello ', name ! 表示
47
48 return
49 endsubroutine hello
これまで関数やサブルーチンを呼び出す際には定義時の引数並びの順番通りに与えなければならなかった．しかしキーワード引数という機能を用いて，順番を気にせず引数を与えることも可能である．( open 文の使い方を思い出そう.) すなわち，上で定義された hello を呼び出す際に
1call hello(unit=0, name='Albert') ! 標準エラー出力へ
2call hello(name='Einstein') ! 標準出力へ
のように引数を"仮引数名 = 値"という形式で渡すことで，引数の順番を意識せずに呼び出しが出来る．
なお，この例では unit は optional 属性付きで宣言されているので省略することも出来る. optional 属性が無い引数については，キーワード引数を用いれば順番は気にしなくて良いが，全ての引数を指定する必要がある．
なお，このようにキーワード引数の機能を用いるには 内部手続き として宣言するか， `外部手続き`_ の場合には interface 宣言で明示的に仮引数名を呼び出し側に知らせてやらなければならない．(やはり外部手続きは面倒である.)
7.5. 再帰呼び出し(recursive) †
再帰呼び出しとは，関数やサブルーチンの中で自分自身を呼び出すことである．このような再帰手続は明示的に recursive を用いて関数やサブルーチンを定義しなければならない．なお，何も考えずに自分自身を呼び出すと簡単に無限ループになってしまうので，そうならないように注意しよう．例えば以下は階乗の計算をする例である．
10 !
11 ! <<< recursive (再帰的呼び出し) >>>
12 !
13 ! ちなみに以下の実装で正しい値が得られるのは n = 12 までである．なぜか？
14 !
15 recursive function fact(n) result(m) ! recursiveを指定する
16 implicit none
17 integer, intent(in) :: n
18 integer :: m
19
20 if(n == 1) then ! 無限ループにならないように
21 m = 1
22 else
23 m = n * fact(n - 1) ! 自分自身を(異なる引数で)呼び出す
24 endif
25
26 endfunction fact
ここで \(n ! = n \times (n-1) !\) という漸化式を用いている．ある種のアルゴリズムは再帰を使うと非常にスッキリと書くことができるので重宝することも多いだろう．ただし関数やサブルーチンの呼び出しそのものにもコスト(時間)がかかるので，不用意に用いるとパフォーマンスのボトルネックになることもあるため注意して欲しい．(基本的には再帰呼び出しを使わない方がパフォーマンスは良くなる場合が多い．)
7.6. 内部手続きと外部手続き †
定義場所 では，簡単のため複数あるサブプログラムの定義方法の1つのみを扱った．このように定義されたサブプログラムは「内部手続き」と呼ばれるが，これとは別に「外部手続き」なるものも存在する．実際には外部手続きは後で学ぶモジュールを用いてモジュールの内部手続として実装する方が良いのだが，このやり方も知っておいて損は無い．特にモジュール化されていない外部ライブラリを用いる場合や，古いFortran 77仕様のプログラムを扱う際には(内部手続きが存在しないため)外部手続きの理解が必須となる．
7.6.1. 内部手続き
定義の仕方は 定義場所 を参照して欲しい．ここで注意しなければならないのは変数のスコープについてである．じつは， メインプログラム中で宣言した変数には内部手続きからアクセスすることが出来る (ただし逆は出来ない)ことに注意して欲しい．このことを表すサンプルが sample7.f90 である．これはそれほど長くないので，以下に全体を示す．
1program sample
2 implicit none
3
4 ! 16進数変換のためのテーブル (内部手続きからも参照される)
5 character(len=1), parameter :: hex_char(0:15) = &
6 & (/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', &
7 & 'A', 'B', 'C', 'D', 'E', 'F'/)
8
9 integer :: n = 10
10 character(len=8) :: hexstr
11
12 ! どちらの n を参照するか？
13 call sub()
14
15 ! 整数を16進数に変換して表示
16 n = 15 * 16**6 + 4 * 16**4 + 3 * 16**3 + 16**2 + 1
17 call decimal2hex(n, hexstr)
18 write(*, *) 'decimal = ', n, ' ===> hex = ', hexstr
19
20 stop
21contains
22 !
23 ! 内部手続きのスコープについて (1)
24 !
25 ! 内部手続きからはメインプログラムで宣言された変数を参照可能．ただし逆は不可．
26 !
27 ! nという名前の変数をサブルーチン内で宣言するかどうかで挙動が変わる
28 subroutine sub()
29 implicit none
30 ! もし以下の行があればメインプログラムのnとサブプログラムのnは独立
31 !integer :: n
32
33 write(*, *) n ! メインプログラム中の変数nにアクセス
34 endsubroutine sub
35
36 !
37 ! 内部手続きのスコープについて (2)
38 !
39 ! メインプログラムで定義された変数は内部手続きから参照出来るが，一般論としては
40 ! 引数として渡すようにした方が安全である．以下の例のようにプログラム全体で共通
41 ! に用いる定数であれば問題は起こらないことが多い．
42 !
43 ! 10進数を16進数に変換
44 subroutine decimal2hex(decimal, hex)
45 implicit none
46 integer :: decimal
47 character(len=*) :: hex
48
49 integer :: i, n, d
50
51 d = decimal
52 do i = 1, 8
53 n = d / 16**(8 - i)
54 d = d - n * 16**(8 - i)
55 ! メインプログラムで宣言された変数(hex_char)を参照
56 hex(i:i) = hex_char(n)
57 enddo
58
59 endsubroutine decimal2hex
60
61endprogram sample
33行目で内部手続き sub からメインプログラム中に定義された変数 n にアクセスしている．しかし，もし内部手続き sub 中で変数 n が定義されている場合(31行目のコメントを外した場合)には，この変数は sub 内部のみで有効な(メインプログラム中の n とは独立な)変数になる．一般的には，サブプログラムからメインプログラム中の変数を不用意に直接参照するのは間違いのもとになりやすい．それよりは，引数や返値を通じて値のやり取りを明示的に行う方が分かりやすいプログラムとなることが多い．
例外としてはプログラム全体で共通に用いる定数が挙げられる．定数の場合は参照されるだけなので，問題にならない場合が多い．例えば5-7行目で宣言している定数配列 hex_char を内部手続き decimal2hex の56行目で参照している．この場合は定数を参照するだけで，不用意に値が変更されることがないため，問題になることはあまりないであろう．
まとめると，メインプログラムで定義された変数に内部手続きからアクセスすることができるものの，気をつけて使わないと思わぬ動作を引き起こす可能性があるため， 基本的には関数の返値や関数・サブルーチンの引数を介してデータをやり取りする 方が間違いが少ないであろう．
7.6.2. 外部手続き(非推奨)
外部手続きは program から end program で囲まれた範囲 以外 に定義される．以下に示す例のように，定義する場所はメインプログラムの前でも後でもどちらでも良い．適切にコンパイル・リンクすれば別ファイルで定義したサブプログラムを用いることも可能である．
1!!!!!!!!!!! 外部手続きの定義場所 (1) !!!!!!!!!!
2function square_ext1(x) result(y)
3 implicit none
4 real(8) :: x
5 real(8) :: y
6
7 y = x**2
8
9 return
10endfunction square_ext1
11!!!!!!!!!!!
12
13! メインプログラム
14program sample
15 implicit none
16 !
17 ! 外部関数のinterface宣言
18 ! (本当は後で学ぶモジュールを使うほうがスマート)
19 !
20 interface
21 real(8) function square_ext1(x)
22 real(8) :: x
23 endfunction square_ext1
24 endinterface
25
26 !
27 ! 外部サブルーチンのinterface宣言
28 !
29 interface
30 subroutine sub_ext()
31 endsubroutine sub_ext
32 endinterface
33
34 ! 外部関数を呼び出すには実はこの書き方でも良いが，色々と問題が多いので非推奨
35 real(8), external :: square_ext2
36
37 ! 外部関数呼び出し
38 write(*, *) square_ext1(2.0_8), square_ext2(4.0_8)
39
40 ! 外部サブルーチン呼び出し
41 call sub_ext()
42
43 stop
44endprogram sample
45
46!!!!!!!!!!! 外部手続きの定義場所 (2) !!!!!!!!!!
47function square_ext2(x) result(y)
48 implicit none
49 real(8) :: x
50 real(8) :: y
51
52 y = x**2
53
54 return
55endfunction square_ext2
56
57subroutine sub_ext()
58 implicit none
59
60 write(*, *) 'sub_ext'
61
62 return
63endsubroutine sub_ext
64!!!!!!!!!!!
内部手続きとは異なり，外部手続きを使う場合にはメインプログラムで使用する外部手続きを明示的に宣言する必要がある．これは square_ext1 については20-24行目， sub_ext については29-32行目のように interface を用いて行なっている． interface では関数やサブルーチンの呼び出し形式のみ(引数のデータ型やその順番，関数の返値など)を宣言する [6]．
実際にはもう少しサボることが出来てしまう．関数については squre_ext2 は36行目で返値のみを external 属性をつけて宣言することで呼び出しができる．サブルーチンについては実は特に何も宣言しなくても呼び出しが可能である．しかし，ここでは interface を使って明示的に外部手続きを宣言すること強く推奨しておきたい．なぜなら，メインプログラムの外で定義された関数やサブルーチンについてはコンパイラが(引数の数や型などの)呼び出し形式を知る方法が無いため，間違った呼び出し方をしていてもコンパイルが通ってしまう．しかし不正な呼び出しをしているため，当然実行時にはエラーが発生してプログラムが異常終了することになる．一般的に実行時のエラーの方がコンパイルエラーよりも厄介でデバッグにも時間がかかるため，コンパイル時にチェックが可能な interface による宣言の方が良いのである [7]．(内部手続きでは文字通りメインプログラムの内部に定義されているので，コンパイラがメインプログラムをコンパイルする際に呼び出し形式のチェックが可能である.)
とにかく外部手続きは(行儀よく使おうと思うと)面倒なので，特に理由が無い限りは内部手続きを用いる方が良い．どうしてもメインプログラムの外で手続きを定義する必要がある場合には後で学ぶモジュールを用いる方が間違いが圧倒的に少ないのである．
7.7. 第7章 演習課題
7.7.1. 課題1
サンプルプログラムをコンパイル・実行して動作を確認せよ．さらに，適宜修正してその実行結果を確認せよ．
7.7.2. 課題2
与えられた倍精度実数 \(a (> 0)\) の平方根の近似値を返す関数を実装せよ．ただし平方根は以下のような逐次近似で計算するものとする．
ここで \(x_{n}\) は \(\sqrt{a}\) の \(n\) 番目の近似値である．初期値としては \(x_{0} = a\) を与え，反復は例えば \(\epsilon = 10^{-5}\) に対して， \(\|x_{n+1} - x_{n}\| < \epsilon \|x_{n}\|\) となるまで繰り返せば良い．（反復の途中結果は必要ないので \(x_{n}\) に配列を使う必要はないことに注意せよ．）
実装した関数と組み込み関数 sqrt の結果を比較し， \(\epsilon\) で与えた精度の範囲内で正しいことを確認すること．例えば，以下は標準入力から与えられた数値(この例では2.0)の平方根を計算して表示する例である．
1$ ./a.out
22.0
3sqrt(x) = 1.4142135623730951
4approx = 1.4142135623746899
7.7.3. 課題3
テストの点数が整数配列から与えられた時にヒストグラムを作成するサブルーチン histogram を実装せよ．例えば点数配列とビン幅を入力とし，作成されたヒストグラムの各ビンの中央値，各ビン内の人数を出力とする以下の様な形式のサブルーチンを作成すればよい．ただし，与える整数は \(0 \leq n < 100\) とするが，もしこの範囲を超えた入力があった場合にはエラーを表示して終了すること．
1subroutine histogram(score, binw, binc, hist)
2 implicit none
3 integer, intent(in) :: score(:) ! 点数(人数分)
4 integer, intent(in) :: binw ! ビンの幅(例えば10点)
5 real(8), intent(out) :: binc(:) ! ビンの中央値(例えば5, 15, ..., 95)
6 integer, intent(out) :: hist(:) ! 各ビン内の人数
7
8 ! ここでヒストグラムを作成
9
10end subroutine histogram
この場合のように固定幅のビンでヒストグラムを作成するのは簡単である． i 番目の点数がヒストグラムの j 番目の要素に入るとすると， j は
1j = score(i) / binw + 1
のように求められる(複雑なif文による分岐は必要ない!)．このインデックスが配列 hist の上限と下限の間に収まっていなければエラーとすれば良い．
このサブルーチンを用いて， score2.dat からデータを読み込みヒストグラム作成するプログラムを作成せよ．またその結果をgnuplotで図示せよ( with boxes を用いると良い)．出力結果は例えばビン幅を10とした場合には以下のようになる．
1$ ./a.out < score2.dat
2 5.0000000000000000 0
3 15.000000000000000 3
4 25.000000000000000 26
5 35.000000000000000 63
6 45.000000000000000 152
7 55.000000000000000 248
8 65.000000000000000 254
9 75.000000000000000 159
10 85.000000000000000 80
11 95.000000000000000 15
なお score2.dat の形式は 5章の課題3 の score1.dat と同じであるが，データは異なるものになっている．
同様に範囲外のデータを含む score3.dat を読みこませると
1$ ./a.out < score3.dat
2 Invalid input
のようなエラーを表示するように実装せよ．
7.7.4. 課題4
データ \(x_1, x_2, \ldots, x_n\) を大きさの順に並べ替える処理をソートという．ソートのアルゴリズムの中でも一番簡単なのがバブルソートと呼ばれるものである．これは以下の処理を \(n-1\) 回繰り返すことで実現される．
\(i=1, \ldots, n-1\) まで順に \(x_{i}\) と \(x_{i+1}\) の大小を比較し, \(x_{i} > x_{i+1}\) なら順番を入れ替える
この処理を \(k\) 回行うと \(x_{i}\) のうち \(k\) 番目に大きい要素が \(x_{n-k+1}\) に配置されるので， \(n-1\) 回繰り返すことで全ての要素を並び替えることが出来る．
ただし1回目の処理が終了した時点で最大値が \(x_n\) になっており，この要素については並べ替えの必要がない．従って2回目は \(x_1, \ldots, x_{n-2}\) までを処理すれば十分である．同様に \(m\) 回目の処理が終了した時点で \(x_{n-m+1}, \ldots, x_{n}\) までの位置は確定しているので， \(m+1\) 回目には \(x_1, \ldots, x_{n-m}\) までの処理を行えばよい．これによって比較回数を \((n-1)^2\) 回から半分に減らすことができる．
このバブルソートによって整数配列をソートするサブルーチン bsort を実装せよ．これは例えば以下のような形になるだろう．
1subroutine bsort(array)
2 implicit none
3 integer, intent(inout) :: array(:) ! 配列にはソートされた結果が代入される
4
5 ! バブルソート
6
7end subroutine bsort
作成したプログラムを用いて rand.dat のデータをソートし，結果が正しいことを確認せよ．ここでもデータファイルの形式は上の score2.dat と同一である．従って同様にリダイレクトで
1$ ./a.out < rand.dat
のように読みこませれば良い．
7.7.5. 課題5
データ \(x_1, x_2, \ldots, x_n\) からある値 \(X\) に等しいものを探しだす処理を探索という．その中でも以下は二分探索と呼ばれるアルゴリズムである．
まず \(x_i (i=1,\ldots,n)\) をソートする．
\(l = 1, r = n\) として以下の処理を繰り返す．
\(l > r\) ならば失敗(見つからなかったので処理を終える)．
\(m = (l + r)/2\) とし, \(X = x_{m}\) なら成功(見つかったので処理を終える)．
\(X > x_{m}\) なら \(l = m + 1\) ， \(X < x_{m}\) なら \(r = m - 1\) として[1]に戻る．
これをサブルーチン bsearch として実装せよ．
また open 文で rand.dat からデータを読み込み， bsort で配列をソートした後に bsearch で実際に適当な値の探索を行うプログラムを作成せよ．例えば標準入力から与えられた値をソート後の配列から探索を行い，結果を表示するようにすれば良い．この時の実行結果は以下のようになる．(この例では10や2004が入力値である．)
1$ ./a.out
2Input an integer : 10 # キーボード入力
3 10 was not found !
4
5$ ./a.out
6Input an integer : 2004 # キーボード入力
7 2004 was found at index 100
なおサブルーチン bsearch は例えば以下のような形式とすればよい．見つかった場合には探しだす値と等しい値が格納されている配列のインデックスを，見つからなかった場合には -1 を仮引数 idx に代入して返すようにするとよいだろう．
1subroutine bsearch(array, var, idx)
2 implicit none
3 integer, intent(in) :: array(:) ! ソートされた配列
4 integer, intent(in) :: var ! 探したい値
5 integer, intent(out) :: idx ! 見つかった要素へのインデックス
6
7 ! 二分探索
8
9end subroutine bsearch
7.7.6. 課題6
1行に1つの英単語が記述されているファイルを読み込み，英単語を辞書順にソートして出力するプログラムを作成せよ．空白行かファイルの終端に達するまで全ての単語を読み込むこと．ただし英単語は10文字以下，また読み込む単語の数は100個以下であると仮定してよい． words.txt を入力として作成したプログラムの動作を確認せよ．
実行結果は例えば以下のようになるだろう．(空白行の後の NOT_TO_BE_READ は読まれていないことに注意!)
1$ ./a.out < words.txt
2*** before sort ***
3tomato banana apple grape watermelon plum
4*** after sort ***
5apple banana grape plum tomato watermelon
なお，文字列の大小比較は辞書順となる(例えば 'apple' < 'banana')ので， bsort を一部修正するだけで文字列のソートが出来る．文字型変数の配列は例えば
1character(len=10) :: char_array(100)
などのように宣言すれば良い．ここで char_array は10文字分の文字列を保持する長さ100の配列である．
7.7.7. 課題7 †
与えられたファイルに含まれるアルファベットの各文字の出現回数をヒストグラムとして標準出力に表示するプログラムを作成せよ．ファイルの長さは任意とする．即ちファイルの終端まで全ての文字を読み込まなければならない(空白行があっても読み込み続ける)．ただし以下の条件を満たすこと．
アルファベット
a-z, A-Z以外は無視してよい．大文字と小文字は区別しない．
文字数があまりに多い時にはヒストグラムを適当に規格化すること．
wikipedia.txt (Wikipedia より引用)を処理した時には例えば以下のような出力となるだろう．この例では最大で60個の 'o' が出力されるように規格化してある．各アルファベットの後の括弧内の数字は文字数を表している．( wikipedia.txt で試す前に words.txt などの小さいデータで試した方が良いだろう．)
1$ ./a.out < wikipedia.txt
2A( 110):oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
3B( 9):ooooo
4C( 44):oooooooooooooooooooooooo
5D( 32):ooooooooooooooooo
6E( 99):oooooooooooooooooooooooooooooooooooooooooooooooooooooo
7F( 31):ooooooooooooooooo
8G( 38):ooooooooooooooooooooo
9H( 32):ooooooooooooooooo
10I( 87):ooooooooooooooooooooooooooooooooooooooooooooooo
11J( 2):o
12K( 2):o
13L( 38):ooooooooooooooooooooo
14M( 49):ooooooooooooooooooooooooooo
15N( 91):oooooooooooooooooooooooooooooooooooooooooooooooooo
16O( 81):oooooooooooooooooooooooooooooooooooooooooooo
17P( 39):ooooooooooooooooooooo
18Q( 0):
19R( 93):ooooooooooooooooooooooooooooooooooooooooooooooooooo
20S( 63):oooooooooooooooooooooooooooooooooo
21T( 79):ooooooooooooooooooooooooooooooooooooooooooo
22U( 37):oooooooooooooooooooo
23V( 14):oooooooo
24W( 5):ooo
25X( 1):o
26Y( 15):oooooooo
27Z( 0):
組込み関数 ichar ， char を用いるとよい．また規格化の際には四捨五入をする関数 nint を用いることが出来る．
特に仕様は限定しないが，関数やサブルーチンを使った分り易いプログラムを目指すこと．