はじめに
前回の記事では、ソースコードをコンパイルするためのコンパイラのインストールを行い、「Hello, World!」をターミナル上に出力しましたが、これだけでは全く使い物になりません。
これからいくつかの記事に分けてFortranプログラミングを紹介していきますが、今回は、必ず知っておかなければならない繰り返し処理のための doループおよびファイルの入出力について説明します。
説明だけ見てもよくわからないと思いますので、1/4円の積分を使って円周率を求めるプログラムを例にとって説明します。
macOS Monterey Version 12.7.4
Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中
計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています!
理論と実装の溝を埋めてくれる良書です!
doループとopen文
doループ
doループの基本的な形は次のとおりです。
do カウンタ = 最初の整数, 最後の整数, ステップ
なにか処理を書く
end do
do と end do ではさまれた部分に繰り返したい処理を記述します。インデントは不要ですが、可読性向上のために、インデントするようにします。
open文
ファイルの入出力の前には、open文を使ってファイルをopenする必要があります。よく使う形を示します。
open(unit=入出力ファイル番号, iostat=open文のステータス, file=ファイル名 &
& action=writeやreadなどのファイルに対して行う操作)
各引数の意味は下で示すプログラムで明らかになると思います。
1/4円の積分により円周率を求めるプログラム
円の方程式の第一象限部分
$$ y = \sqrt{1-x^2} $$
と $x$ 軸で囲まれた部分の面積は $\pi/4$ です。この性質を使って、数値積分により円周率の近似値を求めてみたいと思います。
分割数 $n$ を決めて、区間 $0\le x\le 1$ を $n$ 分割し、各分割区間の左端の値と分割の大きさをかけ合わせたもの(長方形の面積)を足し合わせます。
$$ S = \sum_{i=0}^{n-1} y_i\Delta x $$
ここで、
$$ \Delta x=1/n,\quad y_i = \sqrt{1-(i\Delta x)^2} $$
です。これよりも優れた方法(台形則、シンプソン則)はよく知られているところですが、最も簡単なのでこれを使います。絵を書いてみるとわかるように、真の円周率よりも大きな値が計算されるだろうと予想されます。
入力ファイルとして分割数 $n$ が書かれたファイルを読み、上の方法で積分を計算し、結果をファイルに出力するプログラムを作成します。
プログラムの実装
上で説明したもののFortranでの実装例を以下に示します。エディタで次のようなファイルを作成し、適当な場所(~/Desktop/LabCode/fortran/02_compute-pi
とします)に適当な名前(ここでは、compute-pi.f90とします)で保存してください。
program main
implicit none
! declare variables
integer(4) :: i
integer(4) :: n
integer(4) :: lin = 60, lout = 61
integer(4) :: iost
real(8) :: x, y
real(8) :: dx
real(8) :: sum, pi
! read the input file
open(lin, file="input.dat", action="read", iostat=iost)
if (iost /= 0) then
write(*, *) "ERROR: cannot open file!"
stop
end if
read(lin, *) n
if (n < 0) then
write(*, *) "ERROR: the number of divisions specified is not positive!"
stop
end if
write(*, *) "the number of divisions is ", n
close(lin)
! rough estimation of pi
dx = 1.0d0 / dble(n)
sum = 0.0d0
do i = 0, n-1
x = dx * dble(i)
y = sqrt(1.0d0 - x*x)
sum = sum + dx*y
end do
pi = sum * 4.0d0
! save result
open(lout, file="output.dat", action="write")
write(lout, "(i8, x, es12.6)") n, pi
close(lout)
! end of program
stop
end program main
コンパイル
保存したら、コンパイルしましょう。まず、プログラムを保存した場所に移動します:
$ cd ~/Desktop/LabCode/fortran/02_compute-pi
続いて、ソースファイルをコンパイルします。次のように打ち込んでください:
$ gfortran compute-pi.f90 -o compute-pi
コンパイルが正常に終了すると、実行ファイル compute-pi ができているはずです!
入力ファイルの準備
次に、入力ファイル input.dat を作っておきましょう。中身は、分割数だけが書かれたファイルです。n=1000としたければ、
10000
と一行だけ書かれたファイルを作成し、名前を input.dat とし、プログラムを実行するディレクトリ(今回はプログラムを置いた場所と同じ~/Desktop/LabCode/fortran/02_compute-pi
)に保存してください。
プログラムの実行
プログラムを実行してみましょう。次のようにターミナルに入力してください:
$ ./compute-pi
ターミナル上に 「the number of divisions is 1000」と表示され、output.dat が生成されて中身が
1000 3.143555E+00
となっていれば成功です!
予想通り、円周率に近い値が計算されているものの、やや大きな値となっています。
input.dat がない場合、 ターミナル上に「ERROR: cannot open file!」と表示され、プログラムが終了します。input.datに0以下の値が書かれている場合には、「ERROR: the number of divisions specified is not positive!」と表示されます。試してみてください。
コードの解説
上に書いたソースプログラムの program main
から end program main
で囲まれた部分の解説をします。
implicit none
! declare variables
integer(4) :: i
integer(4) :: n
integer(4) :: lin = 60, lout = 61
integer(4) :: iost
real(8) :: x, y
real(8) :: dx
real(8) :: sum, pi
変数の宣言をします。implicit none
は暗黙の型宣言を無効化するもので、ほとんど必ず使用します。
ここで示しているように、同じ型の変数は同時に宣言することもできますし、変数宣言と同時に値を代入しておくこともできます。
可読性のために、同じような役割のものをブロック状に配置しています。
! read the input file
open(lin, file="input.dat", action="read", iostat=iost)
if (iost /= 0) then
write(*, *) "ERROR: cannot open file!"
stop
end if
read(lin, *) n
if (n < 0) then
write(*, *) "ERROR: the number of divisions specified is not positive!"
stop
end if
write(*, *) "the number of divisions is ", n
close(lin)
この部分ではファイルの入力をしています。lin
はファイル番号で、変数宣言のところで 60
が代入されています。このファイルに関する操作 (read
, close
) はこのファイル番号を指定して行います。ファイル名は input.dat であることを想定しています。このファイルに対して行いたいことは読み込みですから、action
には “read"
を指定します。
まず、ファイルが存在するかどうかのチェックを行っています。ファイルが正しく開けたら iost
に 0
が代入されます。存在しない場合には 0
以外の値が入りますので、メッセージをターミナルに出力してプログラムを終了します。
無事にファイルを開くことができれば、中身を読んで n
に代入します。これは read 文で行います。read 文の第一引数にはファイル番号、第二引数は書式を入れます。n
は正の整数でなければならないので、そうでない場合にはメッセージをターミナルに出力してプログラムを終了します。
すべてうまく行っている場合には、読み込んだ n
の値を出力してファイルを閉じます。
! rough estimation of pi
dx = 1.0d0 / dble(n)
sum = 0.0d0
do i = 0, n-1
x = dx * dble(i)
y = sqrt(1.0d0 - x*x)
sum = sum + dx*y
end do
pi = sum * 4.0d0
この部分は、上で述べた長方形の面積の足し合わせで積分を近似し、円周率の近似値を求めているところです。
「doループとopen文」のセクションで述べたように、do
と end do
で囲まれた部分が繰り返されます。<ステップ>
は省略可能で、省略した場合 1
になります。0からn-1まで繰り返し、sum
に足してゆきます。
ところどころ数値に d0
がつけられていますが、これは数値が倍精度実数で $\times 10^0$ を表しています$\times 10^3$ を表したければ d3
とします。これは倍精度実数の数値であるときは必ずつけなければなりません。
また、dble(n)
や dble(i)
は () の中身を倍精度実数に変換する関数です。「倍精度実数 / 整数 = 倍精度実数」、「倍精度実数 × 整数 = 倍精度実数」となるため、なくても結果は変わりません(少なくとも私の処理系では)が、型を意識したプログラミングのためにはあったほうが良いでしょう。
! save result
open(lout, file="output.dat", action="write")
write(lout, "(i8, x, es12.6)") n, pi
close(lout)
結果をファイルに出力します。ファイル番号は変数宣言のところで代入した 61
です。ファイル名は output.dat とし、書き込みなので action
を “write”
とします。
write文の第一引数にはファイル番号、第二引数は書式を入れます。"(i8, x, es12.6)"
の意味は「8桁分の文字数を確保して整数を左詰め、スペース、12桁分のスペースを確保して小数点以下6桁を科学表記の指数形式で表示する」というものです。科学表記というのは仮数が1以上10未満になるような指数表記方法です。指定したファイルに、指定した書式で n
と pi
を出力します。
書き込みが終わったらファイルを閉じます。
! end of program
stop
end program
が stop
を兼ねることが多いらしいのですが、そうでない場合のために stop 文を最後に書いておきます。
最後に
今回は、数値計算をする上で必須のdoループとファイルの入出力について紹介しました。今後もFortranを使ったプログラミングの基礎について紹介していきたいと思います。
Pythonで学ぶ偏微分方程式の数値シミュレーションの技術書を販売中
計算方法、グラフやアニメーションの作成方法、計算ログの残し方を惜しみなく解説しています!
理論と実装の溝を埋めてくれる良書です!