現在(2011-04-21 (木) 10:28:29)作成中です。 既に書いている内容も変わる可能性が高いので、注意。



はじめに

特徴

C言語との差分に注目したFortran90/95言語の入門

目標

この「計算科学演習I」で扱うFotranソースコードが自由に読めること。 Fortran90/95の特性を活かした綺麗なコードが自由に書けること。

エディタとコンパイラ

  • Emacsのmajorモードをf90に設定すると便利(f95モードは用意されていない)
    • ミニバッファでf90-modeと打つ (Esc-x f90-mode)
    • tabキーでソースコード(現在行)を整形。
    • tabキーで適切なend構文を自動入力。
    • ^jで整形改行
  • scalar上でのコンパイルコマンドは pgf95 またはpgf90

Fortran90/95の歴史

  • FORTRAN(全部大文字の)という言語は存在しない。
    • Fortran90Fortran95という言語ならある。
  • FORTRAN66。
    • 1966年に標準化。
  • FORTRAN77
    • 1977年に標準化。
    • if/then/else
    • 広まった
  • Fortran90
    • 1991年に標準化。
    • 大幅な改訂
  • Fortran95
    • F90からのマイナーバージョンアップ

Fortran90はFORTRAN77とは大きく違う。違う言語と考えるべき。

| f95 - f90 | << | f90 - f77 |

なぜFortran90/95か?

  1. 計算速度が速い
    • スーパーコンピューティングは速さが命。
    • C/C++でもFortranと同じくらい速いコードは書けるが、遅いコードも書けてしまう。
      • 言語としての自由度が高いために、コンパイラが困る。最適化ができなくなる。
  2. 便利
    • 道具(言語)は目的にあったものを
    • 数学的計算にはFortran90/Fotran95が適している
      • For-mula Tran-slator
    • 数値計算ライブラリの抱負な蓄積。

Fortran Legacy

LegacyなコードはFortranの財産

Legacyな人間は負の遺産。

『FORTRAN77に何の不満もない。これでxx年やってきた。』

・・・時代に33年(2010-1977)遅れている。

定番、hello, world プログラム

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言語と違う)

【演習】

計算機「scalar」 で、上のFortran90/95プログラムをエディタで入力し、 ファイル名hello_world.f95として保存せよ。

  • ls -l コマンドでファイルを確認せよ。
  • hello_world.f95 をコンパイルせよ。
    • pgf95 hello_world.f95
  • 実行せよ。
    • ./a.out

そのプログラム(hello_world.f95)のメッセージ("hello, world.")を自由に変え、 コンパイル&実行せよ。

サンプルコードによるFortran90/95入門

まずはC言語の復習も兼ねて、線形合同法による乱数生成プログラムのC言語版を見てみよう。

C言語版

