計算物理学II (第６回)

今回の内容• 前回の復習(変数・代⼊⽂・四則演算・ループ・組込関数)• 配列• 変数を複数並べたもの。ベクトルや⾏列を表現• if⽂• プログラムは上から順番に⼀⽂ずつ実⾏するが、if⽂によって特定の条件を満たしたときだけ実⾏することが可能。
• do⽂(for⽂)と組み合わせて使うことで条件が満たされるまで繰り返す、などの制御も可能。

• 詳細は講義資料の⽂法のまとめを参照

プログラムの基本構造
Fortranprogram プログラム名
宣⾔部 (implicit none宣⾔、変数宣⾔)
実⾏部
end program プログラム名
プログラム名は英字から始まる適当な名前

C#include<ヘッダファイル名.h>int main(){
宣⾔部 (変数宣⾔など)
実⾏部 (計算など)}
プログラムはmain関数として実⾏される。
ヘッダファイルにはプログラム内でライブラリ関数などを⽤いる場合に読み込む。
数値計算で必要となるのはstdio.h,math.h,stdlib.h,complex.hなど変数• 変数：数値・⽂字などのデータを格納する領域。名前をつけて区別。
• 型があり、FortranやCでは宣⾔部で全ての型を明⽰的に指定する型

Fortran

C

サイズ

整数型

integerint

4バイト (= 32ビット)単精度実数型(物理の数値計算ではあまり使わない)

real, real*4, real(4)

float

4バイト

倍精度実数型

double precision, real*8,real(8)double8バイト

複素数型(物理の数値計算ではあまり使わない)

complex, complex*8

float complex(complex.hが必要)

8バイト

倍精度複素数型

complex(kind(0d0)),
complex*16, double complex

double complex(complex.hが必要)

16バイト

論理型

logicalbool (stdbool.hが必要)

trueかfalseか

⽂字型

characterchar

1⽂字1バイト

太字はよく使うもの。
変数名：英数字とアンダーラインが使える。英字から始まる四則演算と代⼊⽂• 数学の=(イコール)は左辺と右辺が等しいことを表すがプログラムの=は右辺の計算結果を左辺の変数に代⼊する、の意味• 四則演算は +, -, *, /• 累乗は Fortranでは** で表す。Cにはない (何度も掛けるか、数学関数のpowを使う)• 同じ型同⼠の演算は同じ型になる。
• 整数と実数の演算→実数、実数と複素数の演算→複素数
FortranC

a = 1.0d0
i =i + 5
a = 4**3
a= 5/2

a = 1.0;i = i + 5; または i += 5;
a = 4*4*4; または a = pow(4.0,3.0);
a = 5 / 2;
Cは⽂の終わりに;が必要

倍精度実数型のaに1.0を代⼊整数型のiの値を5増やす
aに4の3乗を代⼊aは整数型でも実数型でも2となる(整数と整数の演算は整数)
2.5にしたければ5.0/2や5.0/2.0、5/2.0などを代⼊実⾏⽂• プログラムの実⾏⽂は上から順番に⼀⾏ずつ実⾏される(重要)FortranCinteger :: a, b, c

int a, b, c;

整数変数a,b,cを宣⾔

a=2
b=5
c=a+b

a = 2;
b = 5;
c = a + b;

変数aに2を代⼊
変数bに5を代⼊
変数cにa + bを代⼊。この時点でaとbに値が⼊っているので2+5が計算されて変数cに7が代⼊される

integer :: a, b, c

int a, b, c;

c=a+b
a=2
b=5

順番を変えると結果が変わる。
cにa+bを代⼊。aとbに値がまだ⼊っていない。システム依存になるが
c = a + b; おそらくaにもbにも初期値の0が⼊っているのでcもおそらく0となるaに2を代⼊
a = 2;
bに5を代⼊b = 5;
cの値はその後変更していないので0のまま。
数学的には同じに⾒えるがプログラムでは結果は異なる。
プログラムを書くときは右辺に使う変数に値が⼊っているかを確認しながら書くループ(doループ/forループ)• 実⾏⽂は上から⼀⾏ずつ実⾏するのが原則だがよく使う例外は2つ(ループ、条件⽂)• doループ(Fortran)：do~end doで囲まれた領域をカウンタを変えて繰り返し実⾏• forループ(C) for( 式1(実⾏⽂); 式2(条件⽂); 式3(実⾏⽂)){繰り返し⽂}• 式1を実⾏(初期化)、式2の条件を判定、満たされていれば繰り返し⽂を実⾏• 式3を実⾏してから式2の条件判定に戻る

