現在(2018-04-05 (木) 13:09:14)作成中です。 既に書いている内容も大幅に変わる可能性が高いので、注意。 as of 2019-09-23 (月) 06:08:52 (222)



担当教員

陰山 聡

演習日

  • 2010.05.06
  • 2010.05.13

概要と達成目標

  1. 概要
    1. why f? (なぜFortran95言語か)
    2. f && c (Fotran95とCの同じところ)
    3. [演習]Hello, world. コンパイルと実行
    4. diff f c (Fotran95とCの違うところ)
    5. wonderful f (Fotran95の素晴らしいところ)
    6. practical f (Fotran95によるコーディングの実際)
  2. 目標
    1. この「計算科学演習I」で扱うFotranソースコードが自由に読めること。
    2. 2次元拡散方程式を解くコードをFotran95言語の特性を活かして書けるようになること。

はじめに

エディタとコンパイラ

  • 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
    • 数値計算ライブラリの抱負な蓄積。Legacyなコード。財産。

間違った理由

『FORTRAN77に何の不満もない。これでxx年やってきた。 (新しい言語を勉強するのは面倒だ。)』

・・・Legacyなコードはいいが、Legacyな人間は良くない。 時代に30年以上(2019-09-23-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

このようにCもFotranも手続き型のプログラミング言語なので、 基本的には同じプログラムになる。 いくつかの細かい違いは、

  • 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に対する誤解

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

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

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

まだ納得できないない?では例を一つ

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

【演習】

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

  • scalar上でのCコンパイラ=cc, gcc

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

このプログラムの意味は後で説明する。

数式の表現が簡単な例(複素数)

$e^{i\pi} = -1$ つまり

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)

と一行で書ける。

【演習】

以下のプログラムをmatmul_a_transpose_b.f95というファイルに保存し、 コンパイル&実行せよ。

program matmul_a_transpose_b
  implicit none
  real, dimension(10,10) :: A, B, C
  A = 1.0  ! A(i,j) = 1	for all	i and j.
  B = 2.0  ! B(i,j) = 2
  C = matmul(transpose(A),B)
  print *, 'max element of C is ', maxval(C)
end program matmul_a_transpose_b

例(級数)

級数

series_one_forth.jpg

の最初の1000項の和を求めるプログラムも、式をそのまま書けばいい。まさにFormula Translation。

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)

例(3次元ベクトル場のエネルギー)

空間中に分布する磁場(3成分のベクトル場)

Bx(100,100,100),  By(100,100,100),  Bz(100,100,100)

の全磁気エネルギーを計算して出力するには、

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

と一行(定義式そのもの)を書けば良い。

C言語との比較によるFortran90/95入門

ソースコード

大文字・小文字

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

【演習】

次のCプログラムはコンパイルエラーになる。 どこにバグがあるが指摘せよ。

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

予約語

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

コメント行の書き方

C言語では

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

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

Fortran90/95では

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

定数

C言語

#define NX 100

Fortran90/95

integer, parameter :: NX = 100

数値演算

四則演算はC言語と同じ

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

べき乗は

a**b

倍精度浮動小数点数の指定方法

kindパラメータ

計算科学では実数の物理量を数値的に表現するのに単精度ではなく、 倍精度の浮動小数点数を使うのが普通である。 (単精度実数では精度が不十分なため。)

まは、スーパーコンピュータは単精度浮動小数点ではなく、 倍精度浮動小数点数の演算(四則演算)を高速に処理できるよう設計されている。 従ってFortran90/95言語での倍精度浮動小数点数の表現方法に早めに慣れておくことは重要である。

Fortran90/95で倍精度浮動小数点数を表す方法は幾つかある。 最も簡単でportabilityの高い方法は以下の方法であろう。 まず、単精度浮動小数点数の精度に関係したある整数定数を定義する。 ここではその整数をSPという名前にする。 SPとはSingle Precision(単精度)の頭文字をとったものである。

integer, parameter ::  SP = kind(1.0)

次に、このSPを使って単精度として指定した浮動小数点数の2倍の精度を持つ浮動小数点数、 つまり倍精度浮動小数点数に関係したある整数定数(ここではDPという名前)を以下のようにして定義する。

integer, parameter ::  DP = selected_real_kind(2*precision(1.0_SP))

この整数定数DPはDouble Precisionを意味する。

こうして定義した二つの整数SPとDPを使って単精度浮動小数点や倍精度浮動小数点数を宣言する。

real(kind=SP) :: a
real(kind=DP) :: b

とすると、aは単精度浮動小数点、bは倍精度浮動小数点数である。

