Topics
    ソリューション
    Menu

    技術コラム:MPI並列プログラミング入門

    目次

    1. はじめに
    2. 逐次プログラムの解析と並列化指針の決定
    3. MPIの初期設定
    4. 入出力の取り扱い
    5. エラー/終了処理
    6. 並列計算のためのプロセス分割
    7. データの送受信
    8. 並列性能測定
    9. 参考文献

    1. はじめに

    ここでは、MPI(Message Passing Interface)を用いて、ハウスコードなどユーザ自身が作成した逐次計算プログラムを並列化する手順について述べる。プログラミング言語はFortran言語で解説しているが、CやC++言語で書かれたプログラムの場合も、同様の考え方で行える。

    一般に、プログラムを並列化する目的としては、

    1. 計算の高速化
      多数のCPUを用いて並列に計算を進め、大規模計算に要する時間を短縮する。
    2. 大容量メモリの確保
       計算に必要なメモリ空間を多数のPCに分散し、大容量のメモリ領域を確保する。

    の2点が挙げられるが、ここでは、ひとまずメモリの分散確保については考えず、並列計算を行う全CPUにおいてそれぞれ逐次版プログラムと同じだけのメモリ領域を確保し、計算の高速化のみを図る。並列計算がうまく行けば、メモリ節約(メモリ分散)の手だてを考えることは、それほど難しくはない。
    MPIは、MPIフォーラムにより開発・規格化されたCPU間のデータ通信(メッセージパッシングという)を行うためのライブラリであり、FortranやCなどのプログラム中からサブルーチン(関数)呼び出しを行うことによって、並列計算に必要なデータ交換を実現する。したがって、MPIを用いて作成された並列プログラムは、PCクラスタはもちろん、SMP機など様々なプラットホームで利用できる移植性の高いといえる。MPIで提供されているサブルーチンは計100を超えるが、ここでは、

    • プログラムの初期化・終了サブルーチン
      MPI_INIT(), MPI_COMM_SIZE(), MPI_COMM_RANK(), MPI_ABORT(), MPI_FINALIZE()
    • 並列計算を行う全CPUに同じデータを送信するブロードキャスト通信
      MPI_BCAST()
    • 特定の1対1のCPU間における通信
      MPI_SEND(), MPI_RECV()
    • 並列性能を計測するための時間計測サブルーチン
      MPI_WTIME()

    に話題を絞り、基本的な(最低限の)並列化プログラムの開発方法について述べる。以上のサブルーチンの使用法を理解するだけで、基本的なMPI並列プログラムを開発することができる。
    典型的な(逐次)数値計算プログラムの構造は、たいていの場合、逐次計算プログラムを起動後、

    1. キーボードやあらかじめ用意されたファイルから、計算に必要なデータ読み込み
    2. 計算実行
    3. 計算結果を画面やファイルへ出力

    という手順となる。この流れの通り、MPIを用いて並列化(3つのCPUで実行)するには、図1のようプログラム構成となる。次節以降、この流れに沿って解説を進める。

    mpi_parallel_programming_guide01

    図1 MPI並列プログラムの流れ

    2. 逐次プログラムの解析と並列化指針の決定

    既存のプログラムを並列化する場合、逐次プログラムのどの部分が支配的な計算なのかを知る必要がある。並列化を行う際、全ての処理部分を並列化することはできないため、プログラム中で最も支配的な計算部分を見つけ出して並列化する必要がある。この部分を特定するには、LinuxやUnixのprofコマンドやgprofコマンド等が便利である。これは、プログラム実行時のサブルーチンごとの計算時間割合を表示する。

    gprofコマンドの使用方法

    コンパイルオプションとして-pg(profコマンドの場合には-p)をつけてコンパイルを行う。プログラムを実行すると、自動的にgmon.outファイルが生成される。その状態で

    $ gprof 実行ファイル名

    を実行すると、図2のように、サブルーチンごとにCPU時間、呼び出し回数などの情報が出力される。

    CPU時間
    の比率
    CPU時間
    の累計
    CPU時間サブルーチン・関数名
    %comulativeself selftotal 
    timesecondssecondscallss/calls/callname
    73.24348.41348.3180100.040.04axcal_
    10.97400.5852.17152.1752.17angle2_
    4.05419.8519.27   cvtas_a_to_t
    3.02434.2414.3981.801.80updatephi_
    1.61441.927.6880.960.96updater_
    0.86446.024.11   rs_get_field
    0.76449.653.63   cvt_text_to_ieee_t
    0.75453.233.57   for_read_seq_lis_xmit
    0.71456.623.39122880.000.00rotmat2_
    0.40458.511.90   for__ignore_space
    0.34460.121.61   cvtas_t_to_a
    0.33461.691.56   for__desc_ret_item
    0.25462.891.21   for__cvt_value
    0.25464.071.1890.130.13updatechi_
    0.23465.171.09   for__get_su
    0.22466.191.03310680.000.00vecinp_
    0.20467.170.97   rs_find_field


    図2 gprofコマンドの実行結果

    図2の例では、サブルーチンaxcalの計算部分が総計算時間の70%以上を占めており、計8,000回以上も呼び出されていることが確認できる。このサブルーチン内の計算、あるいはサブルーチン呼び出しを並列化すれば、計算の高速化をはかることができる。

    3.MPIの初期設定

    MPIプログラムは、同一のプログラム実行ファイルが並列計算に用いる全てのCPU上でそれぞれ実行される。したがって、どのCPUに何の計算部分を担当させるのか、どのCPUからどのCPUへ何のデータを転送するのかを明示的に記述する。
    MPIプログラムを初期化するための関数呼び出しについて説明する。

    mpi_parallel_programming_guide03

    図3 MPI初期設定部分

    並列化を行うための初期設定、環境設定として、MPIのサブルーチンを呼び出すプログラムの宣言部分に

    include ‘mpif’ integer istatus(MPI_STATUS_SIZE)

    を指定する。MPIを初期化するため、図3のように、はじめに次の3つのMPIルーチンを呼び出す。

    call MPI_INIT(ierr) call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr) call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)

    MPI_INI()はMPI初期化のためのサブルーチンであり、MPI_COMM_SIZE()とMPI_COMM_RANK()は、このプログラムにおいて並列計算に使用しているCPU数(プロセス数)をnprocsに、自分のCPU番号(ランク値)をmyrankにそれぞれ格納する。並列計算のCPU数は、MPIプログラムの起動時(mpirunコマンドの引数)に設定され、プログラム中でその値を取得し変数へ格納する必要がある。CPUには0, 1, 2, …とそれぞれ番号(ランク)が割り振られ、この番号を使って通信の相手先、計算内容などを特定する。
    図4に実際の並列プログラムにおけるMPIの初期設定部分を示す。3行目ではMPI用インクルードファイルmpif.hを指定している。mpif.hの部分は必ず小文字で記述しなければならない。4行目は、MPIルーチンを呼び出すときに必要な情報格納配列 istatus を宣言している。この2つは、プログラムの先頭(宣言部分、ステートメントの前)に記述する。
    13行目において、並列計算に用いるCPUの総数NPROCSと、自分のCPU番号MYRANKは、並列化作業において重要であり使用頻度が高いため、どのサブルーチンからでも参照できるようにCOMMON変数として設定している。NPROCSとMYRANKは、サブルーチンMPI_COMM_SIZE()、MPI_COMM_RANK()をそれぞれ呼び出すことで得られる。
    15行目は、MPI_INIT()を呼び出し、MPI環境の初期化処理を行っている。このMPIルーチンは、他のどのMPIルーチンよりも前にコールされなければならない。 16、17行目では、 NPROCS、MYRANKの値を取得するためMPI_COMM_SIZE()とMPI_COMM_RANK()をコールしている。

    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

          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    C
          include ‘mpif.h’
          integer istatus(MPI_STATUS_SIZE)
    C
          INTEGER,ALLOCATABLE,DIMENSION(:) ::
         &      NN, NE, NPD, NFM
          REAL*8,ALLOCATABLE,DIMENSION(:) ::
         &      A, B, U, X, TMP
    C
          REAL*8 time1, time2
    C
          COMMON /MPIVAL/ MYRANK, NPROCS
    C
          call MPI_INIT(ierr)
          call MPI_COMM_SIZE(MPI_COMM_WORLD, NPROCS, ierr)
          call MPI_COMM_RANK(MPI_COMM_WORLD, MYRANK, ierr)
    C
          nums=NPROCS-1
          if(myrank.eq.0) then
            write(6,*) ‘MPI: number of slaves =’,nums
          end if
    C
    C——————-
    C START COMPUTING
    C——————-



    図4 MPI初期設定のプログラム例

     図5に示すように、3つのCPUを用いて並列計算を行う場合、MPI_COMM_SIZE()を呼び出すと、どのCPUにおいても変数NPROCSに3が格納されて戻る。また、MPI_COMM_RANK()をコールすると、MYRANKには0~2のCPU番号(ランク値)がそれぞれ割り当てられて戻ってくる。この各CPUに振られたランク値(MYRANK)を用いて、自分あるいは相手のCPUを特定し並列計算の制御を行う。

    mpi_parallel_programming_guide05

    図5 合計3 CPUでMPIプログラムを起動したの場合の制御変数値

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18

          SUBROUTINE GAUSS(A, B, U, X, TMP, NP, NPDF, NU, NN)
    C
    C——————————————-
          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    C
          include ‘mpif.h’
          integer istatus(MPI_STATUS_SIZE)
    C
          DIMENSION A(NU),B(NPDF,3),U(NPDF,3,X(NPDF<3)
    C
          COMMON /MPIVAL/ MYRANK, NPROCS
    C——————————————-
          EPSI=1.0D-6
    C
    C
          DO NITE=1,1000 !ITERATION NUMBER
    C
             DO I=1,3



    図6 サブルーチンプログラムにおいて必要な設定

    図6に示すように並列化を行う際には関連するサブルーチンにおいても6行目、7行目のようにインクルードファイルmpif.hの指定と、istatus()配列の設定を行うとともに、11行目のようにCOMMON変数を設定し、MYRANK、NPROCSの参照が可能なようにしておく。

    4. 入出力の取り扱い

    計算に必要なデータの入力は、OPEN文で指定されたファイル入力と、キーボード(標準)入力がある。このうちファイル入力は、NFS(Network File System)などによって全てのCPUが同一のファイルへアクセス可能になっていれば、それぞれのCPUが入力ファイルに接続し、データを読み込むことができるため、図7に示すように、変更の必要はない。

    open (unit=10,file=’Initial_data.dat’)

    read (10,*) data1, data2

    図7 ファイル入力(変更の必要なし)

    mpi_parallel_programming_guide08

    図8 入出力部分

    計算結果を出力ファイルに書き込む場合、複数のCPUが同時に同じファイルへ書き出しを行うと、ファイル内容の整合性がとれなくなるなどの問題が発生するため、図8に示すように、並列計算を行うCPU群のうち1台のみ(これをマスタと呼ぶ)がファイル出力を行う。たとえば、マスタのCPUをランク値=0とすると、出力ファイルをOPENし、WRITEする部分のプログラムmyrank=0の時だけ実行するよう図9のように変更する。

    if (myrank .eq. 0) then   !myrank=0 の時のみ以下の処理が実行される

      open(unit=20,file=’Solution_data.dat’)

      write(20,*) answer1, answer2

    end if

    図9 マスタ(myrank=0)のみがファイル出力するプログラム

    標準入出力(キーボード入力、ディスプレイ出力)においても、ファイル出力と同様にCPU群の1台(マスタ)のみが行うように指定しなければならない。さもなければ、CPUの数に応じて同じデータ入力を繰り返さなければならず、出力も重複してしまう。普通、キーボードから入力された変数値は全てのCPUにおいて必要となるため、入力後、変数の値をブロードキャスト通信し、全てのCPUに転送する。以上の操作のプログラム例を図10に示す。

          if (myrank .eq. 0) then   !myrank=0 の時のみ以下の処理が実行される

            write(*,*) ‘Please in put Nodal Points number.’

            read(*,*) nn

          endif

    C   変数 nn を他の全ての CPU に送信する(ブロードキャスト通信)

          call (MPI_BCAST(nn,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

    図10 キーボード入力とブロードキャスト通信

    画面出力およびキーボード入力をマスタのみで行うように指定し、キーボードから入力された値(変数nn)をブロードキャスト通信MPI_BCAST()する。
    ブロードキャスト通信は、図11に示すように、ある1つのCPU(プロセス)から並列計算に用いる全CPUに同じデータ(メッセージ)を一斉送信する。先のプログラム例は、マスタ(myrank=0)において入力された変数nnをほかの全CPUに送信している。

    mpi_parallel_programming_guide11

    図11 ブロードキャスト通信

    MPI_BCAST()の使用方法

    call MPI_BCAST(message, count, datatype, root, comm, ierror)

    入出力変数

    message 送信元(マスタ)では送信バッファ(変数、あるいは配列)の先頭アドレスを指定する。受信側では受信バッファの先頭アドレスを指定する。
    ierror(整数型)完了コードが戻る。

    入力変数

    count(整数型)送信または受信するメッセージの長さ(データの個数)を指定する。
    datatype(整数型)送信または受信するデータの型(表 1参照)を指定する。
    root(整数型)送信元のランク値を指定する。
    comm(整数型)MPI_COMM_WORLDと指定する。

    表1 MPIのデータ型

     FORTRANのデータ型MPIの基本データ型名称
    整数型INTEGERMPI_INTEGER
    実数型REALMPI_REAL
    DOUBLE PRECISIONMPI_DOUBLE_PRECISION
    複素数型COMPLEXMPI_COMPLEX

    実際のプログラムには、READ文やWRITE文は数多くあるため、grepコマンドやエディタの検索機能を用いて一括検索すると便利である。

    grepコマンド

    指定した文字列を検索する。

    $ grep -in  検索する文字列 ファイル名(ワイルドカードを使って複数指定可)

    オプション -i は文字列の大文字・小文字を無視し、-nはマッチした文字列の行番号を表示させる。

    mpi_parallel_programming_guide12

    図12 grepコマンドの実行例

    5. エラー/終了処理

    mpi_parallel_programming_guide13

    図13 MPI終了処理部分

    プログラムの終了時、MPI環境の終了処理を行うサブルーチン、MPI_FINALIZE()をコールしなければならない。また、プログラムの途中で強制的に実行を中止させたい場合には、MPI_ABORT()をコールする。全てのCPUでこのサブルーチンがコールされなければ、プログラムが正常に終了できなくなってしまう。
    実際のプログラムにおいては、エラー処理などのためSTOP文が複数存在する。そこで、図14のサブルーチンSTOP_PROGRAMを作成し、STOP文と置き換えると便利である。実際には、図15のように、プログラムの最後にサブルーチンSTOP_PROGRAMを呼び出して処理を終了させれば良い。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13

          SUBROUTINE STOP_PROGRAM
    C
          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    C
          include ‘mpif.h’
    C
          COMMON /MPIVAL/ MYRANK,NUMS
    C
          WRITE(*,*) ‘STOP MYRANK=’,MYRANK
          CALL MPI_FINALIZE(IERR)
    C
          STOP
          END

    図14 MPI終了処理のサブルーチン

    1

    2
    3
    4
    5



          WRITE(*,*) ‘All calculations successfully finished!’
    C
          CALL STOP_PROGRAM
    C
          END

    図15 終了処理サブルーチンの呼び出し例

    6. 並列計算のためのプロセス分割

    mpi_parallel_programming_guide16

    図16 並列計算部分

    DOループを分割し、並列計算する場合について、説明する。例として、係数マトリックスが全て同一である3つの連立1次方程式

    [A]{x1}={b1}
    [A]{x2}={b2}
    [A]{x3}={b3}

    を解いて3つの解を求める場合を考える。逐次プログラムではベクトル{b1}, {b2}, {b3}を2次元配列 b に格納し、図17のプログラムのように、DOループにより3つの連立1次方程式を解いている。ここで、{b1}は2次元配列bのb(*,1)部分に、{b2}はb(*,2)部分、{b3}はb(*,3)部分にそれぞれ格納さている。また、計算後の解{x1}, {x2}, {x3}は2次元配列b(*,3)に格納される。

    1
    2
    3
    4
    5
    9

          do i = 1,3
    C
            call CGSOLV         !連立1次方程式を解く
    C
         &       (A,b(1,i),g,g1,w,u(1,1),tmp,npdf,nite,nsol,i)
    C
          enddo

    図17 並列処理前のプログラム

    ここで、並列化の方法として、図19に示すように3つの連立1次方程式をそれぞれCPUに割り当て計3 CPUで並列計算を行う場合について考える。なお、全てのCPUは係数マトリックス[A]と右辺ベクトル{b1}, {b2}, {b3}の格納された全ての配列を保有しているものとする。

    mpi_parallel_programming_guide18

    図18 並列処理のイメージ

    各CPUの番号(ランク値)と、右辺ベクトルが格納された2次元配列bとを関連づける。つまり、myrankが0のCPUは[A]{x1}={b1}を計算をするために、右辺ベクトル{b1}が格納されたb(*,1)を引数としてサブルーチンCGSOLVに渡す。myrankが1ならば、右辺ベクトル{b2}が格納されたb(*,2)を、myrankが3なら右辺ベクトル{b3}が格納されたb(*,3)を引数に指定する。この方法により、図19のように変更される。

    1
    2
    3
    4
    5
    6
    10
    11

          do i = 1,3
            if (myrank .eq. i-1) then   !各 CPU で異なる作業を行なう
    c
              call CGSOLV               !連立1次方程式を解く
         &        (A,b(1,i),g,g1,w,u(1,1),tmp,npdf,nite,nsol,i)
    c
            endif
          enddo

    図19 並列処理のプログラム

    7. データの送受信

    mpi_parallel_programming_guide20

    図20 データ通信部分

    係数マトリックスが全て同じ3つの連立1次方程式

    [A]{x1}={b1}
    [A]{x2}={b2}
    [A]{x3}={b3}

    を計3 CPUで解く例より、データの送受信について説明する。
    各CPUが担当する連立1次方程式が解き終わると、得られた解を1つのCPUに集約しなければならない。各CPUは解を格納するため2次元配列bを持っており、例えばCPU 2番 (myrank=1)においては、求められた解はb(*,2)に格納されている。
    CPU 1 (myrank=0)にCPU 2、CPU 3で得られた解{x2}, {x3}を集約するため、1対1ブロッキング通信サブルーチンであるMPI_SEND()とMPI_RECV()を用いる。MPI_SEND()はデータを送信するサブルーチン、MPI_RECV()は送信されたデータを受け取るサブルーチンである。{x2}, {x3}を受け取るCPU 1 (myrank=0)はMPI_RECV()をコールし、{x2}, {x3}を引き渡すCPU 2 (myrank=1)、CPU 3 (myrank=2)はMPI_SEND()をコールする。この通信のイメージを図21に示す。

    mpi_parallel_programming_guide21

    図21 1対1MPI通信のイメージ

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18

          if (myrank .eq. 0) then
    c
            do i = 1,2
              call MPI_RECV(b(1,i+1), npdf, MPI_DOUBLE_PRECISION,
         &                i, 1, MPI_COMM_WORLD, ISTATUS, ierr)
    c
            end do
    c
          else if (myrank .eq. 1) then
    c
            call MPI_SEND(b(1,2), npdf, MPI_DOUBLE_PRECISION, 0 ,1,
         &               MPI_COMM_WORLD, ierr)
    c
          else if (myrank .eq. 2) then
    c
            call MPI_SEND(b1,3), npdf, MPI_DOUBLE_PRECISION, 0 ,1,
         &               MPI_COMM_WORLD, ierr)
    c
          end if

    図22 1対1 MPI通信のプログラム

    図22に実際のプログラムを示す。CPU 2、CPU 3はそれぞれ配列b(*,2)、b(*,3)をCPU 1 (myrank=0)に送信する。そして、CPU 1ではDOループにより送信先のランク値を指定して受信し、配列Xに解を格納する。この通信により3つの連立1次方程式の解はCPU1に全て集約することができる。

    以下にMPI_SEND()とMPI_RECV()の使用方法を示す。

    MPI_SEND()

    call MPI_SEND(message count, datatype, dest, tag, comm, ierror)

    出力

    ierror(整数型)完了コードが戻る。

    入力

    message 送信バッファの先頭アドレスを指定する。
    count(整数型)送信するメッセージの長さを指定。
    datatype(整数型)送信するデータのデータ型を指定する。(データ型については表 1を参照)
    dest(整数型)送信先のCPUのランク値(myrank)を指定する。宛先プロセスに対してメッセージを送信しない場合、MPI_PROC_NULLを指定。
    tag(整数型)送信するメッセージの種類を区別する場合に指定。特に区別しない場合、適当な数値を入れる。
    comm(整数型)特別の場合を除いて、MPI_COMM_WORLDと指定すればよい。

    MPI_RECV()

    call MPI_RECV(message, count, datatype, source, tag, comm, status, ierror)

    出力

    message 受信バッファ(配列データ)の先頭アドレスを指定する。
    ierror(整数型)完了コードが戻る。

    入力

    count(整数型)受信するメッセージの長さを指定する。
    datatype(整数型)受信するデータのデータ型を指定する。(データ型については表1を参照)
    source(整数型)受信したいメッセージを送信するプロセスのcomm内でのランクを指定する。
    tag(整数型)受信したいメッセージに付けられているタグの値を指定する。
    comm(整数型)特別の場合を除いて、MPI_COMM_WORLDと指定すれば良い。
    status 状況オブジェクトの配列を指定する。プログラムの初めに指定した

    integer istatus(MPI_STATUS_SIZE)

    のistatusを使用する。

    以上でMPI並列プログラムの完成である。

    8. 並列性能測定

    並列性能を測定するための実行時間の計測方法について述べる。
    実行時間の測定方法のため、現在の時刻を返す関数MPI_WTIME()を使用し、経過時間を計測したい処理部分の前後でMPI_WTIME()を実行し、得られた2つの時間差がその部分の経過時間となる。

    MPI_WTIME()の使用方法

    MPI_WTITE()を呼び出すと、時間(単位は秒)戻る。戻り値は倍精度実数である。

    ATAI =MPI_WTIME()

    図23にMPI_WTIME()を使用して実行時間を計測するプログラム例を示す。計算の前後で得られた時間の差が計測したい部分の経過時間となる。

          READ*8 START, END
    C
          START = MPI_WRITE()         !時間を測定

    測定する計算部分

          END = MPI_WRITE()           !時間を測定
    C
          TIME = END – START          !差により経過時間がわかる

    図23 MPI_WTIMEの使用例

    9. 参考文献

    1. 青山幸也、並列プログラミング入門(MPI編)、高度情報科学技術研究機構(RIST)神戸センター(2013)。
    2. 青山幸也、チューニング技法虎の巻、高度情報科学技術研究機構(RIST)神戸センター(2013)。
      上記2つは、http://www.hpci-office.jp/pages/seminar_textより入手可能である。
    3. P.パチェコ、(秋葉訳)、MPIプログラミング、培風館、(2001)。
    4. 樫山和男、西村直志、牛島省、並列計算法入門、丸善、(2003)。
    5. 湯浅太一、安村通晃、中田登志之、はじめての並列プログラミング、共立出版、(1998)。
    6. 飯塚肇、緑川博子、並列プログラミング入門、丸善、(2000)。
    7. 石川裕、佐藤三久、堀敦史、住元真司、原田浩、高橋俊行、Linuxで並列計算をしよう、共立出版、(2002)。