FortranCsum = 0
do i = 1, n
sum = sum + i
end dosum = 0;
for(i = 1;i<n ; i++){sum = sum + i;}を計算していることになっている
ループは和や積の計算、ベクトル、⾏列などindexを変えて⾏う計算、パラメータを変えて同じ計算を繰り返すときなどに使われる

標準⼊出⼒• 標準⼊⼒(キーボードからの⼊⼒)
FortranCread (*,*) var
read (5,*) var1, var2

scanf(“%d”, &var):scanf(“%lf %lf”, &var1, &var2);• 標準出⼒(画⾯への出⼒)• 標準エラー出⼒(画⾯への出⼒)
write (*,*) var
write (6,*) var
write (0,*) var
print *, var

printf(“%d”, var);printf(”var1 = %f,

C⾔語では⼊出⼒の型を明⽰する必要がある%d：整数型、%lf倍精度実数型
scanfのときには変数名の前に&をつける(変数のアドレス)var2 = %f¥n”, var1,var2);
C⾔語の単精度・倍精度実数の出⼒は%fを使う

fprintf(stderr, “var1 = %f¥n”, var1); //エラー出⼒

Fortran: read, writeの１つ⽬の引数でファイル装置番号を指定する。*とすると標準⼊出⼒になるが、Fortranでは5が標準⼊⼒、6が標準出⼒、0が標準エラー出⼒となっている
シェル側でリダイレクトを使うことで標準⼊出⼒先をファイルなどに変更することができる計算した結果はwriteかprintで必ず出⼒する

組込関数• 数学関数 (log, exp, sqrt, sin, cos, tanなど)は⽤意されている• ⾓度の単位はラジアン• C⾔語で数学関数を使うときの注意：• math.hをインクルードする• コンパイル時にオプション ‒lm をつける
• gcc ex4.c ‒lm

配列• 配列とは：変数が複数並んだもの(メモリ上でも並んで配置される)。
ベクトルや⾏列が表現できる。
配列の宣⾔の基本
Fortrandouble precision :: a(5)
a(1)

a(2)

a(3)

! aは5つの要素を持つ倍精度実数変数の配列

a(4)

a(5)

←メモリ上にこのように5つの倍精度実数変数の領域が連続して確保されるC
double a[5];
a[0]

a[1]

// aは5つの要素を持つ倍精度実数変数の配列
a[2]

a[3]

a[4]

←C⾔語では配列のラベルは0から始まる。
メモリ上にこのように5つの倍精度実数変数の領域が連続して確保される宣⾔部で設定

配列の初期値設定

FortranC

double precision :: a(1:5) = 0.0d0
integer :: k(1:5) = (/1, 2, 3, 4, 5/)double a[5] = {0.,0.,0.,0.,0.};int k[5] = {1, 2, 3, 4, 5};

実⾏⽂として設定してもよい(こちらのほうが⾃由度が⾼い)。ループを使う
do i = 1, 5a(i) = 0.0d0
end dofor (i=0;i<5;i++){a[i] = 0.0;}

異なる配列の違う添字を持つ要素の代⼊・演算も可能
do i = 1, 5a(i) = b(i+3)
end dofor(i=0;i<5;i++){a[i] = b[i+3];}

1次元配列の計算例
FortranC

real*8, dimension(1:5) :: a, b, c
real*8 :: adotb

double a[5], b[5], c[5];double adotb;! ここでaとbには値を代⼊しておく(省略)// ここでaとbには値を代⼊しておく(省略)ベクトル和do i = 1, 5c(i) = a(i) + b(i)end dodo i = 1, 5
write(*,*) “c(“, i, “)=“, c(i)end dofor( i = 0; i<5; i++){c[i] = a[i] + b[i];}for( i = 0;i<5;i++){printf(“c(%d)=%f¥n”, i, c[i]);}ベクトルの内積

adotb = 0.0d0
do i = 1, 5
adotb = adotb + a(i)*b(i)end doprint *, “a dot b = “, adotb

adotb = 0.0;
for( i = 0; i<5; i++){adotb += a[i]*b[i];}printf(“a dot b = %f¥n“, adotb);

1次元配列の組み込み関数(Fortran)
和c = sum(a(1:10))! a(1)からa(10)までの和をcに代⼊

最⼤値・最⼩値
d = maxval(a(1:n))e = minval(a(1:n))! a(1)からa(n)の中の最⼤値をdに代⼊! a(1)からa(n)の中の最⼩値をeに代⼊

2つの1次元配列の内積
f = dot_product(a(1:n), b(1:n)) ! n要素の配列aとbの内積をfに代⼊a, bが複素数型の場合は aの複素共役とbの積が計算される

