【Fortran】doループ、ファイルの入出力

【Fortran】doループ、ファイルの入出力

はじめに


前回の記事では、ソースコードをコンパイルするためのコンパイラのインストールを行い、「Hello, World!」をターミナル上に出力しましたが、これだけでは全く使い物になりません。

これからいくつかの記事に分けてFortranプログラミングを紹介していきますが、今回は、必ず知っておかなければならない繰り返し処理のための doループおよびファイルの入出力について説明します。

説明だけ見てもよくわからないと思いますので、1/4円の積分を使って円周率を求めるプログラムを例にとって説明します。

動作検証済み環境

macOS Monterey Version 12.7.4

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" を指定します。

まず、ファイルが存在するかどうかのチェックを行っています。ファイルが正しく開けたら iost0 が代入されます。存在しない場合には 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文」のセクションで述べたように、doend 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未満になるような指数表記方法です。指定したファイルに、指定した書式で npi を出力します。

書き込みが終わったらファイルを閉じます。

  ! end of program
  stop

end programstopを兼ねることが多いらしいのですが、そうでない場合のために stop 文を最後に書いておきます。

最後に


今回は、数値計算をする上で必須のdoループとファイルの入出力について紹介しました。今後もFortranを使ったプログラミングの基礎について紹介していきたいと思います。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です