/*
 * 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版

次に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に対応。

【演習】

上のFortran90/95プログラムをrandom_number.f95という名前のファイルに保存し、 コンパイル&実行せよ。

Fortran90/95に対する誤解

以下全てFORTRAN66/77と誤解している。

  • 変数名は6文字までなんでしょう?・・・そんなことはありません。
  • ソースコードは全部大文字で書くんでしょう?・・・それも可能ですが、わざわざそんな読みにくいことをする必要はありません。
  • ソースコードは固定形式(7列目から72列目まで)なんでしょう?・・・そんなことはありません。
  • 構造体がないんでしょう?・・・あります。よく使います。
  • ポインタがないんでしょう?・・・あります。あまり使う必要はありませんが。
  • 演算子を自分で定義することができなんでしょう?・・・できます。よく使います。
  • 関数や演算子の多重定義ができなんでしょう?・・・できます。よく使います。
  • 再帰呼び出しができなんでしょう?・・・できます。あまり使いませんが。
  • 関数とデータをまとめてひとかたまりにする(クラス化する)ことなんてできないんでしょう?・・・できます。よく使います。上のサンプルプログラムで使ったmoduleがそれです。
  • 情報の隠蔽(カプセル化)ができないんでしょう?・・・できます。よく使います。

実際的なサンプルコード

上の最後に述べた情報隠蔽(カプセル化)は名前空間を分離し、 バグの入りにくいコードを書くためには不可欠な機能である。 上のサンプルプログラムにカプセル化を施したコードを下に示す。 冒頭近くに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の特徴

数値演算や計算機シミュレーションに適した言語である。 特にスーパーコンピュータ向けの言語としては最適。 オブジェクト指向プログラミングや並列計算機への対応も意識された言語である。

  • これらの意味を納得できるようになることが目標。

C言語と極端に差がつく例

部屋の中の温度場の分布から平均気温を求めよう。 温度場を3次元float配列 f(nx,ny,nz)で表す。3つの整数nx, ny, nzは実行時まで不定とする。 平均気温を求めるには、 全ての格子点(i,j,k)上でのfの値を足して格子点の総数で割ればよいものとする。

任意サイズの3次元単精度実数(浮動小数点数)配列を受け取り、 その平均値を返す関数を作ろう。

Fortran90/95ではこう書ける。

このようにわずか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言語にもともと含まれている関数である。

【演習】

上の関数をC言語で作れ。

  • scalar上でのCコンパイラはcc, gccである。

サンプルプログラム:複素数

e_i_pi.jpg

をFortran95で書くと、

complex :: i = (0.0,1.0)
real :: pi = 3.141593
print *,' exp(i*pi) = ', exp(i*pi)

行列計算

行列のかけ算のためにFortran90/95ではmatmulという組み込み関数が用意されている。 また、行列の転置をとる組み込み関数transposeも用意されている。 したがって、行列A(例えば10行10列の2次元配列)の転置と別の行列Bの積を計算しそれを行列Cとする計算、つまり

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

と書けばよい。

【演習】

以下のプログラムを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

サンプルプログラム: 級数

級数

series_one_forth.jpg

の最初の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

【演習】

上のプログラムを series_one_forth.f95という名前に保存し、 コンパイル&実行せよ。

部分配列

上の例では第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次元ベクトル場のエネルギー積分

空間中に分布する磁場(3成分ベクトル場)の全磁気エネルギーは

emag01.jpg

で与えられる。 計算機シミュレーションコードでこの磁場が 3次元配列 Bx(nx,ny,nz), By(nx,ny,nz), Bz(nx,ny,nz) で書かれているとき、 全磁気エネルギーは

emag02.jpg

と書ける。この量を計算して出力するFortran90/95コードは、 定数倍のファクターの違いは無視して

print *,' energy = ', sum(Bx**2+By**2+Bz**2)

とわずか1行で書ける。 配列のサイズ(nx,ny,nz)が実行時まで不定でもよい。

ソースコードの書き方

大文字・小文字

  • Fortran90/95のソースコードではアルファベットの大文字と小文字は区別しない。
    • kobe と Kobe と KOBE は同じ
    • 文字列変数の中では当然区別される。

自由形式。C言語とほぼ一緒。

  • 1行は132文字まで

予約語

Fortran90/95には予約語が存在しない。 ユーザ定義の変数名であればコンパイラが賢く判断してくれる。 一方、C言語には32個も予約語がある。 C99ではさらに5個増えた。

たとえば・・・ 次のCプログラムはコンパイルエラーになる。

main()
{
  int dee, daa, doo;
  int de,  da,  do;
  dee = de*de;
}

【演習】

下のプログラムをreserved_words.f95というファイルに保存し、 予約語がないことを確認せよ。

program reserved_words
  implicit none
  integer :: id=2, ie=3
  integer :: if
  if ( id==ie ) if=0
end program reserved_words

コメント行の書き方

C言語では

 /* この間がコメント */

である。(C99では // から行末までもコメントとなった。)

Fortran90/95では

  ! 一行の中でこの文字以降がコメント( C99やC++、JAVAの//と同じ)

【演習】

これまでに作ったプログラムに自由にコメントを入れよ。

定数

C言語ではプリプロセッサを使って

#define NX 100

とするが、

Fortran90/95ではC++のconstと同様

integer, parameter :: NX = 100

と書ける。コンパイラが型チェックをしてくれるのでバグが入りにくい。

【演習】

上で書いたseries_one_forth.f95の定数ntermsを設定(integer, parameter nterms=100)した後に、 変更(nterms=200など)するとエラーになることを確認せよ。

数値演算

四則演算はC言語と同じ

a+b, a-b, a*b, a/b

優先順位も同じ。

Fortran90/95でのべき乗は

a**b

余りは(C言語ではa % b)はFortran90/95では、

mod(a,b)

【演習】

以下のプログラムを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 2019-08-20 (火) 10:14:18 (1646)