上のプログラムのkind=の部分は省略可能である。 つまり

real(SP) :: a
real(DP) :: b

としてもよい。

数値データの精度指定にもこのSPやDPを使う。

a = 1.0_SP    ! 単精度の1.0
b = 1.0_DP    ! 倍精度の1.0

サンプルプログラムを見てみよう。

program double_precision_real
  implicit none
  integer, parameter ::  SP = kind(1.0)
  integer, parameter ::  DP = selected_real_kind(2*precision(1.0_SP))
  real(SP) :: a
  real(DP) :: b
  a = 3.1415926535897323846264338327950288_SP
  b = 3.1415926535897323846264338327950288_DP
  print *,' a, b = ', a, b
end program double_precision_real

【演習】

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

【演習】

double_precision_real.f95に現れる二つのkindパラメータSPとDPには実際にはどんな値が入っているか?

数字表現

  • 浮動小数点数
    • 1. ! デフォルト精度
    • 0.1_SP ! kindパラメタSPの浮動小数点数
    • 3.141593_DP ! kindパラメタDPの浮動小数点数
    • 0.31415926535897932e-1_DP ! 上に同じ

型名

C ( C99 )Fortran90/95補註
文字charcharacter
文字列char[n+1]character(len=n)nは文字長
整数intintegerinteger(kind=4)でも可
実数floatrealreal(kind=SP)でも可
倍精度実数doublereal(kind=DP)
「長い」整数longinteger(kind=8でも可)kindの整数はシステム依存
bool( _Bool )logical値は .true. または .false.
複素数( _Complex )complex
構造体structtype詳しくは後述

構造体 (derived type)

C言語の構造体とほとんど同じである。下の例をみよ。

program type
  implicit none

  type student
     character(len=20) :: first_name, last_name
     integer :: age
  end type student

  type(student) :: st

  st = student("Albert", "Einsten", 19)

  print *,st%first_name, st%last_name, st%age
end program type

構造体のメンバにアクセスするのに%を使う。

【演習】

上のプログラムをtype.f90に保存し、コンパイル&実行せよ。

【演習】

student型の構造体変数をもう一つ(例えばst2という名前)を作り、 stのデータをst2にコピーせよ。

宣言

C言語では変数の宣言はブロックの先頭にまとめて置く。(C99やC++はもっと自由だが。) Fortran90/95でも同じ。

暗黙の型宣言

