今時の Fortran 入門 ( Introduction to Modern Fortran )

(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 で可能になった機能の一部を紹介しました。また新機能の背後にある発想も述べることで入門としてみました。