数値解析基礎・演習 Fortran90-(4) subroutine の利用 ver.2011.05

数値解析基礎・演習 Fortran90-(4) subroutine の利用 ver.2011.05
情報環境学専攻（社会・環境系） 斎藤隆泰
1.
はじめに
17
end program sub1
これまでに if 文や do 文も学び、少しずつであ
------------------------------------------------------------------
るが Fortran90 の文法の知識も広がってきた。本
このプログラムは 1 つのプログラム中で、1～5 ま
講義名は数値解析基礎・演習である。本格的な数
での整数の和と 1～10 までの整数の和を計算し、
値計算を Fortran で行うには今回学ぶ subroutine
その結果を出力するものである。6～10 行目で 1
や次回学ぶ配列の知識は特に重要な文法の一つで
～5 までの和を、12～16 行目で 1～10 までの和を
ある。そこで、今回はその subroutine について講
求めている。また、8 行目と 14 行目に登場する
義の補足をしておく。性能の良い subroutine を開
dble は dble(i)で整数 i を倍精度実数へと変換して
発することは、大規模な 1 つのプログラムをグル
いる。
ープで開発するときにも都合が良いし、作成した
さて、do 文の説明で述べたが、できるだけ同じ
subroutine はプログラマー個人の大切な財産とも
ような繰り返し計算はプログラムで自動化した方
成り得る。
が、手間が省けて良い。そこで、6～10 行目と 12
～16 行目のように同じような総和の計算をして
2. subroutine の役割
subroutine とは何なのか？利用するとどういう
いる部分を 1 つにまとめることはできないだろう
か？subroutine を用いれば、先のプログラムを次
メリットがあるのかを、この節で考えてみよう。
のように書き換えることが可能である。
そのために、まず下記のプログラムを作成し、実
------------------------------------------------------------------
行してみよう。
1
program sub2
------------------------------------------------------------------
2
implicit none
1
program sub1
3
integer::i
2
implicit none
4
real(8)::x,y
3
integer::i
5
4
real(8)::x,y
6
5
7
6
x=0.0d0
8
7
do i=1,5
9
8
call sum(5,x)
x=x+dble(i)
10
write(*,*)'x=',x
write(*,*)'y=',y
9
end do
11
10
write(*,*)x
12
11
call sum(10,y)
13
end program sub2
12
y=0.0d0
14
13
do i=1,10
15
subroutine sum(n,x)
16
implicit none
14
y=y+dble(i)
15
end do
17
integer::n,i
16
write(*,*)y
18
real(8)::x
どのようにとってもよい。上の例では最初の引数
19
x=0.0d0
のペアは 5 と n、2 番目の引数のペアは x と x で
20
do i=1,n
それぞれ整数、倍精度実数同士が引き渡しされて
21
22
x=x+dble(i)
end do
23
24
いる。
そのため、
6 行目の第一引数である 5 を 5.0d0
としてはならない。5.0d0 は倍精度実数型である
が、対応する subroutine の引数 n は整数で宣言し
end subroutine sum
-----------------------------------------------------------------ここで初めて出てくる文法は 6，8 行目の call
文と 15 行目の subroutine 文である。program
ているからである。また、引数は原則 1 対 1 に対
応しなければならないため、call 文と subroutine
文の引数の数は同じでなければならない（ただし、
一部の例外はあり得る）
。
sub1 との違いは、1～5 まで、1～10 までの和を
つまり 6 行目の call 文で整数 5 と倍精度実数 x
計算するプログラムを直接作るのではなく、1～
を subroutine へ引数として渡し、15 行目で n=5
n(n はある整数)までの和を計算する subroutine
を受け取り、20 行目からの do 文で i=1,5 の和を
プログラムを作成し、それを n=5 or 10 として 2
計算した後、計算結果である x の値が引数の 1 対
回用いている点にある。つまり、15 行目からの
1 の対応関係により program sub2 の x として引
subroutine は、subroutine に n=5 or 10 を与えて
き渡されているのである。同じ事は 8 行目でも行
n までの和を計算させる、いわば関数のような役
われている。8 行目では第一引数を 10 としている
目を担っている。ここで、subroutine と call 文の
ので、subroutine の中では n=10 における和の計
使い方について説明しよう。例えば 6 行目と 15
算が実行され、実行結果は、引数の引き渡しによ
行目の使用例に注目する。
り第二引数の y へ引き渡されて program sub2 へ
戻る。
6
call (5, x)
このように一度 subroutine を作成すれば、その
subroutine をプログラム中のどこからでも呼び出
15
subroutine(n, x)
し、適切な引数を与えることで、所望の計算結果
を得ることができる。
6 行目や 15 行目の call 文や subroutine 文の括
ここでは program sub2 から subroutine sum を
弧内の変数等を引数と呼ぶ。引数は整数、実数、
呼び出しているが、program sub2 をメインプログ
複素数、具体的な数値を直接カンマで区切りなが
ラムと呼ぶ。メインプログラムは通常 program
ら連続して書くことができる。6 行目では 5 と x
sub2 や sub1 のように program 文から始まる（モ
を引数として subroutine を呼び出している。この
ジュール等を除く。モジュールについては講義で
ように、subroutine を利用するには call 文を用い
触れない）呼び出し側のプログラムを指す。メイ
れば良い。ただし、call 文の引数は subroutine の
ンプログラムから subroutine を呼び出すことは
引数へ、それぞれ上の例の矢印のように 1 対 1 に
できるが、subroutine からメインプログラムを呼
対応し引き渡されることに注意する。このとき、
び出すことはできない。ただし、subroutine の中
それぞれ対応する引数の型（integer, real 等）は
から、さらに別の subroutine を呼び出すこともで
同じでなければならない。このように、call 文と
きる。
subroutine 文の引数が 1 対 1 で対応していれば、
このように subroutine を用いれば、同じような
引数の並びは call(x,5),subroutine(x,n)のように、
処理が必要となる計算を 1 つの subroutine を使う
ことで何回も計算することができ、非常に便利で
18
ある。大規模なプログラム開発では、例えば、あ
19
z=dcmplx(a,b)
る担当者は 1 から任意の整数までの和を求める
20
write(*,*)cdabs(z)
subroutine を作成し、別の担当者は任意の整数 n!
21
を求める subroutine を作成する・・・等、作業を
22
分担することも多い。しかも一度正しい
23
subroutine を作れば、別の program で必要とな
24
subroutine absolute(a,b,c)
った時に、再度その subroutine を使えば良く、
25
implicit none
様々な点で多くのメリットがある。
26
real(8)::a,b,c
end program comp
27
3．subroutine を作成する例題
subroutine について少し知識がついたところで、
28
c=dsqrt(a**2+b**2)
29
別の例題を解いてみよう。次の例題は、ある複素
30
end subroutine absolute
数の実部と虚部を引数とし、その複素数の絶対値
------------------------------------------------------------------
を計算する subroutine を作成し、メインから呼び
出すプログラムである。できるだけ何も見ないで
ここでは複素数の扱いについても触れる。4 行
自分で考えながらやってみよう。失敗しても、間
目の complex*16 は実部、虚部共に倍精度で表さ
違った文法を書いたとしても、単に修正してコン
れた複素数を利用するための宣言文である。ここ
パイルすれば良いだけの話である。
では任意の複素数を扱うために、複素変数 z を使
うための宣言をしている。8 行目で実部として a
------------------------------------------------------------------
を、12 行目で虚部として b を読む。14 行目で実
1
program comp
部 a と虚部 b、そして subroutine で計算させる絶
2
implicit none
対値 c を引数として call 文を用いて subroutine
3
real(8)::a,b,c
を呼び出す。ここでは、subroutine 名を absolute
4
complex*16::z
とした。引数の受け渡しは、14 行目の call 文と
5
24 行目の subroutine 文の引数を 1 対 1 に対応さ
6 !<----real part---->
せ、それらの引数の型を合わせる。28 行目で絶対
7
write(*,*)'input real part'
値を計算し、その結果を左辺の c に代入する。c
8
read(*,*)a
は再び 14 行目の引数 c としてメインプログラムに
9
戻 り 、 そ の 結 果 を 17 行 目 で 出 力 し て い る 。
10 !<----imaginary part---->
subroutine の仕組みの理解が深まったであろう
11
write(*,*)'input imaginary part'
か？
12
read(*,*)b
13
14
subroutine の使い方に加えて、ここでは複素数
の 扱 い に つ い て も 学 ぼ う 。 19 行 目 で
call absolute(a,b,c)
15
z=dcmplx(a,b)としているが、この文の意味は数学
で言う z=a+bi であり、実数 a,b を実部、虚部とし
16
write(*,*)'absolute value='
た複素数を複素数で宣言した変数 z に代入すると
17
write(*,*)c
いう意味である。このように実部と虚部をカッコ
でくくり、dcmplx を付けることによって倍精度の
subroutine 名前で終わる。
複素数を定義出来る(Appendix A も参照)。定義し
実は、
Fortran90 で subroutine を呼び出す場合、
た複素数を左辺の倍精度複素変数 z へ代入してい
intent 属性なるものを付ける必要性がでる場合が
る。さらに、20 行目であるが、cdabs(z)で複素変
ある。intent 属性は Fortran77 ではなかった文法
数 z の絶対値を計算できる。実はこうすれば、わ
であり、ここでは詳しく述べない。ここで説明し
ざ わ ざ 24 ～ 30 行 目 ま で の 絶 対 値 を 計 算 す る
たことをとりあえずきちっと自分のものにするこ
subroutine を 作 る 必 要 が な い 。 こ の よ う に 、
とが大事である。
Fortran では、数値計算で用いるいくつかの重要
な関数が、最初から用意されている。これは
Appendix A 変数の型変換
「Fortran90-(1)基本的な使い方」で少し紹介した
次のプログラムを作成し、実行してみよう。次
組み込み関数の 1 つである。組み込み関数にはベ
のプログラムは π=3. 1415926535 8979323846
クトルの内積等、様々なものが準備されているの
264338327・・・・・を Fortran の組み込み関数
で、興味がある学生は、自分で調べてみて欲しい。
を使って求めるプログラムである。
3
------------------------------------------------------------------
subroutine 利用のまとめ
今回は subroutine の利用方法について簡単に
1
program variable
説明した。subroutine を利用する方法、気をつけ
2
implicit none
ねばならない基本的な点をまとめておく。
3
integer::i
4
real::pi,a
5
real(8)::pid,ad
call 名前（a, b, c, ・・・・）
6
で利用する subroutine を呼ぶ。名前の部分は呼び
7
a=1
出したい subroutine 名をつける。
8
ad=1.0d0
subroutine の 1 行目は必ず、
9
write(*,*)a
10
write(*,*)ad
subroutine 名前（a, b, c, ・・・・）
11
12
pid= atan(1.0)*4.0
のようにし、subroutine の名前は call 文における
13
write(*,*)pid
名前と等しいかを確認する。また、引数も１対１
14
に対応させる。そのため、call 文と subroutine 文
15
pid=datan(1.0d0)*4.0d0
では引数の数は原則同じとなる。また、対応する
16
write(*,*)pid
引数の型はそれぞれ同じものでなければならない
17
（ 例 え ば 、 先 の program comp で は call
18
文,subroutine 文の先頭はいずれも a であるが、メ
------------------------------------------------------------------
end program variable
イン、subroutine どちらにおいても倍精度実数の
real(8)で宣言されている）。subroutine でもメイ
このプログラムを著者の環境（intel fortran コン
ンプログラムと同様、implicit none から始め、必
パイラー）でコンパイルし、実行すると、次の結
要 な 変 数 を 宣 言 し 、 subroutine の 最 後は end
果を得る。
結果を左辺の倍精度実数で定義された変数に代入
実行結果 -----------------------------------------------------
したとしても、もとの計算が単精度実数であるた
1.000000
めに、9 桁以降の計算が正しく計算されている保
1.00000000000000
証がないのである。そのため、円周率の計算結果
3.14159274101257
が単精度分の 8 桁までしか合わないのである。
3.14159265358979
------------------------------------------------------------------
先に組み込み関数 atan の倍精度を datan と述
べたが、Fortran では組み込み関数名の前に倍精
度(double)の頭文字 d を付けると、倍精度実数型
12、15 行目の atan や datan は逆正接(arctan：ア
で計算される組み込み関数、c を付けると複素数
ークタンジェント)を計算する組み込み関数であ
の計算を表す場合が多い。つまり、組み込み関数
り、これを利用することによって、円周率πを計
atan に d をつけた datan は倍精度実数で計算され
算している(atan と datan の違いは後に説明する)。
た結果を返すのである。
さて、実行結果を見てみよう。7 行目で 1 を、8
例えば、
行目で 1.0d0 をそれぞれ a, ad に代入し、その結
果を 9,10 行目で書かせているが、結果は ad の方
（単精度）実数
倍精度実数
が 16 桁表示されていることに気付く。この違いは、
sin(x)
dsin(x)
宣言文にある。4 行目では a を実数(real)で宣言し
cos(x)
dcos(x)
ているが、5 行目では ad を倍精度実数(real(8))で
sqrt(x)
dsqrt(x)
宣言している。つまり倍精度実数では通常の実数
の倍の桁数を扱えることわかる。実数では 8 桁、
のように、左側の組み込み関数は単精度、右側の
倍精度実数は 16 桁分 1.00・・・・が表示されて
組み込み関数は倍精度で行われた演算結果を返す。
いる。なお、ここで単に実数とは単精度実数のこ
他の例を挙げよう。program comp の 19 行目は
とを指す。
z=dcmplx(a,b)
これを踏まえて、12、15 行目の結果を見てみよ
う。13，16 行目で計算結果が代入される pid は倍
精度実数型(real(8))で宣言されているので、結果
であるが、これは z を倍精度複素数型としている
はどちらも 16 桁表示されているが、それぞれの結
からである。z を単精度とし、a,b を単精度実数と
果を見比べると 8 桁までしか両者は合っていない。
するなら、dcmplx の d を取って z=cmplx(a,b)と
しかも、13 行目の結果では 9 桁以降は、問題の最
すれば良い。
初で示した実際の円周率と合っていない。この違
このような事実を考えると、本講義では、倍精
いはなぜ生じたのであろうか？それは、右辺の計
度実数や倍精度複素数で宣言した変数に関係する
算が、単精度実数で行われているか、倍精度実数
演算には、単精度の計算と区別するために、定数
で行われているかによる。12 行目の右辺は実数で
には 2.0d0 等、d0 を付けること、組み込み関数を
計算されている。そのため、atan の引数は 1.0、
伴う場合は d を付けて倍精度演算を行うことを薦
その後の乗算も 4.0 としている。しかし 15 行目は
める。
datan（atan の倍精度計算）を用いて引数を 1.0d0、
乗算も 4.0d0 としている。つまり、12 行目の右辺
は通常の単精度実数で計算しているために、計算