Fortran90/95ではimplicit noneを省略すると「暗黙の型宣言」をしたことになる。 暗黙の型宣言とは、次の6つのアルファベット i,j,k,l,m,n で始まる変数は整数である等のルールである。 バグが入りやすいので、暗黙の型宣言は使わない方が良い。 &color(#ff0000){Fortran90/95プログラムでは常にimplicit none宣言をすること。}

暗黙の型宣言はなぜバグが入りやすいか

宣言したつもりのない変数を間違って使っていても気がつかない可能性がある。 例えば、次のプログラムにはバグがある。すぐに見つけられるか?

program use_implicit_none
  integer :: fresh_meat
  flesh_meat = 100 ! yen
  print *, "today's price = ", fresh_meat
end program use_implicit_none

【演習】

上のプログラムをuse_implicit_nene.f95というファイルに書き込み、コンパイル&実行せよ。

  • バグがあるのにコンパイラは教えてくれないだけでなく、 何事もないかのように正常終了する。しかも、それらしい答えを返す。 バグがあることに気がつかない!

【演習】

今作ったuse_implicit_nene.f95の2行目(program ...の次の行)に implicit noneと書いてから、コンパイルせよ。

行の継続

行末に&をつけると継続行

a = b + c + &
       d  + f

セミコロン

セミコロンは改行と同じ

tmp = right
right = left
left = tmp

tmp = right; right = left; left = tmp

は同じ。

1次元配列

C言語では

array01[NX];

が大きさNXの1次元配列である。 メモリ上に、array01[0], array01[1], array01[2], ..., array01[NX-1]と並んでいる。

Fortran90/95では

integer, dimension(NX) :: array01

と宣言するとメモリ上に、array01(1), array01(2), array01(3), ..., array01(NX)と並ぶ。

配列のインデックスは1から始まる。

【演習】

大きさNXの1次元整数配列array_01を作り、 各要素に1,2,3,...,NXを代入した上で、 全要素の和をとるFortran90/95プログラムを書け。

2次元配列

厳密に言えばC言語では2次元配列はない。 あるのは「配列の配列」である。 つまり二つの整数インデックスiとjを使って

array02[j][i]

という形でアクセスできる量はある。

array02のサイズ(インデックスiとjの範囲)は実行時に不定として、 array02をまるごと引数で受け取る関数がC言語では作れない。 従ってこれを2次元配列とは呼び難い。

メモリ空間中での2次元配列の各要素の位置

上記のarray02は、

array02[0][0], array02[0][1], array02[0][2], ..., array02[0][NX-1]

と並び、そして (それに続く場合もあるし、 まったく別の場所にあるかもしれない)ある位置から、

array02[1][1], array02[1][1], array02[1][2], ..., array02[1][NX-1]

という順番で並ぶ。

一方、 Fortran90/95ではarray02(i,j)という形でアクセスできる量は本物の2次元配列である。 サイズがわかっている場合には

integer, dimension(NX,NY) :: array02

と宣言する。 不明の場合には、

integer, dimension(:,:) :: array02

と宣言し、後ほど(実行時に)メモリをallocateする。 メモリ上では、array02の要素は

array02(1,1), array02(2,1), array02(3,1), ..., array02(NX,1), array02(1,2),...

に全ての要素(NX*NY個)が連続して並ぶ。

配列の要素がメモリ空間中に連続してあることは、 メモリへのアクセス速度が極めて重要となるスーパーコンピューティングでは極めて重要である。

3次元配列

C言語では

array03[k][j][i]

が上の意味での「擬似」3次元配列である。

Fortran90/95では

array03[i][j][k]

が3次元配列で、 宣言は

integer, dimension(NX,NY,NZ) :: array03

または

 integer, dimension(:,:,:) :: array03

とする。

if構文

if文はC言語と似ている。下の例を見れば一目瞭然であろう。

if ('''criterion''') then
  '''action'''
end if

actionが1行で書けるなら

if ( '''criteron''' ) '''action'''

ともできる。

else-if構文も簡単である。

if ('''criterion_1''') then
  '''action_1'''
else if ('''criterion_2''' ) then
  '''action_2'''
else if ('''criterion_3''' ) then
  '''action_3'''
end if

論理式

これも下の例を見ればわかるであろう。

if (a==b) then   ! aとbが等しい
if (a>b) then
if (a>=b) then
if (a<b) then
if (a<=b) then
if (a/=b) then    ! aとbが等しくない

論理値

bool変数はlogicalという型名を持ち、

.true.

.false.

の二つをとる。

logical :: a, b, c
a = .true.
b = .false.
c = a .and. b
c = a .or. b

case 構文

C言語でいうswitch文である。

select  case ('''case expression''')
case ('''case selector 1''')
  '''action 1'''
case ('''case selector 2''')
  '''action 2'''
case ('''case selector 3''')
  '''action 3'''
  .
  .
case default   ! なくても良い。
  '''default action'''
end select

C言語では普通breakが必要だがFortran90/95のcase文には不要である。

【演習】

次のプログラムをcase.f95に保存し、実行せよ。

program case
  implicit none
  integer :: month

  month = 6

  select case (month)
  case (1)
    print *,'January.'
  case (2)
    print *,'February'
  case (3:5)
    print *,'Spring'
  case default
    print *,'Other than spring.'
 end select
end program case

【演習】

後述するが、Fortran90/95は文字列の取り扱いがC言語よりもずっと簡単である。 例えば、長さの違う文字列をselect case構文で簡単に場合分けできる。

program case_string
  implicit none
  character(len=*), parameter :: color = 'white'

  select case (color)
  case ('aka')
    print *, 'red'
  case ('ao')
    print *, 'blue'
  case ('midori')
    print *, 'green'
  case default
    print *, "I don't know the color."
 end select
end program case_string

繰り返し構文(doループ)

C言語でいうfor文である。

incrementが1の場合のdo-loop

program do_loop
  implicit none
  integer :: i
  do i = 1 , 100
     print *, ' i = ', i
  end do
end program do_loop

incrementが1以外の場合のdo-loop

program do_loop
  implicit none
  integer :: i
  do i = 1 , 100 , 2
     print *, ' i = ', i
  end do
end program do_loop

【演習】

  • 上のプログラムをdo_loop.f95に保存し、increment値(上の場合2)を変えて実行せよ。
  • do_loop.f95のdo文を
    do i = 100 , -100 , -2 
    として実行せよ。

exitとcycle

整数のカウンターのないdo-loopも可能である。

do 
   ...
end do 

このままだと無限ループになるので、通常はある条件を満たせばループから抜け出すようにする。

C言語だとbreak文でループから抜け出す。(正確には最も内側のループから抜け出す。) Fortran90/95でbreakに相当するのはexitである。

do 
  ...
  if ('''condition''') exit !---+
end do                      !   | ループから抜け出す。
                            !<--+

C言語ではループの先頭に戻る時にはcontinue文を使うが、 それに相当するのはFortran90/95ではcycleである。

do                           !<--+
  ...                        !   | 戻る
  if ('''condition''') cycle !---+
  ...
end do                      

【演習】

if文、select case文、do文を使って自由にプログラムを作ってみよ。

関数

C言語と同様に関数(function)という手続き(procedure)は、 Fortran90/95プログラムの重要な構成部品である。 関数の使い方は以下の例を見ればわかるであろう。

【演習】

  • 次のプログラムをfunction01.f95として保存し、実行せよ。
    integer function next_int(i)
      integer, intent(in) :: i
    
      next_int = i+1
    end function next_int
    
    program function01
      implicit none
      integer, external :: next_int
    
      print *, 'ans = ', next_int(-100)
    
    end program function01
    intent(in)は引数iが入力引数であり、関数next_intの内部で変更されることがないことを保証するものである。
  • 関数を以下のように内部副プログラムとしてメインの中に含むこともできる。 次のプログラムをfunction02.f95として保存し、実行せよ。 (内部副プログラムについては後述。)
    program function02
     implicit none
    
      print *, 'ans = ', next_int(-100)
    
    contains
    
      integer function next_int(i)
        integer, intent(in) :: i
    
        next_int = i+1
      end function next_int
    
    end program function02

サブルーチン

関数の本来の機能は入力に基づいて出力を返すものであるが、 特に出力(返り値)が不要な場合も多い。 そのようなときC言語ではvoid関数とする。 Fortran90/95ではvoid関数をsubroutineと呼ぶ。

サブルーチンを呼び出すにはcall subroutine名とする。 以下の例をみよ。

【演習】

subroutine next_int(input, output)
  integer, intent(in)  :: input
  integer, intent(out) :: output

  output = input + 1
end subroutine next_int

program subroutine01
  implicit none

  integer :: i, j
  i = 10
  call next_int(i,j)
  print *, 'ans = ', j

end program subroutine01

値渡しと参照渡し

C言語は値渡し、Fortran90/95は参照渡しである。 C言語が値渡しであることは、以下の例でkの値が変わっていないことから確認できる。

void increment(int i)
{
  i = 0;
  printf("i=%d\n", i);
}

main()
{
  int k = 2010;
  printf("before: k = %d\n", k);
  increment(k);
  printf("after:  k = %d\n", k);
}

これと似たプログラムをFortran90/95で書くと以下のようになる。

subroutine increment(i)    ! このプログラムは危険。
  integer :: i
  i = 0;
  print *, "i=", i
end subroutine increment

program call_by_reference
  implicit none
  integer :: k = 2010

  print *, "before: k = ", k
  call increment(k)
  print *, "after:  k = ", k
end program call_by_reference

このプログラムではkの値はサブルーチンincrementを読んだ後に変更されている。

【演習】

上の二つのプログラムをそれぞれcall_by_value.c、call_by_reference.f95という名前で保存し、実行せよ。

入出力属性

値渡しに比べて参照渡しは値をコピーする必要がないので速い。 従って計算速度を重視するFortran90/95言語では値渡しを採用していないのは当然である。 だが、上のcall_by_reference.f95プログラムで見たように、 引数として渡した変数がそのサブルーチン内部で勝手に(誤って)変更されると見つけにくいバグになり、危険である。 このような事故を防ぐためにFortran90/95には引数に入出力属性をつけることができる。 サブルーチンincrementの引数iは入力属性をもつと指定すると(下のプログラムの2行目)、 コンパイラがエラーを出してくれる。

subroutine increment(i)    ! このプログラムはコンパイルエラー
  integer, intent(in) :: i
  i = 0;
  print *, "i=", i
end subroutine increment

program call_by_reference
  implicit none
  integer :: k = 2010

  print *, "before: k = ", k
  call increment(k)
  print *, "after:  k = ", k
end program call_by_reference

【演習】

call_by_reference.f95を上のように変更してコンパイルせよ。

入出力属性

入出力属性には以下の3種類がある。

intent(in)入力その手続き内で値は変更されない変数
intent(out)出力その手続き内で値が設定される変数
intent(inout)入出力両者の混合。デフォルト。

バグの混入を防ぐために、Fortran90/95プログラムではすべての引数に入出力属性をつけることを強く勧める。

注釈


*1 Fortran90/95では、if (end==else) if=0と書いても文法的に問題ない。