⾏列• 2次元配列で⾏列が表現できる• 2次元配列の宣⾔と初期化
real*8, dimension(1:5,1:5) :: b
または real*8 :: b(5,5)
または real*8 :: b(1:5,1:5)など
do j = 1, 5
do i = 1, 5b(i,j) = 0.0d0end do
end dodouble b[5][5];for(i=0;i<5;i++){for(j=0;j<5;j++){b[i][j] = 0.0;}}

Fortran: b(1:5,1:5) = 0.0d0 !0を代⼊するだけなら1⾏でもかける。

2次元配列のメモリ上での配置
Fortranでの多次元配列 (b(5,5)の場合) はb(1,1)b(2,1)b(3,1)

b(4,1)b(5,1)b(1,2)b(2,2)b(3,2)b(4,2)

・・・

b[0][4]b[1][0]b[1][1]b[1][2]b[1][3]・・・

C⾔語での多次元配列(b[5][5]の場合)はb[0][0]b[0][1]b[0][2]b[0][3]の順でメモリ上に配置される。
doループ/forループで配列の値を使って計算する場合はメモリ上で連続的にアクセスするほうがよい(キャッシュミスが少なくなる)
例えばFortranではdoループは１つ⽬のインデックスを先に回すほうがよい
do j = 1, 5
do i = 1, 5b(i,j) = 0.0d0
end doend do⾏列の和・値の出⼒
FortranC

real*8, dimension(1:5,1:5) :: a, b, c

double a[5][5], b[5][5], c[5][5];! a, bに値をここで代⼊(省略)/* a, bに値をここで代⼊(省略) */do j = 1, 5
do i = 1, 5c(i,j) = a(i,j) + b(i,j)
end doend dofor(i = 0; i<5; i++){for(j = 0; j<5; j++){c[i][j] = a[i][j] + b[i][j];}}do i = 1, 5
write(*,*) (c(i,j), j = 1, 5)
end dofor (i = 0; i<5; i++){for(j = 0; j<5; j++){printf(“ %f “, c[i][j]);}printf(“¥n”);}

⾏列の積
⼀つの要素を計算するのにループが必要
C

Fortranreal*8, dimension(1:5,1:5) :: amat, bmat, cmat

double amat[5][5], bmat[5][5], cmat[5][5];! amatとbmatの値を代⼊/* amatとbmatの値をここで代⼊ */do j = 1, 5
do i = 1, 5
cmat(i,j) = 0.0d0
end doend dofor(i = 0; i<5; i++){for(j = 0; j<5; j++){cmat[i][j] = 0.0;

do j = 1, 5
do i = 1, 5
do k = 1, 5
cmat(i,j) = cmat(i,j) + amat(i,k) * bmat(k,j)
end doend doend do}}for (k = 0; k<5; k++){cmat[i][j] += amat[i][k] * bmat[k][j];}2次元配列の組み込み関数(Fortran)
転置⾏列
bmat(1:n,1:n) = transpose(amat(1:n,1:n))! bmatにamatの転置⾏列を代⼊

