現在(&lastmod;)作成中です。 既に書いている内容も&color(#ff0000){変わる};可能性が高いので、注意。 ------- #contents ------- * はじめに [#sec0c77d] ** 特徴 [#s2563d6c] C言語との差分に注目したFortran90/95言語の入門 ** 目標 [#k7f271a6] この「計算科学演習I」で扱うFotranソースコードが自由に読めること。 Fortran90/95の特性を活かした綺麗なコードが自由に書けること。 * エディタとコンパイラ [#e4275a16] - Emacsのmajorモードをf90に設定すると便利(f95モードは用意されていない) -- ミニバッファでf90-modeと打つ (Esc-x f90-mode) -- tabキーでソースコード(現在行)を整形。 -- tabキーで適切なend構文を自動入力。 -- ^jで整形改行 - scalar上でのコンパイルコマンドは pgf95 またはpgf90 - フリーのFortran95コンパイラ -- gfortran (gnu Fortran): http://gcc.gnu.org/fortran/ -- g95: http://www.g95.org/ * Fortran90/95の歴史 [#mbcaa859] - ''FORTRAN''(全部大文字の)という言語は存在しない。 --''Fortran90''や''Fortran95''という言語ならある。 - FORTRAN66。 --1966年に標準化。 - FORTRAN77 --1977年に標準化。 -- if/then/else -- 広まった - Fortran90 -- 1991年に標準化。 -- 大幅な改訂 - Fortran95 --F90からのマイナーバージョンアップ Fortran90はFORTRAN77とは大きく違う。違う言語と考えるべき。 ~| f95 - f90 | << | f90 - f77 | * なぜFortran90/95か? [#l4a65bd7] + 計算速度が速い -- スーパーコンピューティングは速さが命。 -- C/C++でもFortranと同じくらい速いコードは書けるが、遅いコードも書けてしまう。 --- 言語としての自由度が高いために、コンパイラが困る。最適化ができなくなる。 + 便利 -- 道具(言語)は目的にあったものを -- 数学的計算にはFortran90/Fotran95が適している --- '''For'''-mula '''Tran'''-slator -- 数値計算ライブラリの抱負な蓄積。 * Fortran Legacy [#m1ff3a26] ** LegacyなコードはFortranの財産 [#b1d8f826] ** Legacyな人間は負の遺産。 [#a208ae5d] 『FORTRAN77に何の不満もない。これでxx年やってきた。』 ・・・時代に33年(2010-1977)遅れている。 * 定番、hello, world プログラム [#r780ea81] K&Rに出てくるC言語でのhello, worldプログラムは、 #include <stdio.h> main () { printf("hello, world.\n"); } である。 Fortran90/95ではこうなる。 program hello_world implicit none print *, "hello, world." end program hello_world このようにほとんど同じである。 いくつかの細かい違いは、 - Fortran90/95ではmainプログラムは program '''name'''として宣言する。'''name'''は任意。 - implicit noneの意味は後述。 - Fortran90/95では標準入出力機能は言語に組み込まれている(ヘッダファイルをインクルードする必要はない(そもそもヘッダファイルというものがない))。 - Fortran90/95では print *で標準出力へ。 - 行末にセミコロン不要。 - 文字列はダブルクォーテーションマーク(")で囲む。(C言語と同じ) - シングルクウォーテーションマーク(')でも同じ意味。(これはC言語と違う) * &color(#0000ff){【演習】}; [#bf1c1111] 計算機「scalar」 で、上のFortran90/95プログラムをエディタで入力し、 ファイル名hello_world.f95として保存せよ。 - ls -l コマンドでファイルを確認せよ。 - hello_world.f95 をコンパイルせよ。 -- pgf95 hello_world.f95 - 実行せよ。 -- ./a.out そのプログラム(hello_world.f95)のメッセージ("hello, world.")を自由に変え、 コンパイル&実行せよ。 * サンプルコードによるFortran90/95入門 [#t53ce9d9] まずはC言語の復習も兼ねて、線形合同法による乱数生成プログラムのC言語版を見てみよう。 ** C言語版 [#q5b62857] /* * random number by * the linear congruential generator (lcg): * * n_new = (n_old*A+C) % M * */ #include <stdio.h> #define A 109 #define C 1021 #define M 32768 /* = 2^15 */ int lcg_n = 10; /* initial value */ double random_0_to_1(void) { lcg_n = (lcg_n*A+C) % M; return (double)lcg_n / (M-1); } main(void) { int i; int loop_max = 100; int n = 10; // initial condition for (i=0; i<loop_max; i++) { printf("%d %f\n",i, random_0_to_1()); } } ** Fortran90/95版 [#a4f54c2b] 次にFortran90/95言語で同じことをするコードを見てみよう。 ! ! random number by ! the linear congruential generator (lcg): ! ! n_new = mod(n_old*A+C,M) ! ! by Akira Kageyama, Kobe Univ. ! on 2010.05.05. ! module lcg implicit none integer, parameter :: SP = kind(1.0) integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP)) integer, parameter :: A = 109 integer, parameter :: C = 1021 integer, parameter :: M = 2**15 integer :: lcg_n = 10 ! initial value contains function random_0_to_1() real(DP) :: random_0_to_1 lcg_n = mod(lcg_n*A+C,M) random_0_to_1 = real(lcg_n,DP) / (M-1) end function random_0_to_1 end module lcg program random_number use lcg implicit none integer :: i integer :: loop_max = 100 do i = 1 , loop_max print *,i, random_0_to_1() end do end program random_number moduleはデータと関数をひとまとめにしたもの。 C++やJavaのclassに対応。 * &color(#0000ff){【演習】}; [#ff61e3dd] 上のFortran90/95プログラムをrandom_number.f95という名前のファイルに保存し、 コンパイル&実行せよ。 * Fortran90/95に対する誤解 [#nca1bce0] 以下全てFORTRAN66/77と誤解している。 - 変数名は6文字までなんでしょう?・・・そんなことはありません。 - ソースコードは全部大文字で書くんでしょう?・・・それも可能ですが、わざわざそんな読みにくいことをする必要はありません。 - ソースコードは固定形式(7列目から72列目まで)なんでしょう?・・・そんなことはありません。 - 構造体がないんでしょう?・・・あります。よく使います。 - ポインタがないんでしょう?・・・あります。あまり使う必要はありませんが。 - 演算子を自分で定義することができなんでしょう?・・・できます。よく使います。 - 関数や演算子の多重定義ができなんでしょう?・・・できます。よく使います。 - 再帰呼び出しができなんでしょう?・・・できます。あまり使いませんが。 - 関数とデータをまとめてひとかたまりにする(クラス化する)ことなんてできないんでしょう?・・・できます。よく使います。上のサンプルプログラムで使ったmoduleがそれです。 - 情報の隠蔽(カプセル化)ができないんでしょう?・・・できます。よく使います。 * 実際的なサンプルコード [#pcddf2ad] 上の最後に述べた情報隠蔽(カプセル化)は名前空間を分離し、 バグの入りにくいコードを書くためには不可欠な機能である。 上のサンプルプログラムにカプセル化を施したコードを下に示す。 冒頭近くに3行ほど追加しただけである。 諸君は既にJavaを知っているから、この3行の意味は明白であろう。 ! ! random number by ! the linear congruential generator (lcg): ! ! n_new = mod(n_old*A+C,M) ! ! by Akira Kageyama, Kobe Univ. ! on 2010.05.05. ! module lcg implicit none private ! default public :: SP, DP ! constants public :: random_0_to_1 ! function integer, parameter :: SP = kind(1.0) integer, parameter :: DP = selected_real_kind(2*precision(1.0_SP)) integer, parameter :: A = 109 integer, parameter :: C = 1021 integer, parameter :: M = 2**15 integer :: lcg_n = 10 ! initial value contains function random_0_to_1() real(DP) :: random_0_to_1 lcg_n = mod(lcg_n*A+C,M) random_0_to_1 = real(lcg_n,DP) / (M-1) end function random_0_to_1 end module lcg program random_number_capsule use lcg implicit none integer :: i integer :: loop_max = 100 do i = 1 , loop_max print *,i, random_0_to_1() end do end program random_number_capsule * Fortran90/95の特徴 [#y74761e3] 数値演算や計算機シミュレーションに適した言語である。 特にスーパーコンピュータ向けの言語としては最適。 オブジェクト指向プログラミングや並列計算機への対応も意識された言語である。 - これらの意味を納得できるようになることが目標。 * C言語と極端に差がつく例 [#gbfec419] 部屋の中の温度場の分布から平均気温を求めよう。 温度場を3次元float配列 f(nx,ny,nz)で表す。3つの整数nx, ny, nzは実行時まで不定とする。 平均気温を求めるには、 全ての格子点(i,j,k)上でのfの値を足して格子点の総数で割ればよいものとする。 任意サイズの3次元単精度実数(浮動小数点数)配列を受け取り、 その平均値を返す関数を作ろう。 ** Fortran90/95ではこう書ける。 [#f2c59871] このようにわずか4行で書ける。 real function mean_value(f) real, dimension(:,:,:) :: f mean_value = sum(f) / (size(f,1)*size(f,2)*size(f,3)) end function mean_value このプログラムの意味は後で説明するが、 ここに出てくるsumやsizeはライブラリではなく組み込み関数、 つまりFortran90/95言語にもともと含まれている関数である。 ** &color(#0000ff){【演習】}; [#eab3a00a] 上の関数をC言語で作れ。 - scalar上でのCコンパイラはcc, gccである。 * サンプルプログラム:複素数 [#gba646e1] //$e^{i\pi} = -1$ //つまり #ref(e_i_pi.jpg) をFortran95で書くと、 complex :: i = (0.0,1.0) real :: pi = 3.141593 print *,' exp(i*pi) = ', exp(i*pi) * 行列計算 [#e6951de8] 行列のかけ算のためにFortran90/95ではmatmulという組み込み関数が用意されている。 また、行列の転置をとる組み込み関数transposeも用意されている。 したがって、行列A(例えば10行10列の2次元配列)の転置と別の行列Bの積を計算しそれを行列Cとする計算、つまり // C = {}^t\!A\, B #ref(matmul_a_transpose_b.jpg) をFortran90/95で書くと、 real, dimension(10,10) :: A, B, C C = matmul(transpose(A),B) と一行で書ける。 行列の足し算も簡単である。 上の式の右辺に別の行列Dを加える演算はそのまま C = matmul(transpose(A),B) + D と書けばよい。 * &color(#0000ff){【演習】}; [#hc79163b] 以下のプログラムをmatmul_a_transpose_b.f95というファイルに保存し、 コンパイル&実行せよ。 program matmul_a_transpose_b implicit none real, dimension(10,10) :: A, B, C, D A = 1.0 ! A(i,j) = 1 for all i and j. B = 2.0 ! B(i,j) = 2 D = 3.0 ! D(i,j) = 3 C = matmul(transpose(A),B) + D print *, 'max element of C is ', maxval(C) end program matmul_a_transpose_b * サンプルプログラム: 級数 [#zbbbb62a] 級数 #ref(series_one_forth.jpg) // \sum_{n=1}^\infty\,\frac{1}{i}\cdot\frac{1}{i+1}\cdot\frac{1}{i+2} = \frac{1}{1}\cdot\frac{1}{2}\cdot\frac{1}{3} + \frac{1}{2}\cdot\frac{1}{3}\cdot\frac{1}{4} + \frac{1}{3}\cdot\frac{1}{4}\cdot\frac{1}{5} + \cdots = \frac{1}{4} の最初の1000項の和を求めるプログラムも、式をそのまま書けばいい。まさにFormula Translation。 program series_one_forth implicit none integer, parameter :: nterms = 1000 real, dimension(nterms) :: x, y, z integer :: i do i = 1 , nterms x(i) = 1.0 / i y(i) = 1.0 / (i+1) z(i) = 1.0 / (i+2) end do print *,'ans = ', sum(x*y*z) end program series_one_forth * &color(#0000ff){【演習】}; [#f13dd55b] 上のプログラムを series_one_forth.f95という名前に保存し、 コンパイル&実行せよ。 * 部分配列 [#y64cdf79] 上の例では第1行目から第1000項目までの和をとっていた。 もしも第50項目から第90行目までの和が欲しければ以下のようにすればよい。 program series_one_forth_02 implicit none integer, parameter :: nterms = 1000 real, dimension(nterms) :: x, y, z integer :: i do i = 1 , nterms x(i) = 1.0 / i y(i) = 1.0 / (i+1) z(i) = 1.0 / (i+2) end do print *,'ans02 = ', sum(x(50:90)*y(50:90)*z(50:90)) end program series_one_forth_02 ここに出てきた x(50:90)などは部分配列と呼ぶ。 部分配列も含めてFortran90/95には配列を扱うための便利な機能が多数用意されている。 配列機能については次回詳しく解説する。 * サンプルプログラム: 3次元ベクトル場のエネルギー積分 [#oc44c8ef] 空間中に分布する磁場(3成分ベクトル場)の全磁気エネルギーは #ref(emag01.jpg) で与えられる。 計算機シミュレーションコードでこの磁場が 3次元配列 Bx(nx,ny,nz), By(nx,ny,nz), Bz(nx,ny,nz) で書かれているとき、 全磁気エネルギーは #ref(emag02.jpg) と書ける。この量を計算して出力するFortran90/95コードは、 定数倍のファクターの違いは無視して print *,' energy = ', sum(Bx**2+By**2+Bz**2) とわずか1行で書ける。 配列のサイズ(nx,ny,nz)が実行時まで不定でもよい。 * ソースコードの書き方 [#gc9544f5] ** 大文字・小文字 [#y679c41a] - Fortran90/95のソースコードではアルファベットの大文字と小文字は区別しない。 -- kobe と Kobe と KOBE は同じ -- 文字列変数の中では当然区別される。 ** 自由形式。C言語とほぼ一緒。 [#e45237a0] - 1行は132文字まで * 予約語 [#f9a88182] Fortran90/95には予約語が存在しない。 ユーザ定義の変数名であればコンパイラが賢く判断してくれる。 一方、C言語には32個も予約語がある。 C99ではさらに5個増えた。 たとえば・・・ 次のCプログラムはコンパイルエラーになる。 main() { int dee, daa, doo; int de, da, do; dee = de*de; } * &color(#0000ff){【演習】}; [#wcc364ec] 下のプログラムをreserved_words.f95というファイルに保存し、 予約語がないことを確認せよ。 program reserved_words implicit none integer :: id=2, ie=3 integer :: if if ( id==ie ) if=0 end program reserved_words * コメント行の書き方 [#k20a31b6] C言語では /* この間がコメント */ である。(C99では // から行末までもコメントとなった。) Fortran90/95では ! 一行の中でこの文字以降がコメント( C99やC++、JAVAの//と同じ) * &color(#0000ff){【演習】}; [#i0629ed9] これまでに作ったプログラムに自由にコメントを入れよ。 * 定数 [#z3e31dfd] C言語ではプリプロセッサを使って #define NX 100 とするが、 Fortran90/95ではC++のconstと同様 integer, parameter :: NX = 100 と書ける。コンパイラが型チェックをしてくれるのでバグが入りにくい。 * &color(#0000ff){【演習】}; [#u0ded96b] 上で書いたseries_one_forth.f95の定数ntermsを設定(integer, parameter nterms=100)した後に、 変更(nterms=200など)するとエラーになることを確認せよ。 * 数値演算 [#k4190895] 四則演算はC言語と同じ a+b, a-b, a*b, a/b 優先順位も同じ。 Fortran90/95でのべき乗は a**b 余りは(C言語ではa % b)はFortran90/95では、 mod(a,b) * &color(#0000ff){【演習】}; [#be184ff7] 以下のプログラムをelemental_operators.f95に保存して、 コンパイル&実行せよ。 program elemental_operators implicit none real, parameter :: pi = 3.141593 complex :: z print *, ' pi = ', pi print *, ' pi+pi = ', pi+pi print *, ' pi-pi = ', pi-pi print *, ' pi*pi = ', pi*pi print *, ' pi/pi = ', pi/pi print *, 'pi**pi = ', pi**pi z = (pi,pi) print *,' z = ', z print *,' z*pi = ', z**pi print *,' z**z = ', z**z end program elemental_operators ------------------------------ as of &_now; (&counter;)