CAF と MPI 例題比較

W. Gropp 等の Using MPI の最初の例題を CAF で書いて見ます。

つらつら本を読んでいますと、並列プログラミングの場合、1)通信、2)同期、3)排他を考えなければならないようです。通信がなければ単なる独立な実行と同じになります。通信する場合は順序（タイミング）が定まらねばなりませんが、順序が一意に確定していれば同期によって実現します。順序に任意性がある場合はさらに排他処理が必要になります。

プログラム

CAF の場合単方向通信なので、通信と同期が分離していて、明示的に同期を取らなければなりません。MPI の場合基本がブロッキング型の双方向通信なので、通信が暗黙に同期を含んでいます。

また Fortran 2008 の CAF の場合は集団通信用の命令がないため*1、ギャザー・スキャター型の、順序に任意性のある演算の場合、明示的に排他処理をしなければなりません。 MPI は沢山の集団通信用のルーチンを用意して、排他処理を明示的にしなくてよいようにしてあるようです。

CAF

    program CAF3_2
      implicit none
      integer, parameter :: kd = kind(1.0d0)
      real(kd), parameter :: pi_kd = 4 * atan(1.0_kd)
      real(kd) :: pi[*], mypi, h, sum, x, f, a
      integer :: n[*], i, me, ne
      !
      f(a) = 4.0_kd / (1.0_kd + a * a)
      !
      me = this_image()
      ne = num_images()
      do 
        ! Broadcast n  
        if (me == 1) then 
          print *, 'enter the number of intervals: (0 quits) '
          read *, n
          pi = 0.0_kd
          sync images(*)
        else
          sync images(1)
          n = n[1]
        end if
        if (n == 0) exit
        !
        h = 1.0_kd / n
        sum = 0.0_kd
        do i = me, n, ne
          x = h * (real(i, kd) - 0.5_kd)
          sum = sum + f(x)
        end do
        mypi = h * sum
        !
        ! reduction
        critical 
          pi[1] = pi[1] + mypi
        end critical
        sync all
        if (me == 1) then
          print *, 'pi is', pi, ' Error is', abs(pi - pi_kd)
        end if  
      end do
    end program CAF3_2

CAF 実行結果

 enter the number of intervals: (0 quits)
10
 pi is   3.14242598500110       Error is  8.333314113051493E-004
 enter the number of intervals: (0 quits)
100
 pi is   3.14160098692312       Error is  8.333333331833614E-006
 enter the number of intervals: (0 quits)
1000
 pi is   3.14159273692313       Error is  8.333333401111531E-008
 enter the number of intervals: (0 quits)
10000
 pi is   3.14159265442313       Error is  8.333342904620622E-010
 enter the number of intervals: (0 quits)
0

