こんにちは（@t_kun_kamakiri）
本記事ではFortranの変数とデータの型について解説を行います。
実例を交えながら確認した方が理解しやすいし、何より楽しいと思うのでこちらに書いた「マンデルブロー集合」のコードを例に確認していくことにします。
マンデルブロー集合のFortranコード
全体のコードがこちら。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | program main implicit none integer, parameter :: Nmax=50, Nx=1024, Ny=1024 double precision, parameter :: Lx=4.d0, Ly=3.d0, dx=Lx/dble(Nx), dy=Ly/dble(Ny) complex(kind(0d0)), parameter :: uso=(0.d0,1.d0) integer :: i, j, n, m double precision, allocatable, dimension(:) :: x, y complex(kind(0d0)) , allocatable, dimension(:,:) ::z, c character(len=40) filename allocate( x(0:Nx-1) , y(0:Nx-1), z(0:Nx-1,0:Ny-1), c(0:Nx-1,0:Ny-1)) call system("mkdir data") m = 0 !----------------------------------------------------------------------------------- z = 0 do j=0,Ny-1 do i=0,Nx-1 x(i) = -Lx/2.d0+dx*dble(i) y(j) = -Ly/2.d0+dy*dble(j) c(i,j) = x(i) + uso*y(j) do n=0,Nmax if(cdabs(z(i,j))>2.0d0 .or. n == Nmax) exit z(i,j) = z(i,j)**2 + c(i,j) end do ! 発散しなかったCの集合(Cを上書き) if (n == Nmax) then c(i,j) = 0.0 else ! 発散したCの集合(Cを上書き) c(i,j) = dble(n) end if end do end do deallocate(z) !----------------------------------------------------------------------------------- !output of movie file write(filename,"(a,i5.5,a)") "data/julia_",int(m),".vtk" open(10,file=filename) write(10,"('# vtk DataFile Version 3.0')") write(10,"('test')") write(10,"('ASCII ')") write(10,"('DATASET STRUCTURED_GRID')") write(10,"('DIMENSIONS ',3(1x,i4))") Nx, Ny, 1 write(10,"('POINTS ',i9,' float')") Nx*Ny*1 do j=0,Ny-1 do i=0,Nx-1 write(10,"(3(f9.4,1x))") x(i), y(j), 0.d0 end do end do write(10,"('POINT_DATA ',i9)") Nx*Ny*1 !date input write(10,"('SCALARS julia float')") write(10,"('LOOKUP_TABLE default')") do j=0,Ny-1 do i=0,Nx-1 write(10,*) cdabs(c(i,j)) end do end do close(10) end program main |
こちらを実行することでマンデルブロー集合を再現できます。
美しいジュリア集合
と
マンデルブロー集合 pic.twitter.com/yoxZv8RzPe— カマキリ🐲Python頑張る昆虫 (@t_kun_kamakiri) April 11, 2022
コードはそこまで長くないのでべた書きしていますが、追々改良を進めていきます。
繰り返し文（do文）
Fortranの繰り返し文（do文）の基本構文は以下のように記述します。
1 2 3 | do 変数=初期値,最終値[,刻み幅] 繰り返したい処理 end do |
本コードでの繰り返し文は以下の部分です。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | z = 0 do j=0,Ny-1 do i=0,Nx-1 x(i) = -Lx/2.d0+dx*dble(i) y(j) = -Ly/2.d0+dy*dble(j) c(i,j) = x(i) + uso*y(j) do n=0,Nmax if(cdabs(z(i,j))>2.0d0 .or. n == Nmax) exit z(i,j) = z(i,j)**2 + c(i,j) end do ! 発散しなかったCの集合(Cを上書き) if (n == Nmax) then c(i,j) = 0.0 else ! 発散したCの集合(Cを上書き) c(i,j) = dble(n) end if end do end do |
ここではx,y,z,cに対しては以下の配列を定義しているので、アクセスする際にはdo文を使って繰り返し構文を書いています。
1 2 3 4 5 | double precision, allocatable, dimension(:) :: x, y complex(kind(0d0)) , allocatable, dimension(:,:) ::z, c !配列の割付け allocate( x(0:Nx-1) , y(0:Nx-1), z(0:Nx-1,0:Ny-1), c(0:Nx-1,0:Ny-1)) |
ここで一点注意ですが、Fortranの配列の並び順がarray(1,1),array(2,1),array(3,1),array(4,1),array(1,2)・・・・であるため以下のように1つ目の要素から順番にアクセスをしていった方が効率が良いプログラムになります。
1 2 3 4 5 6 | do j = 1, n do i =1, n a(i,j) = 10* i + j write(*,*)i,j,"",a(i,j) end do end do |
Fortranでの離れたメモリアクセスは実行速度が落ちる
では、実際に「離れたメモリでアクセス」と「メモリの並び順でアクセス」でCPUでの実行速度を計測してみましょう。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 | program main implicit none integer, allocatable ::a(:,:) integer, parameter::n=2 integer::i, j real(kind=8)::t1, t2 allocate(a(n,n)) write(*,*)"========================" write(*,*)"離れたメモリでアクセス" ! 行列の作成 call cpu_time(t1) do i = 1, n do j =1, n a(i,j) = 10* i + j write(*,*)i,j,"",a(i,j) end do end do call cpu_time(t2) write(*,*)"CPU time = ", t2 - t1 ! 全要素を出力 write(*,*)"全要素の出力" write(*,*)a write(*,*)"========================" write(*,*)"メモリの並び順でアクセス" ! 行列の作成 call cpu_time(t1) do j = 1, n do i =1, n a(i,j) = 10* i + j write(*,*)i,j,"",a(i,j) end do end do call cpu_time(t2) write(*,*)"CPU time = ", t2 - t1 ! 全要素を出力 write(*,*)"全要素の出力" write(*,*)a deallocate(a) ! メモリの解放 end program main |
結果はこちらです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | $ gfortran F_list29.f90 $ ./a.out ======================== 離れたメモリでアクセス 1 1 11 1 2 12 2 1 21 2 2 22 CPU time = 1.0399999999999993E-004 全要素の出力 11 21 12 22 ======================== メモリの並び順でアクセス 1 1 11 2 1 21 1 2 12 2 2 22 CPU time = 3.7999999999999839E-005 全要素の出力 11 21 12 22 |
まずwrite(*,*)aで出力した結果が 11 21 12 22となっていることに注意します。
これはarray(1,1),array(2,1),array(1,2),array(2,2)の並び順で配列ができていることを意味しています。
そして結果について、
- 離れたメモリでアクセスした場合：CPU time = 1.0399999999999993E-004
- メモリの並び順でアクセスした場合：CPU time = 3.7999999999999839E-005
となり、計算時間は1/3になっています。
今回のような小さな配列では大した差は出てこないですが、3次元配列、4次元配列となるにつれて離れたメモリへのアクセスは実行速度を落とすことになります。
制御構文（if文）
Fortranでのif文の基本構文は以下のようになります。
1 2 3 4 5 6 7 8 9 | if (条件 1) then 条件 1 が真の場合の処理 else if (条件 2) then 条件 1 が偽で条件 2 が真の場合の処理 else if (条件 n) then 条件 1 から条件 n-1 が偽で条件 n が真の場合の処理 else 条件 1 から条件 n がすべて偽の場合の処理 end if |
本コードでの制御構文は以下の部分です。
1 2 3 4 5 6 7 | ! 発散しなかったCの集合(Cを上書き) if (n == Nmax) then c(i,j) = 0.0 else ! 発散したCの集合(Cを上書き) c(i,j) = dble(n) end if |
参考書
Fortranは日本の書籍がかなり少ないのですが以下のようなサイトがあり、とても参考になります。
その他出版されている参考書を挙げておきます。
こちらの参考書は文法からちょっとした応用（偏微分方程式を解く）まで解説がありFortranを使う人にとって手元に置いておきたい参考書ですね。OpenMPによる並列計算の解説もあります。
しかし、Fortran90/95の記述のみなのでForran2003以降を学ぶにはやはり上記のサイトを参考にした方が良いでしょう。
こちらはFortranのコード集となっています。
Fortranは科学計算のためのプログラミング言語ですので数値計算の参考コードがいっぱい載っているこちらの参考書はとても勉強になります。
Fortran文法の解説はありませんが、数値計算の勉強にはとても良いです。
こちらはC言語とFortranの両方を学びながら数計算の基礎力を身に付けるための参考書です。問題演習の全てにサンプルコードがあるため本書を読みながらサンプルコードをじっくり眺めるだけでも大変力が付きます。
Follow @t_kun_kamakiri 