第2版 update: 2021年10月 葛西 真寿 弘前大学大学院理工学研究科
Jupyter Fortran Kernel を使った Fortran 90/95 プログラミング入門。 はじめての Python プログラミングの内容を Fortran で書いてみたものです。
まずは,以下のような掛け算の答えを出力するプログラムを実行してみます。
program ai1
boku = 2.5
kanojo = 0
love = boku * kanojo
print *, love
end program
エラーなく実行されたが,結果が 0 というのは,あまりに寂しい。
以下のように少し変更して実行すると...
program ai2
watashi = 2.5
kareshi = 0.4
love = watashi * kareshi
print *, love
end program
おかしい... 2.5 に 0.4 をかければ 1.0 になるはずだが...
暗黙の型宣言という歴史的慣習
歴史的慣習により,FORTRAN では,黙っていれば i, j, k, l, m, n
で始まる変数は 整数,それ以外のアルファベットで始まる変数は 実数 として扱われてきた。
上記の例題 program ai2
では,kareshi = 0.4
は実は kareshi = 0
のこと。
このような思いもかけないトラブル(アイにトラブルはつきもの?)を防ぐ意味も兼ねて,最近の教科書では,この暗黙の型宣言を無効にするように,プログラム文の冒頭に
implicit none
と書くことを推奨している。
program ai
! ! から始まる行はコメント行。実行に影響しません。
! 暗黙の型宣言は無効にして,
! プログラムで使う変数は全て宣言します。
implicit none
real :: watashi, kareshi, love
watashi = 2.5
kareshi = 0.4
love = watashi * kareshi
print *, love
end program ai
上のプログラム ai を例に,現代的視点から推奨されている書き方を以下に例示します。
program ...
ではじまります。!
からはじまる行はコメント行です。コメント行はそれを読む人間のための注釈であり,プログラムの実行時には無視されます。implicit none
と宣言したら,以降はプログラム文内で使用する変数は,全て以下の例のように宣言する必要があります。real
(実数)として宣言します。2.5
を変数 watashi
に代入しています。人一倍どころか,人より 2.5
倍の思いがあることを示しています。0.4
を変数 kareshi
に代入しています。watashi
に代入された数値と,変数 kareshi
に代入された数値をかけた答えを,変数 love
に代入しています。love
の値を出力しています。標準出力(画面)への出力は print
に統一します。(print *
の *
については「書式指定」のところで。)end program ...
で終わります。エンゲル係数とは,家計の消費支出に占める飲食費(食料費,食費)の割合(パーセント単位)です。
消費支出 110,000円,飲食費 37,000円の人のエンゲル係数を求めるプログラム例。最低限,以下のようにすれば答えが出ます。
print *, 37000 * 100. / 110000
end
ここからは,単に電卓的に使うだけではなく,応用・発展も見据えて以下のようなプログラム文を作成して実行することにします。
以下の例では,
spending
に消費支出である 110000
(円)という数値を代入し,food
に飲食費である 37000
(円)という数値を代入し,eng
に代入し,eng
を表示させています。program engel
! エンゲル係数を求める
implicit none
! 消費支出 spending, 飲食費 food, エンゲル係数 eng
real :: spending, food, eng
spending = 110000
food = 37000
eng = food * 100 / spending
print *, eng
end program engel
program real
implicit none
real :: r
print *, '1 / 3 = ', 1.0 / 3.0
print *, 'real の最大値は', huge(r)
print *, 'real の最小値は', tiny(r)
end program real
real
では表せない非常に大きい実数や,より有効桁数の多い実数を使う場合は,以下のように real(8)
や real(16)
と宣言します。表示できる最大値・最小値は有効桁数は以下のようにして確認できます。
program real816
implicit none
real(8) :: r8
real(16) :: r16
print *, '1 / 3 = ', 1.0_8 / 3.0_8
print *, 'real(8) の最大値は', huge(r8)
print *, 'real(8) の最小値は', tiny(r8)
print *, ' '
print *, '1 / 3 = ', 1.0_16 / 3.0_16
print *, 'real(16) の最大値は', huge(r16)
print *, 'real(16) の最小値は', tiny(r16)
end program real816
例えば有効桁数 12 桁の実数を表したい場合,real
, real(8)
, real(16)
のどれを使えばいいか調べるのが面倒な場合は
integer, parameter :: dp = selected_real_kind(12)
real(dp)
のようにします。
program printau1
implicit none
! 天文単位は 12 桁の定義定数(誤差無し)
integer, parameter :: dp = selected_real_kind(12)
real(dp), parameter :: au1 = 149597870700.0_dp
print *, "1 天文単位は",au1," m です。"
end program printau1
program int
implicit none
integer :: i
print *, '1 / 3 = ', 1 / 3
print *, 'integer の最大値は', huge(i)
end program int
より桁数の多い整数を扱う場合は以下のように integer(8)
や integer(16)
を使います。表示できる最大値は有効桁数は以下のようにして確認できます。
program int816
implicit none
integer(8) :: i8
integer(16) :: i16
print *, 'integer(8) の最大値は ', huge(i8)
print *, 'integer(16) の最大値は', huge(i16)
end program int816
例えば有効桁数 12 桁の整数を表したい場合,integer
, integer(8)
, integer(16)
のどれを使えばいいか調べるのが面倒な場合は
integer, parameter :: ik = selected_int_kind(12)
integer(ik)
のようにします。
program printau
implicit none
! 天文単位は 12 桁の定義定数(誤差無し)
integer, parameter :: ik = selected_int_kind(12)
integer(ik), parameter :: au = 149597870700_ik
print *, "1 天文単位は",au," m です。"
end program printau
program intchar
implicit none
! 文字型の変数 moji(長さ12)の宣言
character(12) :: moji
moji = 'engel keisu '
print *, moji
end program intchar
優先的に計算したい箇所は丸括弧 ( )
で グループ化します。
program yusen
print *, 2 + 50 - 5 * 6 /2
print *, 2 + (50 - 5 * 6)/2
end program yusen
//
は文字列の連結演算に使います。以下の例を参照。
program charadd
implicit none
character(6) :: mo, ji
mo = 'engel '
ji = 'keisu '
print *, mo // ji
end program charadd
文字列と数値を直接連結することはできません。
文字列と数値を連結などしなくても,print
文は,カンマ (,
) で区切ることによって複数の項目(変数)を表示することができます。
program engel2
implicit none
real :: spending, food, eng
spending = 110000
food = 37000
eng = food * 100 / spending
print *, 'エンゲル係数は ', eng, 'です。'
end program engel2
除算 /
の例を以下に示します。
program divide
print *, '17 / 3 = ', 17 / 3
print *, '17.0 / 3 = ', 17.0 / 3
print *, '17.0_8 / 3 = ', 17.0_8 / 3
print *, '17.0_16 / 3 = ', 17.0_16 / 3
end program divide
2行目では,整数を整数で割るので,小数分を切り捨てた整数値を返します。
3行目では,17.0
と最後に小数点をつけたので,実数型(real
)の値を返します。この例では,小数点以下6桁くらいまでの精度です。
4行目は 17.0_8
で real(8)
の実数としています。この例では,小数点以下15桁くらいまでの精度です。
5行目では,17.0_16
でreal(16)
の実数としています。この例では,小数点以下33桁くらいまでの精度です。
べき乗は **
です。
以下の例では $2^4 = 2 \times 2 \times 2 \times 2$ と $2^{1/2} = \sqrt{2}$ の値を表示します。
program power
print *, '2**4 =', 2**4
print *, '2**0.5 =', 2**0.5
end program power
べき乗に限らず,あまりにも桁数が多くなる計算すると,基本整数型(integer
)の範囲を超えてしまい,おかしな値を出力することがあるので注意。
以下の例では,$2^{31}$ あたりで答えが変になっています。
program oor
! out of range
print *, '2**29 = ', 2**29
print *, '2**30 = ', 2**30
print *, '2**31 = ', 2**30 * 2
print *, '2**32 = ', 2**30 * 2 * 2
end program oor
非常に大きい(または極めて小さい)数値を扱う場合や,real
の有効桁数約6桁では足りない場合,例えば有効桁数を15(以上)とする場合は,以下のようにします。
上の例で整数型の場合に正しい答えが出なかった $2^{31}, 2^{32}$ なども,以下の例のように有効桁数を多くとることによって精度よく出力されます。
program precision
implicit none
! 15桁(以上)の有効桁数を持つ実数・整数の定義用
integer, parameter :: dp = selected_real_kind(p=15)
integer, parameter :: ik = selected_int_kind(15)
real(dp) :: sqrtwo
sqrtwo = sqrt(2.0_dp)
print *, 'sqrt(2.0) = ', sqrt(2.0)
print *, 'sqrt(2.0_dp) = ', sqrtwo
print *, '2**31 = ', 2_ik**31
print *, '2**32 = ', 2_ik**32
end program precision
1 天文単位(1 au)の距離を光速 c で進むのに要する時間はいくらか。ただし,
149597870700
m299792458
m/s である。参考:
上記の例では,プログラムの中で消費支出や飲食費を決めていましたが,今度はプログラムを実行する人が入力できるようにしてみます。
以下のプログラム例では,read *, spending
でキーボードからの入力値を変数 spending
に入れます。
program engelread
implicit none
real :: spending, food, eng
print *, "消費支出? "
read *, spending
print *, "飲食費? "
read *, food
eng = food * 100 / spending
print *, "エンゲル係数は", eng
end program engelread
Fortran kernel は Jupyter Notebook 環境で Fortran を実行できるので大変便利なのですが,残念なことに read
文には対応していないため,上記のようなエラーとなります。そこで,以下のようにコマンドラインでコンパイル・実行します。
Jupyter Notebook 環境そのものは使わずに,弘大 JupyterHub での fortran プログラムを編集・コンパイル・実行する方法については,以下のページをご覧ください。
$ gfortran -o engelread engelread.f95
$ ./engelread
消費支出?
110000
飲食費?
37000
エンゲル係数は 33.6363640
(これまでの例では1回しか計算していませんが)よく使う計算を「関数」として定義する例です。
以下の例では,エンゲル係数を計算する関数 calc()
を定義し,主プログラムでその関数を呼び出しています。
calc
が実数型であることは,主プログラム内でも,関数副プログラム内でも real
として宣言する必要があります。
program engelcalc
implicit none
real :: spending, food, eng, calc
spending = 110000
food = 37000
eng = calc(spending, food)
print *, "エンゲル係数は", eng
end program engelcalc
function calc(spending, food)
! 消費支出 spending と食費 food から
! エンゲル係数を計算
implicit none
real :: calc, spending, food
calc = food * 100 / spending
end function calc
消費支出を入力すると,エンゲル係数が 20 ~ 40 となる飲食費を表示するプログラムをつくりなさい。
関数は上記の例のように自分で定義することもできますが,あらかじめ用意されている「組み込み関数」を使うこともできます。
以下では,組み込み関数の例として,
floor()
nint()
を使用しています。ちなみに,
ceiling()
sqrt(x)
は平方根,sin(x)
, cos(x)
, tan(x)
exp(x)
, 自然対数は log(x)
, 常用対数は log10(x)
int(x)
は実数を整数に変換し,real(n)
は整数を実数に変換します。 dble(x)
は倍精度実数に変換するのでした。program func
implicit none
real :: spending, food, eng
spending = 110000
food = 37000
eng = food * 100 / spending
print *, 'エンゲル係数は', eng
print *, '切り捨てでは ', floor(eng)
print *, '四捨五入では ', nint(eng)
end
Fortran で使える組み込み関数については,インターネット検索すると,いろいろ見つかります。たとえば,以下の JAMSTEC のサイトが参考になります。
小数点以下を切り上げたエンゲル係数の値を表示させるように上記のプログラムを変更しなさい。
program engelif
implicit none
real :: spending, food, eng
print *, "消費支出? "
read *, spending
print *, "飲食費? "
read *, food
eng = food * 100 / spending
if (eng >= 40) then
print *, "エンゲル係数は", eng, "です。高い値です。"
else
print *, "エンゲル係数は", eng, "です。"
end if
end program engelif
Fortran kernel は Jupyter Notebook 環境で Fortran を実行できるので大変便利なのですが,残念なことに read
文には対応していないため,上記のようなエラーとなります。そこで,以下のようにコマンドラインでコンパイル・実行します。
Jupyter Notebook 環境そのものは使わずに,弘大 JupyterHub での fortran プログラムを編集・コンパイル・実行する方法については,以下のページをご覧ください。
$ gfortran -o engelif engelif.f95
$ ./engelif
消費支出?
110000
飲食費?
43000
エンゲル係数は 39.0909081 です。
$ ./engelif
消費支出?
110000
飲食費?
45000
エンゲル係数は 40.9090919 です。高い値です。
条件分岐の if
文の書式は以下のようになっています。
if (条件式1) then
条件式1 が満たされた場合に実行する文
else if (条件式2)
条件式2 が満たされた場合に実行する文
else if (条件式3)
条件式3 が満たされた場合に実行する文
...
else
それ以外の場合に実行する文
end if
上の例では,eng >= 40
つまり eng
が 40
以上のとき,という関係演算子を使いました。他の例は以下のとおりです。
a == b
a
と b
が等しいとき
a /= b
a
と b
が等しくないとき
a < b
a
が b
より小さいとき
a <= b
a
が b
以下のとき
a > b
a
が b
より大きいとき
a >= b
a
が b
以上のとき
以下の例では,1行目で代入された消費支出に対してエンゲル係数を計算します。
1回だけ計算して表示するのではなく,飲食費の初期値(ここでは 30000
円)から,エンゲル係数を計算して表示し,条件式(ここでは 50000
円以下)が成り立つ場合に,5000
円ずつ増やしながら計算と表示を続けます。
表示の際に,実数を整数に変換する関数 int
や小数点以下を四捨五入する関数 nint
を使っています。
program engelloop1
implicit none
real :: spending, food, eng
spending = 110000
food = 30000
print *, "消費支出", int(spending), "円に対するエンゲル係数の値"
do while (food <= 50000)
eng = food * 100 / spending
print *, " 飲食費", int(food), "円の場合: ", nint(eng)
food = food + 5000
end do
end program engelloop1
do while
文の書式は以下の通りです。
do while (条件式)
条件式が成り立つ場合に実行される文
end do
この繰り返し処理の do while
文の中に出てくる
food = food + 5000
という書き方は,最初はびっくりするかも知れません。
数学の式では,両辺が等しくなることはないので意味のない式ですが,Fortran では,=
は左辺と右辺が等しいことを示す「等号」ではなく,右辺で計算した数値を左辺に「代入」することを表します。したがってこの式は,
現在の変数
food
の値に5000
を足した値をあらためて変数food
に代入する
という意味になります。
同様の内容を do
文で書いた例が以下のプログラムです。
program engelloop2
implicit none
integer :: i
real :: spending, food, eng
spending = 110000
print *, "消費支出", int(spending), "円に対するエンゲル係数の値"
do i = 30000, 50000, 5000
food = float(i)
eng = food * 100 / spending
print *, " 飲食費", int(food), "円の場合: ", nint(eng)
end do
end program engelloop2
do
文の書式は以下の通りです。
do ivar = istart, iend, istep
ivar が istart から iend の間,実行される文
end do
ivar
は整数変数で,その初期値は istart
。文が実行されるごとに,ivar
の値は istep
だけ加算され,iend
の値になるとループが終了します。istep
を省略した場合は 1
とみなします。
上記の例では,念のため,基本整数型(integer
)の i
を実数型(real
)に変換する組み込み関数 float()
を使っています。
以下の値を求めるプログラムを作成せよ。
$$\sum_{n = 1}^{10} n^2$$成分が5個の配列を宣言し,各成分に数値を入れたり表示したりする例です。
program array
implicit none
integer :: a(5), b(5), i
! a の各成分を設定
a = [5, 4, 3, 2, 1]
print *, '全成分を一挙に表示'
print *, a
print *, '縦に表示'
do i = 1, 5
print *, a(i)
end do
! b の全成分を 0 に
b = 0
! 第3成分の値を設定
b(3) = 3
print *, '1つだけ 0 でない成分をもつ配列'
print *, b
end program array
以下のようなテキストファイル mydata.txt
を用意します。
$ cat mydata.txt
1
2
3
4
5
6
以下のプログラムでは,ファイル mydata.txt
の数値を成分が6個の配列 a
に読み込み,成分の総和,平均,最大値,最小値を出力します。
program readdata1
implicit none
integer :: row
real :: a(6)
open(10, file = 'mydata.txt')
read(10, *) a
do row = 1, 6
print *, a(row)
end do
! 成分の和は sum(a)
print *, '成分の総和は', sum(a)
! 成分の個数は size(a)
print *, '成分の平均は', sum(a)/size(a)
print *, '最大値は ', maxval(a)
print *, '最小値は ', minval(a)
end program readdata1
次は,以下のような「6行2列」のデータが入ったテキストファイル mydata2.txt
を読み込む例です。
$ cat mydata2.txt
1 1
2 4
3 9
4 16
5 25
6 36
program readdata2ng
implicit none
integer :: row, col
integer :: b(6, 2)
open(10, file = 'mydata2.txt')
read(10, *) b
do row = 1, 6
print *, (b(row, col), col = 1, 2)
end do
end program readdata2ng
上の結果からわかるように,dimension(6, 2)
の配列 b
を用意して read(10, *) b
とするだけでは,期待したようには読み込まれません。
ちょっとトリッキーですが,dimension(2, 6)
の配列 に読み込ませて,その転置をとってみました。以下の例を参照。
program readdata2
implicit none
integer :: row, col
integer :: b(2, 6), c(6, 2)
open(10, file = 'mydata2.txt')
! あえて dimension(2, 6) の配列 b に読み込む
read(10, *) b
! b の転置行列を c に代入
c = transpose(b)
do row = 1, 6
print *, (c(row, col), col = 1, 2)
end do
end program readdata2
ファイルの書き込み例です。この例では,$x$ と $\sin x$ の値を $x = 0.1$ から $x = 1.0$ まで書き込んでいます。
program writedata
implicit none
integer :: i
real :: x
open(12, file = 'myoutput.txt')
do i = 1, 10
x = i * 0.1
write(12, *) x, sin(x)
end do
end program writedata
$ cat myoutput.txt
0.100000001 9.98334214E-02
0.200000003 0.198669329
0.300000012 0.295520216
0.400000006 0.389418334
0.500000000 0.479425550
0.600000024 0.564642489
0.699999988 0.644217670
0.800000012 0.717356086
0.900000036 0.783326924
1.00000000 0.841470957
画面への出力の際,print *, ...
などと書いていたが,この *
の部分について。
*
は(特に書式をしないで)おまかせで出力させる。書式を設定したいときには *
の場所に以下のように書く。
program formatf
implicit none
character(10) :: ten = '1234567890'
! 15桁(以上)の有効桁数を持つ実数の定義
integer, parameter :: dp = selected_real_kind(p=15)
real(dp) :: sqrtwo
sqrtwo = sqrt(2.0E0_dp)
print '(2a10, a)', ten, ten, ' 桁数'
print '(f8.5)', sqrtwo
print '(f18.15)', sqrtwo
! 入れ子にする場合は" と ' で使い分ける
print '("2の平方根は", f18.15, " です。")', sqrtwo
end program formatf
以下の出力を,エンゲル係数の値を(四捨五入して)小数点以下1桁まで表示するように変更しなさい。
program engel2
implicit none
real :: spending, food, eng
spending = 110000
food = 37000
eng = food * 100 / spending
print *, 'エンゲル係数は ', eng, 'です。'
end program engel2
非常に大きい数や小さい数を表示するときには, $ 1.9891 \times 10^{30}$ や $6.67430 \times 10^{-11}$ などのような指数表示をします。
例えば '(e9.2)'
では,全体で 9
桁の表示分を確保し,小数点以下は 2
桁表示します。
このままだと整数部分が 0
で始まるので,'(1pe9.2)'
としてみます。
program formate
implicit none
character(10) :: ten
! 光速は9桁の定義定数(誤差無し)
integer, parameter :: c = 299792458
! 9桁(以上)の有効桁数を持つ実数の定義
integer, parameter :: dp = selected_real_kind(p=9)
real(dp) :: val
val = c
ten = '1234567890'
print '(2a10, a)', ten, ten, ' 桁数'
print '(e9.2)', val
print '(1pe9.2)', val
print '(1pe15.8)', val
! 入れ子にする場合は" と ' で使い分ける
print '("光速は", 1pe15.8, " m/s です。")', val
end program formate
文字型の書式指定は '(a)'
です。'(a5)'
のように表示する文字数を(この場合は 5
に)指定することもできます。'(2a5)'
は '(a5, a5)'
のことです。
program formata
implicit none
character(10) :: ten
ten = '1234567890'
print '(a, a)', ten, ten
print '(2a)', ten, ten
print '(2a5)', ten, ten
print '(a6,a4)', ten, ten
end program formata
整数の書式指定は '(i5)'
のように i
と桁数(今の場合は 5
)をつけて行います。
program formati
implicit none
integer :: v(4), m(4, 4)
integer :: row, col
! v の成分
v = [1, 23, 456, 7890]
do row = 1, 4
print '(i5)', v(row)
end do
! 4行4列の成分を全て 0 に
m = 0
! 対角成分
m(1,1) = 1
m(2,2) = 23
m(3,3) = 456
m(4,4) = 7890
! 書式指定で小綺麗に表示
do row = 1, 4
print '(4i5)', (m(row, col), col = 1, 4)
end do
end program formati