⾏列積
cmat(1:n,1:m) = matmul( amat(1:n,1:k), bmat(1:k,1:m))! (n×k⾏列のamat とk×m⾏列のbmatの積をcmatに代⼊(n×m⾏列)
matmulは経験上計算速度が早くないので⼤⾏列の計算では他のライブラリを使うほうがよい条件⽂• 実⾏⽂は上から⼀⾏ずつ実⾏するのが原則だがよく使う例外は2つだけ• if⽂：条件を論理式で指定し条件を満たすかどうかで処理内容を変える
FortranCinteger :: i

int i;

read (*,*) i

scanf(”%d”, &i);

if ( i<0) i=-i

if (i<0) i=-i;

print *, “i = “, iprintf(“i=%d¥n”,i);条件⽂

FortranC

if (論理式) 実⾏⽂

if (論理式) 実⾏⽂;

論理式が真のときのみ実⾏⽂が実⾏される(実⾏⽂が1⾏の場合)

if (論理式) then
実⾏⽂1
実⾏⽂2
end ifif (論理式){
実⾏⽂1;実⾏⽂2;}

論理式が真のときのみ実⾏⽂1と実⾏⽂2が実⾏される。
実⾏⽂が複数⾏ある場合はこのスタイル

if (論理式) then
実⾏⽂1
else
実⾏⽂2
end ifif (論理式){
実⾏⽂1;}else{
実⾏⽂2;}

論理式が真の場合は実⾏⽂1が実⾏され、偽の場合は実⾏⽂2が実⾏される

if (論理式1) then
実⾏⽂1
else if (論理式2) then
実⾏⽂2
else
実⾏⽂3
end ifif (論理式1){
実⾏⽂1;}else if (論理式2){
実⾏⽂2;}else{
実⾏⽂3;}

論理式1が真のときは実⾏⽂1を実⾏
論理式1が偽かつ論理式2が真の時は実⾏⽂2を実⾏
論理式1も論理式2も偽の場合は 実⾏⽂3を実⾏
else ifは複数あってもよい。

論理式
FortranFortran

Ca == ba .eq. ba == baがbと等しい時に真 (equal)

a /= ba .ne. b

a != baがbと等しくないとき真(not equal)

a >= ba .ge. ba >= b

aがb以上の時に真 (greater equal)a>b

a .gt. ba>baがbより⼤きい時に真(greater than)a <= ba .le. b

a <= baがb以下の時に真(less equal)a<ba .lt. ba<baがbより⼩さい時に真(less than)等しい時は=⼀つではなく2つなので注意

Fortran
意味
論理演算⼦

使⽤例

C
論理演算⼦

意味

使⽤例

.not.

(右)以外

.not. a==b!

(右)以外! a==b.and.

かつa==b .and. b==c&&

かつ

a==b && b==c.or.

もしくはa==b .or. b==c||

もしくは a==b || b==c.eqv.

論理値が等しい a<0 .eqv b<0 (aとbの符号が同じ).neqv.

論理値が異なる a<0 .neqv. b<0 (aとbの符号が違う)

プログラムの終了• プログラムは最終⾏まで到達すれば終了• 条件⽂を使うと途中で終了させることも可能• Fortran: stop⽂で終了。 stop “⽂字列” とすると⽂字列が標準エラー出⼒に出⼒• プログラムの⼀番最後にstopと正常終了のメッセージを書く場合もある• C: return⽂でmain関数から抜けてプログラム終了。異常終了の場合は0でない値を返す。
Fortranread(*,*) xif( x<0.0d0) stop “x must be zero or positive”
print *, “sqrt(x) = “, sqrt(x)! xを標準出⼒から読み込む! xが負であればプログラムを終了! xがゼロか正であれば以下の処理が⾏われるC
double x;
scanf(“%lf”, &x);// xを標準出⼒から読み込む
if( x<0.0){
fprintf(stderr, “x is negative¥n”);
return(-1);// xが負であれば-1を返してプログラムを終了}
printf(“sqrt(x) = %f¥n”, sqrt(x));// xがゼロか正であれば以下の処理が⾏われる条件⽂の例：クロネッカーのデルタ
FortranC

FortranC

integer :: delta, i, j

int delta, i, j;

integer :: delta, i, j

int delta, i, j;read (*,*) i, j

scanf(“%d %d”, &i, &j);read(*,*) i,j

scanf(“%d %d”, &i, &j);

if (i == j) thendelta = 1
else
delta = 0
end ifif (i == j){delta = 1;
}else{
delta = 0;}delta = 0if(i==j) delta = 1

delta = 0;if(i==j) delta = 1;

print *, delta

printf(“delta = %d¥n”, delta);

print *, delta

printf(“delta=%d¥n”, delta);doループとif⽂の組み合わせ
無限ループ

for(x=0.1; ; x+=0.1){// ここでいろいろ計算
a = ....;
if ( a<0 ) break;}do!ここでいろいろ計算
if ( a<0 ) exit
x=x+0.10d0
end doexit/break⽂はループから抜け、ループ⽂の直後に移動
条件が満たされない場合は計算が終わらないため実⾏してみて終わらない場合は強制終了(C-c)する

cycle(Fortran), continue(C) : ループの先頭まで戻り、カウンタを次の値に進める
do i = 1, nfor( i = 0;i<n;i++){a(i) = ....ここで何かを計算する

a[i] = ....; // ここで何かを計算する

if ( a(i)>0 ) cycle!a(i)が正なら戻ってカウンタを⼀つ回す

if (a[i]>0) continue; //a[i]>0なら戻ってカウンタを⼀つ回す! a(i)が0か負の場合だけこの領域で処理が⾏われる
end do/* a[i]が0か負の場合だけこの領域で処理が⾏われる */}演習• ⾃分で書けそうなら資料を⾒て、難しそうならサンプルのプログラムを写して書いてコンパイル、実⾏してください。
• 配列：第５回演習問題の(5)-(7)。条件⽂は第６回の(8)-(11)• 課題(提出任意)をmanabaに出しておきます。
• 第１回のレポート締切は11⽉22⽇(⽉)です。