MPI

    program mpi3_2
      use mpi
      implicit none
      integer, parameter :: kd = kind(0.0d0)
      real(kd), parameter :: pi25dt = 3.141592653589793238462643_kd
      real(kd) :: mypi, pi, h, sum, x, f, a
      integer :: i, n, ierr, myid, numprocs
      ! statement function
      f(a) = 4.0_kd / (1.0_kd + a * a)
      !
      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)      
      do 
        if (myid == 0) then
          print *, 'Enter the number of intervals: (0 quits) '
          read  *, n
        end if    
      ! broadcast n
        call MPI_BCAST(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
      !
        if (n == 0) exit
      !
        h = 1.0_kd / n
        sum = 0.0_kd
        do i = myid + 1, n, numprocs
          x = h * (real(i, kd) - 0.5_kd) 
          sum = sum + f(x)
        end do  
        mypi = h * sum
      !  collect all the partial sums  
        call MPI_REDUCE(mypi, pi, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
      !
        if (myid == 0) then
          print *, 'pi is ', pi, ' Error is', abs(pi - pi25dt)
        end if
      end do  
      call MPI_FINALIZE(ierr)
    end program mpi3_2

MPI 実行結果

C:\Users\O\Documents\Visual Studio 2015\Projects\MPI\MPI\x64\Debug>mpiexec -n 4 mpi.exe
 Enter the number of intervals: (0 quits)
10
 pi is    3.14242598500110       Error is  8.333314113051493E-004
 Enter the number of intervals: (0 quits)
100
 pi is    3.14160098692312       Error is  8.333333331833614E-006
 Enter the number of intervals: (0 quits)
1000
 pi is    3.14159273692313       Error is  8.333333312293689E-008
 Enter the number of intervals: (0 quits)
10000
 pi is    3.14159265442312       Error is  8.333307377483834E-010
 Enter the number of intervals: (0 quits)
0

*1:Fortran 2015 で導入予定です。

CAF の勉強

Fortran 2008 で導入された CAF (CoArray Fortran) 関係の命令を少し勉強したので、まとめをメモしておきます。coarray は PGAS 言語に分類される並列プログラミング用の仕組みです。ここで扱うのは Fortran 2008 の範囲で、Fortran 2015 で拡張される予定の部分については扱いません。勉強中なので間違っていたらゴメン！

コードを適切にセグメントに分けると、その部分集合 {Pi} はシリアル実行なら全順序集合としての構造を持つ。

CAF による並列実行の時、各イメージが独立も実行され相互にやり取りをしないなら、それぞれが全順序集合のコードセグメント {Pi}, {Qi}, {Ri}･･･ のように分かれる。各イメージごとに、CoArray を通じたデータのやり取りをする場合、データの読み込み・書き込みは排他的に順序づけて行なわなければならない。つまりコードセグメント間に順序関係を考えねばならない。

異なるイメージのコードセグメント間に順序関係を考えるとき、コードセグメント全体としては、半順序の構造をもつ。CAF でのデータのやり取りは、順序が定義されているコードセグメント間でしか許されない。（半順序なので、任意の二つのコードセグメント (Pi, Qj) を取り出したとき、この二つに順序が定義されているとは限らない。

CAF でコードセグメント間に順序を定義するときは、sync 命令を用いる。

一つのコードセグメントが複数のコードセグメントとデータをやり取りする時、データをやり取りする二つのコードセグメント間に順序関係が定義されても、データを直接やり取りしないコードセグメント間には順序が定義されていない。この場合、これらの間の順序は任意であるので、順序に依らずにデータ読み書きの排他性を保つ必要がある。この目的で critical..end critical、lock..unlock、atomic の三つの命令がある。

critical と lock はコード実行に関する排他性を保証する機能であり、atomic はデータへのアクセスの排他性を保証する機能である。critical は全イメージがブロッキング型で排他実行される。lock はより詳細に排他実行されるイメージを制御でき、ノンブロッキング型の排他実行も可能である。atomic 命令は、変数へのアクセスを排他的に行うための命令群であり、ノンブロッキング型の排他実行になっている。

d.hatena.ne.jp

critical	no flag	blocking	random order
lock	lock_type	blocking	random order
lock+acquired	lock_type+logical	non-blocking	random order
atom	atomic_logical_kind/atomic_int_kind	non-blocking	user-specified order

プログラム例

独立なイメージの実行　（イメージ間の順序構造なし）

ソース・プログラム

    program CAF00
      implicit none
      print *, 'Hello from', this_image(), 'out of', num_images(), 'images.' 
    end program CAF00

実行結果

 Hello from           3 out of           8 images.
 Hello from           6 out of           8 images.
 Hello from           8 out of           8 images.
 Hello from           1 out of           8 images.
 Hello from           7 out of           8 images.
 Hello from           2 out of           8 images.
 Hello from           5 out of           8 images.
 Hello from           4 out of           8 images.
続行するには何かキーを押してください . . .

イメージの実行　（環状の順序構造：　切れ目が必要）

ソース・プログラム

    program CAF04
      implicit none
      integer :: m[*], itmp[*]
      integer :: i, n, k
      n = num_images()
      i = this_image()
      m = i
      itmp = m
      print *, i, m
      sync all
      
      k = i - 1
      if (k == 0) k = n
      m[k] = itmp

      if (i == 1) print *
      sync all
      print *, i, m
    end program CAF04

実行結果

 3 3
 1 1
 2 2
 4 4

 3 4
 1 2
 2 3
 4 1
続行するには何かキーを押してください . . .

順序を完全な円環にすると、身動きが取れなくなるので、切れ目を入れて直線と端点処理にしている。

Critical による排他実行

ソース・プログラム

    program CAF05
      implicit none
      integer :: me, ni
      integer :: k[*] = 1
      integer :: p(4)[*] 
      
      me = this_image()
      ni = num_images()
      sync all

      critical
        p(k[1])[1] = me
        k[1] = k[1] + 1
      end critical    
      
      sync all
      if (me == 1) print *, p
    end program CAF05

実行結果

           1           3           4           2
続行するには何かキーを押してください . . .

lock による排他実行

ソース・プログラム

    program CAF08
      use, intrinsic :: iso_fortran_env, only : lock_type
      implicit none
      type(lock_type) :: list_lock[*]  
      integer :: i, p[*] = 0
      p = this_image()
      sync all
      if (this_image() == 1) then
        do i = 2, num_images()
           lock(list_lock[1])
             p[1] = p[1] + p[i]
           unlock(list_lock[1])
        end do
        print *, this_image(), num_images(), p
      end if
    end program CAF08    

実行結果

           1          10          55
続行するには何かキーを押してください . . .

atomic による変数への排他アクセス

ソース・プログラム

  program CAF10
    use, intrinsic :: iso_fortran_env
    implicit none
    logical(atomic_logical_kind) :: locked[*] = .true.
    logical :: val
    integer :: me, p = 1, q = 2, i = 0
    me = this_image()
    if (me == p) then
      sync memory
      call atomic_define(locked[q],.false.)
    else if (me == q) then
      val = .true.
      do while (val)
        i = i + 1  
        call atomic_ref(val,locked)
      end do
      sync memory
    end if
    print *, me , i
  end program CAF10    

実行結果

           1           0
           2         188
続行するには何かキーを押してください . . .

イメージ１が初めに lock に排他的にアクセスしているあいだ、イメージ２は実行可能になるまで、188回ループを回って順番を待っている。（ノンブロッキング実行のため待たされることは無い）

Modern Fortran 入門記事 （Intel の広告）

“Fortran, the greatest of the programming languages!”
www.youtube.com

このような紹介記事がありました。入門としてのクオリティーは･･･お、おぅ．．．
Be Curious and Learn: Modern Fortran in a Nutshell | InfoWorld

PGAS でよーがす

www.hpcwire.com

HPC wire に PGAS (Partitioned Global Address Space) がそろそろ来るという記事がありました。Fortran 2008 で導入された coarray Fortran が PGAS 言語の一つとして出てきています。

PGASの定義は、以下の要に与えられています。

PGAS defined

PGAS programming models offer a partitioned global shared memory capability, via a programming language or API, whether special support exists in hardware or not.

来そうな理由として以下のものが挙げられています。

Factor 1: Hardware support for more and more cores connected coherently.

Factor 2: Low latency interconnects.

Factor 3: Software support growing. 

【乞食速報！】Springer 大安売り

Springer　大安売り

www.springer.com

恒例の Springer 大安売りです。一部の電子書籍ですが $9.99 ないし 9,99 EUR になっています。古い Fortran 本など買えます。六月十四日まで。たぶん１会計にしか使えないクーポンコードなので、複数買いは一括でやる必要があります。

秋にもセールはあると思いますが。

Guide to Fortran 2008 Programming　ハードカバーが安い！

なぜか６割引き。

小ネタ

Fortran でテトリス

github.com

カーソル操作に Linux コンソール用のライブラリを利用しているようなので、汎用ではありません。

bash on windows のコンソールで動きました。

日立、メインフレームのハード開発から撤退 ＩＢＭから調達

www.nikkei.com

　日立製作所はメインフレームのハードウエア開発から撤退すると発表した。今後は米ＩＢＭから供給を受ける。ＩＢＭ製のハードに日立のメインフレーム向け基本ソフト（ＯＳ）を搭載したメインフレーム製品の提供を2018年度内に始める。

　メインフレームはメガバンクや企業、行政機関などの基幹業務システムに広く採用されていたが、近年はパソコンサーバーによるオープンシステムの台頭により市場規模が縮小している。1990年代半ばに１兆円規模だったメインフレームの国内出荷金額は15年度には400億円前後まで落ち込んでいる。

　日立は事業の選択と集中に向けて、ハード開発からの撤退を決めた。ただ既存システムの安定稼働を重視する企業を中心にメインフレームの需要は根強い。既存システムのソフトウエア資産の活用を支援するため、メインフレーム向けＯＳやミドルウエアソフトの開発は今後も継続する。

メインフレームの国内出荷金額は15年度には400億円前後まで落ち込んでいる。

浮動小数点数では結合則が成り立たない

IEEE754 の規格では、加減乗除の二項演算について、実数とその浮動小数点数表現の構造が保たれることを要求しています。つまり、加算の場合を例にとると、

a + b = c

fl(a) + fl(b)= fl(a + b) = fl(c)

この性質があるので、実数の可換な性質は保たれます。（overflow とかしなければ）

fl(a) + fl(b)= fl(a + b) = fl(b + a) = fl(b) + fl(a)

一方、実数で成り立つもっと基本的な演算規則である結合則の方は、三項演算に相当していて、規格からの要請もなく（たぶんｗ）成り立たない場合があります。「桁落ち」や「積み残し」が問題化するのは結合則からのズレに還元できると思います。結合則が成り立たない場合、演算の順序を明示的に指示する必要が出てきます。

(a + b) + c = a + (b + c)

(fl(a) + fl(b)) + fl(c) = fl(a + b) + fl(c) 

/=

fl(a) + (fl(b) + fl(c)) = fl(a) + fl(b + c) 

プログラムで、見てみることにしてみます。

プログラム

乱数で実数を生成して、交換子と結合子を計算してして、交換則と結合則からのずれを見てみます。ズレがなければ 0.0 になるはずです。

    program Console1
      implicit none
      integer, parameter :: kd = kind(1.0d0), nmax = 100
      real(kd) :: x(2, nmax), y(3, nmax), z(nmax)
      integer :: i
      call random_seed()
    
      print *, 'commutator'
      call random_number(x)  
      z = comm_plus(x(1, :), x(2, :))   
      print *, 'non-commutative =', count(z /= 0.0_kd), '/', nmax

      print *, 'associator' 
      call random_number(y)  
      z = assoc_plus(y(1, :), y(2, :), y(3, :))
      print *, 'non-associative =', count(z /= 0.0_kd), '/', nmax
      print '(e30.15)', pack(z, z /= 0.0_kd)   
      
      do i = 1, nmax
        if (z(i) /= 0.0_kd) print *,  y(:, i)
      end do  
      
    contains
      real(kd) pure elemental function comm_plus(x, y)
        real(kd), intent(in) :: x, y
        comm_plus = (x + y) - (y + x)
      end function comm_plus
      
      real(kd) pure elemental function comm_mult(x, y)
        real(kd), intent(in) :: x, y
        comm_mult = (x * y) - (y * x)
      end function comm_mult
      
      real(kd) pure elemental function assoc_plus(x, y, z)
        real(kd), intent(in) :: x, y, z
        assoc_plus = ((x + y) + z) - (x + (y + z))
      end function assoc_plus 

      real(kd) pure elemental function assoc_mult(x, y, z)
        real(kd), intent(in) :: x, y, z
        assoc_mult = ((x * y) * z) - (x * (y * z))
      end function assoc_mult 
    end program Console1

計算結果

計算結果から分かるように、交換則の方は確かに満たされていますが、結合則の方は２割程度の確率で満たされていないことが分かります。

ステパノフの本や関数型言語の本などでは、結合則を根底において、モノイド、半群、群などの造を導入して、結合則による演算順序の組み換えを縦横に使って並列動作などさせていますが、整数や文字列のような離散値の場合はよいとして、浮動小数点数の場合は結合則が崩れているので、実数の演算の問題にどこまで適用できるのか慎重に見極める必要がある気がします。

 commutator
 non-commutative =           0 /         100
 associator
 non-associative =          22 /         100
        -0.222044604925031E-15
         0.222044604925031E-15
        -0.444089209850063E-15
         0.222044604925031E-15
         0.222044604925031E-15
         0.222044604925031E-15
        -0.222044604925031E-15
         0.222044604925031E-15
         0.222044604925031E-15
         0.222044604925031E-15
        -0.444089209850063E-15
        -0.222044604925031E-15
         0.222044604925031E-15
        -0.222044604925031E-15
         0.111022302462516E-15
         0.444089209850063E-15
         0.111022302462516E-15
         0.111022302462516E-15
         0.222044604925031E-15
         0.222044604925031E-15
         0.222044604925031E-15
         0.444089209850063E-15
  0.197336647554121       0.583104599529825       0.230515619550751
  0.478456421601024       0.365995662803589       0.695310195955152
  0.765030199209025       0.821438970892834       0.822794756357349
  0.154500948792594       0.776532299353388       0.343291349327081
  0.123168650301795       0.892064457677992       0.624841959267652
  5.340595195978214E-002  0.883762457929461       0.282097598062034
  0.271274877739308       0.990276136982007       0.295042576304925
  0.385678503561146       0.550016479450929       0.368699724012741
  0.649285046006193       0.772531681072521       0.291322963201651
  6.278512502868444E-002  0.484351337035122       0.846490729577704
  0.791263186492663       0.730073920942973       0.590455112135355
  0.599791134233700       0.716285730192570       0.379474641874127
  0.271981838214423       0.686661556999306       0.655546444804149
  9.427993838386348E-002  0.387223376386792       0.590407433074261
  0.295995343550855       0.100678724962143       0.578481472642591
  0.631307369405927       0.834670683810024       0.762999289601546
  0.634957778719910       0.129886627681722       0.170474343695826
  0.131887838808105       6.074681699437982E-002  0.581256123449080
  8.217122777577210E-002  0.437050020391703       0.846554402242005
  0.522359744366527       0.748392434610684       0.567776991641634
  0.555136565205923       0.518809881107340       0.750619304274506
  0.487565338352254       0.577165852794002       0.970361542180519
続行するには何かキーを押してください . . .