(H29. 6. 8 CoArray 部分追加)
(H30.11.24 FORTRAN II 部分追加)
#前置き
ここではマンデルブロ (Mandelbrot) 集合を描くプログラムを、昔風の FORTRAN から今時の Fortran に徐々に書き換えて行くことで、今時の Fortran の新機能や発想の一部を紹介してゆきたいと思います。
###Fortran の現代的意義
Intel や Nvidia のようなプロセッサメーカーは、自社のプロセッサの性能を最大に引き出すコンパイラを独自に開発していて、出荷に合わて利用者に提供しています。そのコンパイラの言語はいずれも C/C++ および Fortran です。（追記 H30：AMD や ARM も、それぞれ C/C++ および Fortran コンパイラを提供しています。）
NEC や Fujitsu, IBM, Cray 等のスーパーコンピュータでも、ハードウェア性能を最大限に引き出すために提供される言語は C/C++ および Fortran となっています。また分散並列処理のための MPI ライブラリも C および Fortran 用のみが開発されています。（参考：CPU/GPU プロセッサ メーカー謹製 Fortran）
つまり、これらの最新ハードウェアで実行時性能を最大限に引き出すためには、C または Fortran を用いることになります。特に Fortran は、C よりも少ない学習コストで容易に高い実行時性能を引き出せるように設計されています。
自然科学や工学分野を専門とする場合、プログラムは電子計算機を思い通りに使役する法術に見なせます。プログラム言語は、専門分野固有の抽象を数値的に実現する際に一時的に存在する中間形式なので、元々の抽象構造をそのまま記述できる方が好ましいことになります。また欲しいのは計算結果なので、ソースプログラムの形式よりも実行プログラムが高速に実行されることが重要視されます。
以上のことに鑑みると、最新ハードウェア上で最大限に性能を活用した計算を、少ないコストで実現したい場合、 Fortran の利用が合理的な選択になります。ここにその存在意義があります。
#今時の Fortran （Modern Fortran)
今時の Fortran すなわち Modern Fortran とは、Fortran90 およびそれ以降の Fortran を指します。具体的には Fortran90/95/2003/2008 を指します。以下にその特徴を記します。
###Fortran90 構造化・並列性
Module、全配列操作(Whole Array Operation)、動的メモリー確保/解放(allocate/deallocate)
###Fortran95 並列命令の強化（共有メモリー）
forall、純粋(pure)/要素毎(elemental)副プログラム、maskの拡張、deallocate 自動化 (allocatable のみ、pointer は手動）
###Fortran2003 動的プログラミング
オブジェクト指向、動的メモリー確保の強化（代入時自動確保）
###Fortran2008 並列命令の強化（分散メモリー）
coarray、do concurrent、block 構造(thread safe)
###並列性
今時の Fortran には並列性という概念があります。並列性は普通、並行動作による処理速度の向上を目的としますが、ここではアルゴリズムの本質的な順序構造の明示的記述を目的とします。（並列化はオーバーヘッドが大きいので、規模の小さい計算では実行速度は向上しないのが普通です。）
コンピュータ言語では、漸化式のように本質的に順序に依存した反復操作も、配列への代入のように本来順序によらない複数の要素に対する操作も、同じ反復構文を用いて記述することが多いのですが、並列性の観点から見ればこれらは明示的に区別されることが望ましいと思われます。今時の Fortran では、順序構造をもった反復と順序に依存しない反復を明示的に区別することが可能です。
#マンデルブロ集合 (Mandelbrot set)
マンデルブロ集合は複素数の漸化式$z_{n+1}=z_n^2+c$、$z_0=0$から得られる集合で、複素平面上の各点$c$毎の発散/収束性をプロットすることで、特徴的な自己相似図形が得られます。$|z_n|>2$となると必ず発散することが分かっているので、発散条件を満たすまでの反復数 $n$ によって色分けするとより印象的な図となります。
なお複素平面上の座標点毎に計算は独立なので、漸化式の反復と座標についての反復操作は独立で順序が入れ替え可能（可換）になっています。つまり並列に実行できる構造になっています。
#ソースプログラム
##昔風の FORTRAN
まずは出発点となる昔風の FORTRAN プログラムを見てみましょう。
ここでは簡単のためまずキャラクターグラフィクスによってマンデルブロ集合を表示することにします。
###FORTRAN77
FORTRAN77 は現在も広く利用されている言語で、LAPACK などの国家予算で維持されるような膨大なライブラリ資産があります。静的なプログラミング言語であるためコンパイラによる自動最適化も容易で、高速な実行プログラムが得られます。
PROGRAM MANDEL
PARAMETER(MAXIT = 90, NX = 61, NY = 31)
PARAMETER(X0 = -2.0, X1 = 2.0, Y0 = -2.0, Y1 = 2.0)
DIMENSION MNDLBR(NX, NY)
COMPLEX C, Z
DO 10 IY = 1, NY
Y = Y0 + (Y1 - Y0) * (IY - 1) / REAL(NY - 1)
DO 20 IX = 1, NX
X = X0 + (X1 - X0) * (IX - 1) / REAL(NX - 1)
Z = (0.0, 0.0)
C = CMPLX(X, Y)
DO 30 ITER = 0, MAXIT
Z = Z * Z + C
IF (ABS(Z) .GT. 2.0) GO TO 99
30 CONTINUE
99 MNDLBR(IX, IY) = ITER
20 CONTINUE
10 CONTINUE
C
DO 40 IY = 1, NY
PRINT '(1H ,61I1)', ((MNDLBR(IX, IY) + 9) / 10, IX = 1, NX)
40 CONTINUE
STOP
END
プログラムの基本的構造は三重ループの反復からなっています。座標に関する反復を外側に置いて、各点毎に順番に漸化式の反復を行って発散性を求めています。ソースコード上では x y 座標の方向に関する等価性が崩れています。また並置的な座標と順序に依存する漸化式が同じ DO LOOP 構文による反復であらわされており、計算が本来もつ対称性の情報が失われています。
FORTRAN77 以降では DO LOOP を抜けた後の DO 変数の値は抜けた時の値が保たれることになっているので発散までの反復数を数えるのにそれを利用しています。FORTRAN66 以前では不定でした。発散しないことは無限回反復しなければ分からないのですが、ここでは90回の反復で止めています。反復数は図形の精度に依存して変える必要があります。
FORTRAN77 までは、出力の１文字目が出力制御文字として利用されており、単なる出力の場合１文字目を空白にする必要があります。ここでは出力 FORMAT を文字列の形で直接 PRINT 文に与えているので、文字列中に文字列を書く目的でホレリス型の表現'1H 'を用いて１文字の空白を表しています。
FORMAT は整数１桁にしてあるため、２桁以上の数がくると桁を溢れて画面に * が出力されることになります。これを利用して 0 は１回の反復で発散確定。1 は２～１１回、2 は１２～２１回・・・、9 は８２～９０回、＊は９０回の反復では発散が確定しなかったことを表すようになっています。
実行結果
MS-FORTRAN POWER STATION 1.0 による実行例
画面の数字は発散確定($|z|>2$)までの反復数を反映しています。半径 2 の円の外側は計算するまでもなく発散が確定しています。
###FORTRAN66
FORTRAN66 で新たにプログラムが書かれることはほぼないのですが、伝説の超古代文明の遺跡のように、あちこちで断片的に生き残って利用されていたりします。
C MANDELBROT SET
DIMENSION MNDLBR(61, 31)
COMPLEX C, Z
NX = 61
NY = 31
MAXIT = 91
X0 = -2.0
X1 = 2.0
Y0 = -2.0
Y1 = 2.0
DO 10 IY = 1, NY
Y = Y0 + (Y1 - Y0) * FLOAT(IY - 1) / FLOAT(NY - 1)
DO 20 IX = 1, NX
X = X0 + (X1 - X0) * FLOAT(IX - 1) / FLOAT(NX - 1)
Z = (0.0, 0.0)
C = CMPLX(X, Y)
NITER = 0
DO 30 ITER = 1, MAXIT
Z = Z * Z + C
IF (REAL(Z)**2 + AIMAG(Z)**2 .GT. 4.0) GO TO 99
NITER = NITER + 1
30 CONTINUE
99 MNDLBR(IX, IY) = (NITER + 9) / 10
20 CONTINUE
10 CONTINUE
C
DO 40 IY = 1, NY
WRITE(6, 100) (MNDLBR(IX, IY), IX = 1, NX)
40 CONTINUE
100 FORMAT(1H , 61I1)
STOP
END
実行結果
IBM 360 エミュレータ上の FORTRAN IV H による実行例
（参考：http://fortran66.hatenablog.com/entry/20090217/1234892880)
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000001111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72681111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*728*1111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000001111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
基本的なプログラム構造は 77 の時と同じですが、DO LOOP 変数が正でなければならいことや、DO LOOP を抜けた後の DO LOOP 変数の値が不定であること、出力文等での演算が制限されること、PARAMETER 文がなく値をハードコーディングする必要があること等、細かな相違点に注意が必要です。
###（付録) FORTRAN II
FORTRAN I/II(/III) は、個別機器毎の言語仕様になっていました。FORTRAN IV からハードウェアの抽象化による、機器と独立した言語仕様を目指しました。
IBM709/7090 のような高級機では、第一カラムに特定文字を書くことで、倍精度や複素演算ができました。プログラムもわずかな修正で実行できます。
実行結果
IBM 7090 エミュレータ上の FORTRAN II による実行例
https://qiita.com/cure_honey/items/9122cc3a21acf4dfe8f3
C MANDELBROT SET
DIMENSION MNDLBR(61, 31)
DIMENSION C(2), Z(2)
NX = 61
NY = 31
MAXIT = 91
X0 = -2.0
X1 = 2.0
Y0 = -2.0
Y1 = 2.0
DO 10 IY = 1, NY
Y = Y0 + (Y1 - Y0) * FLOATF(IY - 1) / FLOATF(NY - 1)
DO 20 IX = 1, NX
X = X0 + (X1 - X0) * FLOATF(IX - 1) / FLOATF(NX - 1)
I Z = (0.0, 0.0)
C(1) = X
C(2) = Y
NITER = 0
DO 30 ITER = 1, MAXIT
I Z = Z * Z + C
IF (Z(1)**2 + Z(2)**2 - 4.0) 30, 40, 40
30 NITER = NITER + 1
40 CONTINUE
99 MNDLBR(IX, IY) = (NITER + 9) / 10
20 CONTINUE
10 CONTINUE
C
C
DO 50 IY = 1, NY
PRINT 100, (MNDLBR(IX, IY), IX = 1, NX)
50 CONTINUE
100 FORMAT(1H ,61I1)
STOP
END
0000000000000000000000000000000000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000001111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
0001111111111111111111111110002111111111111111111111111111000
0011111111111111111113312223072671111111111111111111111111100
0011111111111111111111000000000000061111111111111111111111100
0111111111111111111120000000000000001111111111111111111111110
0111111111112040341200000000000000002111111111111111111111110
0111111112120000000000000000000000005111111111111111111111110
0345555680000000000000000000000000311111111111111111111111111
0111111112120000000000000000000000005111111111111111111111110
0111111111112040341200000000000000002111111111111111111111110
0111111111111111111120000000000000001111111111111111111111110
0011111111111111111111000000000000061111111111111111111111100
0011111111111111111113312223072671111111111111111111111111100
0001111111111111111111111110002111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000001111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
##今時の Fortran
次に今時の Fortran プログラムを見てゆくことにします。
昔風の FORTRAN プログラムをそのまま移した Fortran90 プログラムから始めて、徐々に今風の機能を取り入れてゆき、最後に BMP ファイルにマンデルブロ集合を描き出したいと思います。
###Fortran90
ここでは文法的な書き方を改めた以外は本質的な変更を加えていません。1990年代においては、小文字は変数名などユーザー定義のものに限り、Fortran 固有の命令は大文字で書くことが奨励されていました。
PROGRAM mandel
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 61, ny = 31, maxiter = 90
REAL , PARAMETER :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
INTEGER :: ix, iy, iter, mandelbrot(nx, ny)
REAL :: x, y
COMPLEX :: z, c
DO iy = 1, ny
y = y0 + (y1 - y0) * (iy - 1) / REAL(ny - 1)
DO ix = 1, nx
x = x0 + (x1 - x0) * (ix - 1) / REAL(nx - 1)
c = CMPLX(x, y)
z = (0.0, 0.0)
DO iter = 0, maxiter
z = z * z + c
IF (ABS(z) > 2.0) EXIT
END DO
mandelbrot(ix, iy) = iter
END DO
END DO
!
DO iy = 1, ny
PRINT '(61I1)', (mandelbrot(:, iy) + 9) / 10
END DO
STOP
END PROGRAM mandel
実行結果
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000111111111111111111111111111111111110000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000000111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72671111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*72661111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000000111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
続行するには何かキーを押してください . . .
###Fortran90 Moduleの使用
次に module を使用した例を示します。module とは subroutine や function といった副プログラムを束ねる入れ物のようなものです。この中に副プログラムを収めておくと、インターフェースが自動で生成されるため、コンパイル時に引数のチェックなどが行われます。FORTRAN77 での動作不良の多くに、副プログラムの引数の不整合が多かった事を思うと、非常に便利な仕組みです。
MODULE m_mandel
IMPLICIT NONE
INTEGER, PARAMETER :: maxiter = 90
CONTAINS
INTEGER FUNCTION mandel(c)
COMPLEX, INTENT(IN) :: c
COMPLEX :: z
z = (0.0, 0.0)
DO mandel = 0, maxiter
z = z * z + c
IF (ABS(z) > 2.0) EXIT
END DO
END FUNCTION mandel
END MODULE m_mandel
PROGRAM mandel_main
USE m_mandel
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 61, ny = 31
REAL , PARAMETER :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
INTEGER :: ix, iy, iter, mandelbrot(nx, ny)
REAL :: x, y
DO iy = 1, ny
y = y0 + (y1 - y0) * (iy - 1) / REAL(ny - 1)
DO ix = 1, nx
x = x0 + (x1 - x0) * (ix - 1) / REAL(nx - 1)
mandelbrot(ix, iy) = mandel(CMPLX(x, y))
END DO
END DO
!
DO iy = 1, ny
PRINT '(61i1)', (mandelbrot(:, iy) + 9) / 10
END DO
STOP
END PROGRAM mandel_main
Fortran には subroutine と function という二つの副プログラム形式が存在しています。subroutine は値を返さず、call 文で呼び出します。function は値を返し、代入の右辺などに現れます。
これらの使い分けは、function では暗黙に純粋性が要求されている点にあります。つまり function では引数の値を変更したり、大域変数等を変更したり依存してはならないという点です。
このことは組み込みサブルーチンや関数でも徹底されており、たとえば Fortran90 では、時刻や乱数を求める副プログラムは多くの言語と異なって subroutine になっています。
Fortran では引数が参照型であること、最適化のために演算の順序を入れ替えてよいことがあるため、function に依存性があると最適化が大きく抑止されてしまいます。そのため function の純粋性の要請が最初期からなされていました。（例えば y=sin(x)+sin(x)+sin(x) を y=3*sin(x) にまとめられなくなる。）
ただし文法上では純粋でない function も禁止されておらず、危険性を理解したうえで利用することは可能です。
実行結果
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000111111111111111111111111111111111110000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000000111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72671111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*72661111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000000111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
続行するには何かキーを押してください . . .
###Fortran95
fortran95 では並列計算のための機能が強化されて forall のような新しい構文が導入されました。またこれに関連して、純粋でない依存性のある subroutine や function は、並列動作させることができないので、プログラム上で純粋性の保証をする修飾子 pure が新たに導入されました。
module m_mandel
implicit none
integer, parameter :: maxiter = 90
contains
pure elemental integer function mandel(c)
complex, intent(in) :: c
complex :: z
z = (0.0, 0.0)
do mandel = 0, maxiter
z = z * z + c
if (abs(z) > 2.0) exit
end do
end function mandel
end module m_mandel
program mandel_main
use m_mandel
implicit none
integer, parameter :: nx = 61, ny = 31
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
integer :: ix, iy, iter, mandelbrot(nx, ny)
real :: x(nx), y(ny)
complex :: c(nx, ny)
forall (ix = 1:nx) x(ix) = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
forall (iy = 1:ny) y(iy) = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(x(ix), y(iy))
mandelbrot = mandel(c)
!
do iy = 1, ny
write(*, '(61i1)') (mandelbrot(:, iy) + 9) / 10
end do
stop
end program mandel_main
さて Fortran90 の組み込み関数やサブルーチンの多くは、要素毎演算 (elemental) に対応していました。これは、スカラーに対して定義してある引数に配列を与えることができるというもので、MAP 演算に相当するものです。
Fortran95 では、新たに修飾子 elemental が導入され、ユーザー定義の要素毎演算副プログラムが書けるようになりました。elemental な副プログラムは pure であることが要求されます。これは若干冗長に思われますが、Fortran2015 では、外部ファイルへの I/O を伴う impure な elemental 副プログラムの導入が予定されているので、将来的には意味があります。
forall は HPC ベンダー界の悪夢ともいわれる構文で、Fortran の配列代入の自由度の高さから、実装するのが難しい割に並列性能が大してでません。したがって使用を控えるべきとする意見もあります。ここではごく単純な場合に限って、並列性の明示のために利用することにします。
さて forall 構文は代入文の拡張であって do loop のような一般的反復の並列化とは異なります。forall (ix = 1:nx, iy = 1:ny) の場合、インデックス ix, iy は全く順序によらず実行されてよいことを意味します。つまりこれらの座標の独立性・等価性をあらわしています。
一方で要素毎演算を利用しているので、計算の本体は mandelbrot = mandel(c) という１つの式で表わされるようになり、座標に対する反復をあらわに記述する必要はなくなりました。これにより元々の定式化に近い形でプログラムが記述できています。
実行結果
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000111111111111111111111111111111111110000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000000111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72671111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*72661111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000000111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
続行するには何かキーを押してください . . .
###Fotrtran2008
次に Fortran2003 を飛ばして Fortran2008 で導入された並列化構文を使ってみます。Fortran95 で導入された forall 構文は代入文であり使いにくいものであったため、do loop の並列実行構文として do concurrent 構文が導入されました。
module m_mandel
implicit none
integer, parameter :: maxiter = 90
contains
pure integer function mandel(c)
complex, intent(in) :: c
complex :: z
z = (0.0, 0.0)
do mandel = 0, maxiter
z = z * z + c
if (abs(z) > 2.0) exit
end do
end function mandel
end module m_mandel
program mandel_main
use m_mandel
implicit none
integer, parameter :: nx = 61, ny = 31
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
integer :: ix, iy, iter, mandelbrot(nx, ny)
do concurrent (ix = 1:nx, iy = 1:ny)
block
real :: x, y
x = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
y = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
mandelbrot(ix, iy) = mandel(cmplx(x, y))
end block
end do
!
do iy = 1, ny
write(*, '(*(i1))') (mandelbrot(:, iy) + 9) / 10
end do
stop
end program mandel_main
do concurrent 構文の中では変数がスレッドセーフでなければならないため、それを保証する block 構文も導入されました。ここでは座標 x, y が並列動作する各インデックス毎に独立でなければならないので block 構文内で局所的に宣言する必要があります。
Format には Fortran2008 で導入された無限回反復指定子 *( ) を用いています。
実行結果
0000000000000000000000000000001000000000000000000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000111111111111111111111111111111111110000000000000
0000000000111111111111111111111111111111111111111110000000000
0000000011111111111111111111111111111111111111111111100000000
0000000111111111111111111111111111111111111111111111111000000
0000011111111111111111111111111111111111111111111111111100000
0000111111111111111111111111211111111111111111111111111110000
000111111111111111111111111***2111111111111111111111111111000
0011111111111111111113312223*72671111111111111111111111111100
0011111111111111111111*************61111111111111111111111100
011111111111111111112***************1111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111212************************5111111111111111111111110
**********************************311111111111111111111111111
011111111212************************5111111111111111111111110
0111111111112*4*3412****************2111111111111111111111110
011111111111111111112***************1111111111111111111111110
0011111111111111111111*************61111111111111111111111100
0011111111111111111113312223*72661111111111111111111111111100
000111111111111111111111111***2111111111111111111111111111000
0000111111111111111111111111211111111111111111111111111110000
0000011111111111111111111111111111111111111111111111111100000
0000000111111111111111111111111111111111111111111111111000000
0000000011111111111111111111111111111111111111111111100000000
0000000000111111111111111111111111111111111111111110000000000
0000000000001111111111111111111111111111111111111000000000000
0000000000000000111111111111111111111111111110000000000000000
0000000000000000000011111111111111111111100000000000000000000
0000000000000000000000000000001000000000000000000000000000000
続行するには何かキーを押してください . . .
###Fortran2003
はじめに書いたように、本来、座標に関する反復と漸化式に関する反復は独立で入れ替え可能です。ここでは少し趣向を変えて、漸化式に関する反復を最外ループとしてプログラムを書いていみます。
program mandel
implicit none
integer, parameter :: maxiter = 90, nx = 61, ny = 31
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
character(len = 10), parameter :: text = '|+o0O.@*#-'
integer :: ix, iy, iter, mandelbrot(nx, ny) = 0, m(nx)
real, allocatable :: x(:), y(:)
complex :: c(nx, ny), z(nx, ny) = (0.0, 0.0)
x = [( (x1 - x0) / (nx - 1) * (ix - 1) + x0, ix = 1, nx )]
y = [( (y1 - y0) / (ny - 1) * (iy - 1) + y0, iy = 1, ny )]
forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(x(ix), y(iy))
do iter = 0, maxiter
where (abs(z) <= 2.0)
z = z * z + c
mandelbrot = mandelbrot + 1
end where
end do
do iy = 1, ny
m = (mandelbrot(:, iy) + 8) / 10 + 1
write(*, '(61a1)') (text(m(ix):m(ix)), ix = 1, nx)
end do
stop
end program mandel
where 構文は Fortran90 で導入された構文で配列の代入操作に関するマスク操作を行います。漸化式は全配列操作によって、座標に関する反復を明示することなく行っています。
座標点は内包 do loop と配列構成子を用い生成し、Fortran2003 で導入された割り付け配列への代入による自動割り付けを用いて、適切な大きさの配列にを割り付けながら代入しています。
実行結果
ここでは反復を数字の代わりに文字であらわして作図してみました。
||||||||||||||||||||||||||||||+||||||||||||||||||||||||||||||
||||||||||||||||||||+++++++++++++++++++++||||||||||||||||||||
||||||||||||||||+++++++++++++++++++++++++++++||||||||||||||||
||||||||||||++++++++++++++++++++++++++++++++++++|||||||||||||
||||||||||+++++++++++++++++++++++++++++++++++++++++||||||||||
||||||||+++++++++++++++++++++++++++++++++++++++++++++||||||||
||||||++++++++++++++++++++++++++++++++++++++++++++++++|||||||
|||||+++++++++++++++++++++++++++++++++++++++++++++++++++|||||
||||++++++++++++++++++++++++o++++++++++++++++++++++++++++||||
|||++++++++++++++++++++++++---o+++++++++++++++++++++++++++|||
||+++++++++++++++++++00+ooo0-*o@*++++++++++++++++++++++++++||
||++++++++++++++++++++-------------@+++++++++++++++++++++++||
|+++++++++++++++++++o---------------++++++++++++++++++++++++|
|+++++++++++o-O-0O+o----------------o+++++++++++++++++++++++|
|++++++++o+o------------------------.+++++++++++++++++++++++|
----------------------------------0++++++++++++++++++++++++++
|++++++++o+o------------------------.+++++++++++++++++++++++|
|+++++++++++o-O-0O+o----------------o+++++++++++++++++++++++|
|+++++++++++++++++++o---------------++++++++++++++++++++++++|
||++++++++++++++++++++-------------@+++++++++++++++++++++++||
||+++++++++++++++++++00+ooo0-*o@@++++++++++++++++++++++++++||
|||++++++++++++++++++++++++---o+++++++++++++++++++++++++++|||
||||++++++++++++++++++++++++o++++++++++++++++++++++++++++||||
|||||+++++++++++++++++++++++++++++++++++++++++++++++++++|||||
|||||||+++++++++++++++++++++++++++++++++++++++++++++++|||||||
||||||||+++++++++++++++++++++++++++++++++++++++++++++||||||||
||||||||||+++++++++++++++++++++++++++++++++++++++++||||||||||
|||||||||||||+++++++++++++++++++++++++++++++++++|||||||||||||
||||||||||||||||+++++++++++++++++++++++++++++||||||||||||||||
||||||||||||||||||||+++++++++++++++++++++||||||||||||||||||||
||||||||||||||||||||||||||||||+||||||||||||||||||||||||||||||
続行するには何かキーを押してください . . .
###Fortran2003 OOP
ここではオブジェクト指向を取り入れて Bitmap ファイルとその書き出しルーチンを定義します。またユーザー定義派生型 I/O も定義します。
module m_bmp
implicit none
type :: t_bmp_file_header
sequence ! 14bytes
character(len = 2) :: bfType = 'BM' !integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
integer(4) :: bfSize ! file size in bytes
integer(2) :: bfReserved1 = 0 ! always 0
integer(2) :: bfReserved2 = 0 ! always 0
integer(4) :: bfOffBits
end type t_bmp_file_header
!
type :: t_bmp_info_header
sequence
integer(4) :: biSize = Z'28' ! size of bmp_info_header ; 40bytes
integer(4) :: biWidth
integer(4) :: biHeight
integer(2) :: biPlanes = 1 ! always 1
integer(2) :: biBitCount
integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
integer(4) :: biSizeImage
integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
integer(4) :: biYPelsPerMeter = 3780 ! 96dpi
integer(4) :: biClrUsed = 0
integer(4) :: biClrImportant = 0
end type t_bmp_info_header
!
type :: t_rgb
sequence
character :: b, g, r ! order is b g r
end type t_rgb
!
type :: t_bmp(nx, ny)
integer, len:: nx, ny
type(t_rgb) :: rgb(nx, ny)
contains
procedure :: wr => wr_bmp
procedure :: pr_bmp
generic :: write(formatted) => pr_bmp
end type
contains
subroutine wr_bmp(bmp, fn)
class(t_bmp(*, *)), intent(in) :: bmp
character(len = *), intent(in) :: fn
type(t_bmp_file_header) :: bmp_file_header
type(t_bmp_info_header) :: bmp_info_header
bmp_file_header%bfSize = 14 + bmp_info_header%biSize + 0 + bmp%nx * bmp%ny * 3
bmp_file_header%bfOffBits = 14 + bmp_info_header%biSize
bmp_info_header%biWidth = bmp%nx ! nx shouold be a multiple of 4
bmp_info_header%biHeight = bmp%ny
bmp_info_header%biBitCount = 24 ! color depth 24bits
bmp_info_header%biSizeImage = bmp%nx * bmp%ny * 3
open(9, file = fn//'.bmp', form = 'binary', status = 'unknown')
write(9) bmp_file_header
write(9) bmp_info_header
write(9) bmp%rgb
close(9)
return
end subroutine wr_bmp
! convert to t_RGB
pure elemental type(t_rgb) function to_rgb(ir, ig, ib)
integer, intent(in) :: ir, ig, ib
to_rgb = t_rgb(achar(ib), achar(ig), achar(ir))
end function to_rgb
subroutine pr_bmp(dtv, unit, iotype, vlist, io, iomsg)
class(t_bmp(*, *)), intent(in) :: dtv
integer, intent(in) :: unit
character(len = *), intent(in) :: iotype
integer, intent(in) :: vlist(:)
integer, intent(out) :: io
character(len = *), intent(in out) :: iomsg
character(len = 30) :: fmt
if (iotype == 'LISTDIRECTED') then
write(unit, *, iostat = io) 'nx =', dtv%nx, ', ny =', dtv%ny
end if
end subroutine pr_bmp
end module m_bmp
module m_mandel
implicit none
integer, parameter :: maxiter = 254
contains
pure elemental integer function mandel(c)
complex, intent(in) :: c
complex :: z
z = (0.0, 0.0)
do mandel = 0, maxiter
z = z * z + c
if (abs(z) > 2.0) exit
end do
end function mandel
end module m_mandel
program mandel_main
use m_bmp
use m_mandel
implicit none
integer, parameter :: nx = 640, ny = 640
real , parameter :: x0 = -2.0, x1 = 2.0, y0 = -2.0, y1 = 2.0
integer :: ix, iy, iter, mandelbrot(nx, ny)
real :: x(nx), y(ny)
complex :: c(nx, ny)
type(t_bmp(nx, ny)) :: bmp
!
forall (ix = 1:nx) x(ix) = x0 + (x1 - x0) * (ix - 1) / real(nx - 1)
forall (iy = 1:ny) y(iy) = y0 + (y1 - y0) * (iy - 1) / real(ny - 1)
forall (ix = 1:nx, iy = 1:ny) c(ix, iy) = cmplx(x(ix), y(iy))
mandelbrot = mandel(c)
!
bmp%rgb = to_rgb(255 - mandelbrot, 255 - mandelbrot, 255 - mandelbrot)
call bmp%wr('test')
print *, 'BMP size: ', bmp
stop
end program mandel_main
オブジェクト指向は形式面からのみ見れば、派生型の拡張として、その中にその派生型を引数とする副プログラムが含まれたものと見なせます。人為的な派生型の場合、それを引数としてとる副プログラムは汎用性のない専用のものとなるので、データ構造と副プログラムが一体化されるのは自然なことだと考えられます。これは行列やベクトルのような汎用性のある自然なデータ構造の場合、データ構造と副プログラムが切り離されている方が自然なことと対照されると思います。
BMP ファイルは構造が単純で、定義を調べて定義の通りに Fortran の派生型として変数を並べてゆけば何とかなります。ここで Fortran は実行速度の最適化のために変数の並べ替えを行うので、それを抑止するために sequence 文によって定義通りに変数を並べさせます。
BMP ファイルとしては 24bit の RGB 形式を選ぶことにします。この時、一般の画像サイズに対しては、データに適切な padding をしてやる必要があるので、画像サイズを 4 の倍数に限定して使用することで、これらの処理を省略することにします。本質的ではないので、プログラム上ではチェックのようなものはしないことにします。なお画像サイズは、Fortran2003 で導入されたパラメータ付派生型(parameterized derived type) を利用して宣言時に指定することにします。
BMP ファイルの書き出しは、オブジェクト指向的に派生型に含まれたサブルーチンで行うことにします。バイナリファイルとして派生型を出力すれば、中身の変数が望み通りそのまま書き出されます。
派生型変数を出力するとき、いちいち要素を書き並べるのが面倒になります。ユーザー定義派生型I/O を用いると派生型変数名を I/O 文に書いた時の挙動を制御できます。ここでは並び入出力形式 (List directed I/O) での出力用のサブルーチンを定義して、画像サイズを書かせるようにしています。
実行結果
コンソール
BMP size: nx = 640 , ny = 640
続行するには何かキーを押してください . . .
BMP ファイル
###Fortran 2008 CoArray
ここでは Fortran 2008 で導入された並列プログラミング用の仕組み CoArray を用いて、Mandelbrot 集合の BMP 図を並列に生成してみます。上記のプログラムを基に CoArray 版に改造します。
比較的小さな改変で CoArray による並列プログラムに直せました。
CoArray については別記事を参照してください。
http://qiita.com/cure_honey/items/7694a2824b56de16040e
http://qiita.com/cure_honey/items/3af353390fe913ae86a9
プログラム
module m_bmp
implicit none
type :: t_bmp_file_header
sequence
integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
integer(4) :: bfSize ! file size in bytes
integer(2) :: bfReserved1 = 0 ! always 0
integer(2) :: bfReserved2 = 0 ! always 0
integer(4) :: bfOffBits
end type t_bmp_file_header
!
type :: t_bmp_info_header
sequence
integer(4) :: biSize = Z'28' ! size of bmp_info_header ; 40bytes
integer(4) :: biWidth
integer(4) :: biHeight
integer(2) :: biPlanes = 1 ! always 1
integer(2) :: biBitCount
integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
integer(4) :: biSizeImage
integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
integer(4) :: biYPelsPerMeter = 3780 ! 96dpi
integer(4) :: biClrUsed = 0
integer(4) :: biClrImportant = 0
end type t_bmp_info_header
!
type :: t_rgb
sequence
character :: b, g, r ! order is b g r
end type t_rgb
!
type :: t_bmp
type(t_rgb), allocatable :: rgb(:, :)[:] ! coarray
contains
procedure :: wr => wr_bmp_caf
end type
contains
subroutine wr_bmp_caf(bmp, fn)
class(t_bmp), intent(in) :: bmp
character(len = *), intent(in) :: fn
type(t_bmp_file_header) :: bmp_file_header
type(t_bmp_info_header) :: bmp_info_header
integer :: me, ni, ir, ic
ni = num_images() ! coarray
me = this_image() ! coarray
associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2) * ni)
bmp_file_header%bfSize = 14 + 40 + 0 + nx * ny * 3
bmp_file_header%bfOffBits = 14 + 40
bmp_info_header%biWidth = nx ! nx shouold be a multiple of 4
bmp_info_header%biHeight = ny
bmp_info_header%biBitCount = 24
bmp_info_header%biSizeImage = nx * ny * 3
end associate
if (me == 1) then
open(9, file = fn//'.bmp', access = 'stream', status = 'replace')
write(9) bmp_file_header
write(9) bmp_info_header
endfile(9)
close(9)
end if
if (me > 1) sync images(me - 1) ! order 1, 2, 3...., n ! coarray
open(9, file = fn//'.bmp', access = 'stream', status = 'old', position = 'append')
write(9) bmp%rgb ! coarray
close(9)
if (me < ni) sync images(me + 1) ! coarray
end subroutine wr_bmp_caf
! convert to t_RGB
pure elemental type(t_rgb) function to_rgb(ir, ig, ib)
integer, intent(in) :: ir, ig, ib
to_rgb = t_rgb(achar(ib), achar(ig), achar(ir))
end function to_rgb
end module m_bmp
module m_mandel
integer, parameter :: maxiter = 254
contains
pure elemental integer function imandel(c)
complex, intent(in) :: c
complex :: z
z = c
do imandel = 1, maxiter
if (abs(z) > 2.0) exit
z = z * z + c
end do
end function imandel
end module m_mandel
program Mandel
use m_bmp
use m_mandel
implicit none
integer, parameter :: nx = 1920, ny = 1920
real , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5
complex, allocatable :: c(:, :)[:], z(:, :)[:] ! coarray
integer :: i, j, me, ni, my
integer, allocatable :: niter(:, :)[:] ! coarray
type (t_bmp) :: bmp ! coarray
ni = num_images() ! coarray
me = this_image() ! coarray
if (mod(ny, ni) /= 0) stop 'ny must be a multiple of num_image()'
my = ny / ni
allocate(c(nx, my)[*], z(nx, my)[*], bmp%rgb(nx, my)[*]) ! coarray
!
! make 2D-mesh : c(:, :)
!
block
real :: x(nx), y(ny)
x = [( (x1 - x0) / (nx - 1) * (i - 1) + x0, i = 1, nx )]
y = [( (y1 - y0) / (ny - 1) * (i - 1) + y0, i = 1, ny )]
forall (i = 1:nx, j = 1:my) c(i, j) = cmplx(x(i), y((me - 1) * my + j)) ! coarray
end block
!
! main iteration : niter(:, :)
!
allocate(niter(nx, my)[*]) ! implicitly sync all ! coarray
niter = imandel(c) ! coarray
!
! make bmp file: Mandel.bmp
!
niter(:, 1) = 255 ! line
! coarray
bmp%rgb = to_rgb(255 - niter, 255 - niter, 255 - niter) ! coarray
call bmp%wr('mandel') ! coarray
deallocate(bmp%rgb) ! implicitly sync all ! coarray
end program Mandel
実行結果
ここでは 5 つの横長の短冊状に領域を分割して並列に計算しています。それぞれの領域の境界に線を引いてあります。
#あとがき
以上、昔風の FORTRAN プログラムから出発して、徐々に今時の Fortran に書き換えてゆくことで、今時の Fortran で可能になった機能の一部を紹介しました。また新機能の背後にある発想も述べることで入門としてみました。