例題1: 1次元拡散方程式

問題の説明とシリアルFortranプログラム

周期的な領域 $x \in [- 1,1],$ において、1次元拡散方程式

$\frac{\partial f}{\partial t} = D \frac{\partial^{2} f}{\partial x^{2}}$ ,

を考えます。 $D$ は拡散係数。有限差分法を用いて以下のように表します。

$\frac{f_{i}^{n + 1} - f_{i}^{n}}{Δt} = D \frac{f_{i + 1}^{n} - 2 f_{i}^{n} + f_{i - 1}^{n}}{{(Δx)}^{2}}$ .

$Δt$ は時間刻み幅、 $Δx$ は空間格子の刻み幅です。時間、空間はそれぞれ $t_{n} = n Δt + t_{0}, x_{i} = i Δx + x_{0}$ と離散化されます。ここで $t_{0}$ , $x_{0}$ は、起点を表しますが、後の議論に関係しませんので特に気にする必要はありません。 $f_{i}^{n}$ は $n$ ステップ目の時刻の、 $i$ 番目の格子点における $f$ の値を表します。 $f_{i}^{n + 1}$ について解くと

$f_{i}^{n + 1} = f_{i}^{n} + \frac{D Δt}{{(Δx)}^{2}} (f_{i + 1}^{n} - 2 f_{i}^{n} + f_{i - 1}^{n})$

が得られます。これは時間ステップ $n$ における $f$ の値を用いて、次の時間ステップ $n + 1$ における $f$ の値が計算できることを表しています。この1次元拡散方程式の有限差分解を求めるプログラムをFortran 90によって作成したものが以下です。
(サンプルソースコードをダウンロード: diffuse_serial.f90)

     1	program diffusion_serial

     2	  implicit none

     3	  integer :: NX0=512               ! grid points [input]
     4	  integer :: NT=100000             ! iterations [input]
     5	  character (len=100) :: DATAFILE  ! data file name [input]
     6	  integer, parameter :: NBX=1      ! boundary grid points
     7	  integer, parameter :: NOUT=1000  ! output size
     8	  real, parameter :: DC=1.         ! diffusion coefficient
     9	  real, parameter :: LENGTH=2.     ! domain size [-LENGTH/2:LENGTH/2]
    10	  integer :: nx
    11	  integer :: ix,it
    12	  real, allocatable :: f(:),df(:),x(:)
    13	  real, allocatable :: fp(:),xp(:)
    14	  real :: dx
    15	  real :: t,dt
    16	  real :: pi
    17	  real :: fmin,fmax
    18	  real :: total,total0
    19	  character (len=100), parameter :: MINMAXFILE='minmax.dat'
    20	  real, parameter :: INFINITY=1.e15

    21	  namelist /input/ NX0,NT,DATAFILE
    22	  namelist /minmax/ fmin,fmax

    23	  ! input
    24	  read(5,input)
    25	  DATAFILE=trim(DATAFILE)

    26	  nx=NX0+2*NBX

    27	  allocate(f(nx),df(nx),x(nx))
    28	  allocate(fp(NX0+1),xp(NX0+1))

    29	  pi = 4.*atan(1.)
    30	  dx = LENGTH/NX0

    31	  do ix = 1, nx
    32	     x(ix) = dx*(ix-NBX-.5) - .5*LENGTH
    33	  end do
    34	  xp(1:NX0+1)= .5*(x(NBX:NX0+NBX)+x(NBX+1:NX0+NBX+1))

    35	  dt = .5 * dx * dx / DC
    36	  t = 0.

    37	  !
    38	  ! initialize array
    39	  !
    40	  f(1:nx)=0.
    41	  df(1:nx)=0.

    42	  !
    43	  ! initial condition
    44	  !
    45	  f(1:nx) = .5*(cos(pi*x(1:nx))+1.)
    46	  fp(1:NX0+1)= .5*(f(NBX:NX0+NBX)+f(NBX+1:NX0+NBX+1))

    47	  fmin= INFINITY
    48	  fmax=-INFINITY
    49	  do ix=1,NX0+1
    50	     fmin=min(fmin,fp(ix))
    51	     fmax=max(fmax,fp(ix))
    52	  end do

    53	  !
    54	  ! write initial data
    55	  !
    56	  open(10,file=DATAFILE,form='unformatted')
    57	  write(10) xp
    58	  write(10) t,fp

    59	  total=sum(f(1+NBX:NX0+NBX))*dx
    60	  total0=total
    61	  write(6,'(3(a,f12.4),a)') &
    62	       & 'time = ',t, &
    63	       & ', total = ',total, &
    64	       & ', difference = ',(total0-total)/total0*100,' [%]'

    65	  !
    66	  ! main routine
    67	  !
    68	  do it = 1, NT
    69	     
    70	     df(1+NBX:NX0+NBX) = DC*dt/dx**2 * &
    71	          & ( f(2+NBX:NX0+NBX+1) - 2. * f(1+NBX:NX0+NBX) + f(NBX:NX0+NBX-1) )
    72	     f(1+NBX:NX0+NBX) = f(1+NBX:NX0+NBX) + df(1+NBX:NX0+NBX)

    73	     f(        1:        NBX)=f(NX0+    1:NX0+    NBX)
    74	     f(NX0+NBX+1:NX0+NBX+NBX)=f(    NBX+1:    NBX+NBX)

    75	     !
    76	     ! output data
    77	     !
    78	     t = t + dt
    79	     if(mod(it,NOUT).eq.0) then
    80	        fp(1:NX0+1)= .5*(f(NBX:NX0+NBX)+f(NBX+1:NX0+NBX+1))
    81	        do ix=1,NX0+1
    82	           fmin=min(fmin,fp(ix))
    83	           fmax=max(fmax,fp(ix))
    84	        end do
    85	        write(10) t,fp

    86	        total=sum(f(1+NBX:NX0+NBX))*dx
    87	        write(6,'(3(a,f12.4),a)') &
    88	             & 'time = ',t, &
    89	             & ', total = ',total, &
    90	             & ', difference = ',(total0-total)/total0*100.,' [%]'
    91	     end if

    92	  end do
    93	  close(10)

    94	  open(20,file=MINMAXFILE)
    95	  write(20,minmax)
    96	  close(20)

    97	  stop

    98	end program diffusion_serial

図1にプログラム中における空間格子の配置を示します。領域の大きさは $L$ でプログラム中でLENGTHで与えられ、座標の起点は $x_{0} = - L / 2 - Δx / 2$ となります。領域内、外の格子点はそれぞれNX0, NBXで与えられ、全格子点数は $nx = NX0 + 2 NBX$ となります。領域外の格子点は、境界での値を計算するために必要になります。例えば、領域内左端の格子(1+NBX)の値を計算するには、両隣 (NBXとNBX+2)での値が必要になります。この場合NBXは1で十分ですが、高次の差分では、領域外に1つ以上の格子点が必要になる場合がありますので、図1では、より一般的な場合を示しています。ここで示したような格子点配置では、周期境界条件は(1:NBX) [(NX0+NBX+1:NX0+NBX+NBX)]の格子点のデータを (NX0+1:NX0+NBX) [(NBX+1:NBX+NBX)]の格子点のデータと同一視することで満たされます。

図1 格子点配置

プログラム中、全質量

$F = \int_{- L / 2}^{L / 2} f d x$ ,

も計算されていますが、これは周期境界条件下では保存量であり、時間的に変化しません。保存量 $F$ が正しく保存しているかどうかは、数値解が正しいかどうかのひとつの目安となります。以下はこのプログラムの標準出力です。

# ./diffuse_serial < input.diffuse
time =       0.0000, total =       1.0000, difference =       0.0000 [%]
time =       0.1221, total =       1.0000, difference =       0.0000 [%]
time =       0.2441, total =       1.0000, difference =       0.0000 [%]
time =       0.3662, total =       1.0000, difference =       0.0000 [%]
time =       0.4883, total =       1.0000, difference =       0.0000 [%]
time =       0.6104, total =       1.0000, difference =       0.0000 [%]

"difference [%]"は保存量の初期値との差異を出力していますが、解は正しく求められているようです。図2に解の時間発展を示します。

Animation of Time Evolution of the Solution of the Diffusion Equation

図2 拡散方程式の解の時間発展

並列化

並列化とはプログラム中の互いに依存しない処理を異なるプロセッサ上に分散配置し、複数のプロセッサで同時に計算させることで、理想的には実行時間を $1 / (プロセッサ数)$ に減らすことができます(実際は、プロセッサ間のデータの通信、非均等な処理の分散、並列化できない処理の存在などにより、理想的な並列化は不可能ですが)。ここで考えている問題においては、 $f$ を求める計算は、両隣の格子点における $f$ の値しか必要としないため並列化可能です。そこで、 $f$ をいくつかのプロセッサに分散することを考えてみましょう。 $N_{p}$ をプロセッサ数、 $N_{x}$ を領域内の格子点数とします。 $I_{p}$ 番目( $I_{p} = 0, \dots, N_{p} - 1$ )のプロセッサは、以下で与えられる $N_{min}$ から $N_{max}$ までの格子点における計算を受け持ちます。

$N_{min} = int (\frac{N_{x}}{N_{p}}) I_{p} + 1$ , $N_{max} = max (int (\frac{N_{x}}{N_{p}}) (I_{p} + 1), N_{x})$ ,

ここで $max(a, b)$ は $a$ , $b$ のうち大きい方を返す関数、 $int(a)$ は小数点以下を切り捨てた $a$ の整数部を返す関数とします。次の図は、データがどのように分散されるかを示しています。

図3 複数プロセッサへのデータの分散

$f_{i}^{n + 1}$ を計算するためには、両隣の格子点での値 $f_{i + 1}^{n}$ と $f_{i - 1}^{n}$ が必要になりますが、各プロセッサは受け持っている格子点の両端の格子点での値を計算するために必要なデータを持っていません。例えば $N_{x} = 12$ , $N_{p} = 3$ としましょう。1番目のプロセッサは1番目から4番目の格子を受け持ちますが、4番目の格子点の $f$ を計算するために必要な5番目の格子のデータは持っていません。周期境界条件では、1番目のプロセッサは格子1の $f$ を計算するために必要な格子0と同一視される格子12のデータも持っていないことになります。同様にしてプロセッサ2は格子4,9のデータ、プロセッサ3は格子8,1(=13)のデータが、各時間ステップにおいて必要になります。

図4 プロセッサ間のデータ通信

このようなプロセッサ間のデータの通信はMPIライブラリを用いて実装することになります。以下は、MPIライブラリを用いた並列化例です。
(サンプルソースコードのダウンロード: diffuse_mpi.F90),

     1	program diffusion_mpi

     2	  use mpi

     3	  implicit none

     4	  integer :: NX0=512               ! grid points [input]
     5	  integer :: NT=100000             ! iterations [input]
     6	  character (len=100) :: DATAFILE  ! data file name [input]
     7	  integer, parameter :: NBX=1      ! boundary grid points
     8	  integer, parameter :: NOUT=1000  ! output size
     9	  real, parameter :: DC=1.         ! diffusion coefficient
    10	  real, parameter :: LENGTH=2.     ! domain size [-LENGTH/2:LENGTH/2]
    11	  integer :: nx
    12	  integer :: ix,it
    13	  real, allocatable :: f(:),df(:),x(:)
    14	  real, allocatable :: fp(:),xp(:)
    15	  real :: dx
    16	  real :: t,dt
    17	  real :: pi
    18	  real :: fmin,fmax
    19	  real :: total,total0
    20	  character (len=100), parameter :: MINMAXFILE='minmax.dat'
    21	  real, parameter :: INFINITY=1.e15

    22	  integer, parameter :: rootproc=0
    23	  integer :: myid=rootproc,nps=1,ip,nrange
    24	  integer, allocatable :: npmin(:),npmax(:),npbmin(:),npbmax(:)
    25	  integer :: ierr
    26	  integer :: p_send,p_recv
    27	  integer :: mpirealtype,tag=0,status(mpi_status_size)

    28	  namelist /input/ NX0,NT,DATAFILE
    29	  namelist /minmax/ fmin,fmax

    30	  call mpi_init(ierr)
    31	  call mpi_comm_rank(mpi_comm_world,myid,ierr)
    32	  call mpi_comm_size(mpi_comm_world,nps,ierr)

    33	  mpirealtype=mpi_real
    34	#ifdef DBLE
    35	  mpirealtype=mpi_double_precision
    36	#endif

    37	  allocate(npmin(0:nps-1),npmax(0:nps-1))
    38	  allocate(npbmin(0:nps-1),npbmax(0:nps-1))

    39	  ! input
    40	  if(myid.eq.rootproc) read(5,input)
    41	  call mpi_bcast(NX0,       1,mpi_integer,  rootproc,mpi_comm_world,ierr)
    42	  call mpi_bcast(NT,        1,mpi_integer,  rootproc,mpi_comm_world,ierr)
    43	  call mpi_bcast(DATAFILE,100,mpi_character,rootproc,mpi_comm_world,ierr)
    44	  DATAFILE=trim(DATAFILE)

    45	  nrange=NX0/nps
    46	  do ip=0,nps-1
    47	     npmin(ip)=nrange*ip+1  +NBX
    48	     npmax(ip)=nrange*(ip+1)+NBX
    49	  end do
    50	  npmax(nps-1)=NX0+NBX
    51	  npbmin(0:nps-1)=npmin(0:nps-1)-NBX
    52	  npbmax(0:nps-1)=npmax(0:nps-1)+NBX

    53	  nx=NX0+2*NBX

    54	  allocate(f(nx),df(nx),x(nx))
    55	  allocate(fp(NX0+1),xp(NX0+1))

    56	  pi = 4.*atan(1.)
    57	  dx = LENGTH/NX0

    58	  do ix = 1, nx
    59	     x(ix) = dx*(ix-NBX-.5) - .5*LENGTH
    60	  end do
    61	  xp(1:NX0+1)= .5*(x(NBX:NX0+NBX)+x(NBX+1:NX0+NBX+1))

    62	  dt = .5 * dx * dx / DC
    63	  t = 0.

    64	  !
    65	  ! initialize array
    66	  !
    67	  f(1:nx)=0.
    68	  df(1:nx)=0.

    69	  !
    70	  ! initial condition
    71	  !
    72	  f(1:nx) = .5*(cos(pi*x(1:nx))+1.)
    73	  fp(1:NX0+1)= .5*(f(NBX:NX0+NBX)+f(NBX+1:NX0+NBX+1))

    74	  fmin= INFINITY
    75	  fmax=-INFINITY
    76	  do ix=1,NX0+1
    77	     fmin=min(fmin,fp(ix))
    78	     fmax=max(fmax,fp(ix))
    79	  end do

    80	  !
    81	  ! write initial data
    82	  !
    83	  if(myid.eq.rootproc) then
    84	     open(10,file=DATAFILE,form='unformatted')
    85	     write(10) xp
    86	     write(10) t,fp

    87	     total=sum(f(1+NBX:NX0+NBX))*dx
    88	     total0=total
    89	     write(6,'(3(a,f12.4),a)') &
    90	          & 'time = ',t, &
    91	          & ', total = ',total, &
    92	          & ', difference = ',(total0-total)/total0*100,' [%]'
    93	  endif

    94	  !
    95	  ! main routine
    96	  !
    97	  do it = 1, NT

    98	     df(npmin(myid):npmax(myid)) = DC*dt/dx**2 * &
    99	          & (    f(npmin(myid)+1:npmax(myid)+1) &
   100	          & - 2.*f(npmin(myid)  :npmax(myid)) &
   101	          & +    f(npmin(myid)-1:npmax(myid)-1) )
   102	     f(npmin(myid):npmax(myid)) = &
   103	          & f(npmin(myid):npmax(myid)) + df(npmin(myid):npmax(myid))

   104	     if(nps.eq.1) then

   105	        f(npbmin(myid):npmin(myid)-1)=f(npmax(myid)-NBX+1:npbmax(myid))
   106	        f(npmax(myid)+1:npbmax(myid))=f(npmin(myid):npmin(myid)+NBX-1)

   107	     else

   108	        !
   109	        ! Periodic boundary
   110	        ! Processor 0 sends f(NBX+1:NBX+NBX) and 
   111	        ! Processor nps-1 receives/stores it in f(NX0+NBX+1:NX0+NBX+NBX)
   112	        !
   113	        p_send=mpi_proc_null
   114	        p_recv=mpi_proc_null
   115	        if(myid.eq.0)     p_send=nps-1
   116	        if(myid.eq.nps-1) p_recv=0
   117	        call mpi_send(f(1+NBX),NBX,mpirealtype, &
   118	             & p_send,tag,mpi_comm_world,ierr)
   119	        call mpi_recv(f(NX0+NBX+1),NBX,mpirealtype, &
   120	             & p_recv,tag,mpi_comm_world,status,ierr)

   121	        !
   122	        ! Periodic boundary
   123	        ! Processor nps-1 sends f(NX0+1:NX0+NBX) and
   124	        ! Processor 0 receives/stores it in f(1:NBX)
   125	        !
   126	        p_send=mpi_proc_null
   127	        p_recv=mpi_proc_null
   128	        if(myid.eq.0)     p_recv=nps-1
   129	        if(myid.eq.nps-1) p_send=0
   130	        call mpi_send(f(NX0+1),NBX,mpirealtype, &
   131	             & p_send,tag,mpi_comm_world,ierr)
   132	        call mpi_recv(f(1),NBX,mpirealtype, &
   133	             & p_recv,tag,mpi_comm_world,status,ierr)

   134	        !
   135	        ! Shift communications
   136	        ! send/receive boundary data to right/from left
   137	        !
   138	        p_send=myid+1
   139	        p_recv=myid-1
   140	        if(myid.eq.0)     p_recv=mpi_proc_null
   141	        if(myid.eq.nps-1) p_send=mpi_proc_null
   142	        call mpi_send(f(npmax(myid)-NBX+1),NBX,mpirealtype, &
   143	             & p_send,tag,mpi_comm_world,ierr)
   144	        call mpi_recv(f(npmin(myid)-NBX),NBX,mpirealtype, &
   145	             & p_recv,tag,mpi_comm_world,status,ierr)

   146	        !
   147	        ! Shift communications
   148	        ! send/receive boundary data to left/from right
   149	        !
   150	        p_send=myid-1
   151	        p_recv=myid+1
   152	        if(myid.eq.0)     p_send=mpi_proc_null
   153	        if(myid.eq.nps-1) p_recv=mpi_proc_null
   154	        call mpi_send(f(npmin(myid)),NBX,mpirealtype, &
   155	             & p_send,tag,mpi_comm_world,ierr)
   156	        call mpi_recv(f(npmax(myid)+1),NBX,mpirealtype, &
   157	             & p_recv,tag,mpi_comm_world,status,ierr)

   158	     endif

   159	     !
   160	     ! output data
   161	     !
   162	     t = t + dt
   163	     if(mod(it,NOUT).eq.0) then

   164	        if(nps.gt.1) then
   165	           call mpi_gather( &
   166	                & f(npmin(myid)),npmax(myid)-npmin(myid)+1,mpirealtype, &
   167	                & f(1+NBX),   npmax(myid)-npmin(myid)+1,mpirealtype, &
   168	                & rootproc,mpi_comm_world,ierr)
   169	           p_send=mpi_proc_null
   170	           p_recv=mpi_proc_null
   171	           if(myid.eq.rootproc) p_recv=nps-1
   172	           if(myid.eq.nps-1   ) p_send=rootproc
   173	           call mpi_send(f(NX0+NBX+1),NBX,mpirealtype, &
   174	                & p_send,tag,mpi_comm_world,ierr)
   175	           call mpi_recv(f(NX0+NBX+1),NBX,mpirealtype, &
   176	                & p_recv,tag,mpi_comm_world,status,ierr)
   177	        endif

   178	        if(myid.eq.rootproc) then
   179	           fp(1:NX0+1)= .5*(f(NBX:NX0+NBX)+f(NBX+1:NX0+NBX+1))
   180	           do ix=1,NX0+1
   181	              fmin=min(fmin,fp(ix))
   182	              fmax=max(fmax,fp(ix))
   183	           end do

   184	           write(10) t,fp
   185	           total=sum(f(1+NBX:NX0+NBX))*dx
   186	           write(6,'(3(a,f12.4),a)') &
   187	                & 'time = ',t, &
   188	                & ', total = ',total, &
   189	                & ', difference = ',(total0-total)/total0*100.,' [%]'
   190	        endif
   191	     end if

   192	  end do
   193	  close(10)

   194	  if(myid.eq.rootproc) then
   195	     open(20,file=MINMAXFILE)
   196	     write(20,minmax)
   197	     close(20)
   198	  endif

   199	  call mpi_finalize(ierr)

   200	  stop

   201	end program diffusion_mpi

MPIソースコードの変更点を見てみましょう(重要箇所は赤で示してあります)。まず、第1にMPIによって並列化することを"宣言"しなければなりません(use mpi)。MPIモジュールによってMPIによって提供されるサブルーチンや変数が使用可能になります。MPIソースコードのうち並列化されている箇所はサブルーチンmpi_init, mpi_finalizeで囲まれ、それらの間で複数のプロセッサを使用することができます。並列化された箇所で、各プロセッサがまずするべきことは、自分のプロセッサ番号(ランク)[mpi_comm_rank]と全プロセッサ数[mpi_comm_size]を取得することです。プロセッサランクと全プロセッサ数から、45-52行目で自分が受け持つ格子点の範囲を計算します。 $f$ を計算する処理は、この範囲を用いて各プロセッサに分散されることになります。先に説明したように、各プロセッサは、必要なデータを取得するために両隣のプロセッサと通信する必要がありますが、この通信を行うサブルーチンがmpi_sendとmpi_recvです。サブルーチンmpi_gatherとその直後の送受信は全データを1つのプロセッサ上に集め、出力するための処理です。MPIプログラムでは、データの入出力を1つのプロセッサに押しつけておいた方が簡単で安全です。このため標準入力からの入力も1つのプロセッサだけが行い、読み込んだデータをmpi_bcastによって全プロセッサに送信する方法をとっています。

通信

1対1通信

mpi_send/mpi_recvによってデータ通信を行っている所が4ヶ所あります。周期境界条件のため両端のプロセッサがデータを送受信する所と、各プロセッサが右隣、左隣に受けもっているデータの境界のデータを送受信するシフト通信する所それぞれ2ヶ所づつです。それぞれMPIサブルーチンをコールする前に変数p_send/p_recvにデータを送受信するプロセッサランクを定義しています。最初のブロック(113-120行)では、プロセッサ0はf(NBX+1:NBX+NBX)を送信し、プロセッサnps-1はプロセッサ0からのデータを受け取りf(NX0+NBX+1:NX0+NBX+NBX)に保存します。他のプロセッサは何もしません。mpi_send/mpi_recv引数の送受信先プロセッサランクに"mpi_proc_null"を与えると、これらのサブルーチンは何もすることなく、ただちに返ります。

集団通信

mpi_gatherはデータを1つのプロセッサに集め、mpi_bcastはデータを全プロセッサに送信します。これらのサブルーチンはコミュニケータ(mpi_comm_world)で定義されたプロセッサのグループ(通常すべてのプロセッサ)間で集団的にデータ通信を行います。
Fortranユーザへの注意
C言語版MPIライブラリでは、各サブルーチンは返り値としてエラーコードを返します。Fortranでは、サブルーチンは返り値を持ちませんので、エラーコードはサブルーチンの最後の引数として与えられます。よって、C言語版に比べてFortran版では、サブルーチンの引数が1つ多く必要で、最後に必ずエラーコードのための整数の引数を取ります。通常、エラーが生じるとプログラムはただちに異常終了しますので、このエラーコードを参照することはありません。
mpirealtype
整数変数"mpirealtype"には、実際の浮動小数点精度を与えるMPI型(単精度ならmpi_real、倍精度ならmpi_double_precision)が代入されています。このような変数を用意している理由は、コンパイラに精度を変換する(-r8など)オプションを渡すことによって精度を変換する場合があるためです。仮に、全てのMPIサブルーチンにデータ型としてmpi_realやmpi_double_precisionなどと与えていた場合、コンパイラオプションによる精度の変換によって、全てのデータ型を書き換える必要があり、非常に繁雑です。よって、mpirealtypeを定義しておき、プリプロセッサによってどちらの精度を用いるか選択するようにしています。以下のコマンドは倍精度のプログラムを生成します。
```
# cpp -DDBLE diffuse_mpi.F90
# mpif90 -r8 -o diffuse_mpi diffuse_mpi.f90
```
mpi_comm_world
mpi_comm_worldとは、互いに通信するプロセッサのグループを含むコミュニケータのことです。タグのような意味で使われますが、通常コミュニケータを必要とする所ではmpi_comm_worldと書いておけば十分です。

デッドロック

周期境界条件のための通信を行っている箇所(113-120, 126-133行)は、シフト通信(138-145, 150-157行)と、p_send/p_recvをうまく設定すればまとめることができそうです。141-142行目をp_recv=nps-1/p_send=0, 152-153行目をp_recv=0/p_send=nps-1とすれば、108-133行目を取り除くことができます。しかし、このようなプログラムは、残念ながら、あるシステムにおいては上手く動作しません。mpi_sendはブロッキング通信を行うサブルーチンで、送信するデータをシステムバッファにコピーするまで返りません。図5上図に示すように、システムが十分大きなバッファを持つ場合、プログラムは正しく動作しますが、下図のようにバッファサイズが十分でない場合、mpi_sendは終了せず、プログラムはそれ以降に進めないことになります。このように、すべてのプロセッサが待ちの状態になってしますことを、"デッドロック"と言います。

Figure of Data Flow of Point-to-Point Communication and Deadlock

図5 1対1通信におけるデータの流れとデッドロック

mpi_sendとmpi_recvの順序を入れ換えるとどうでしょうか?状況はさらにひどく、システムバッファのサイズに依存せず、全てのシステムで明らかにデッドロックとなります。なぜなら、mpi_recvもブロッキング通信を行うサブルーチンであり、受信するべきデータを受け取るまで返りませんので、すべてのプロセッサがmpi_recvをコールして時点でプログラムはストップします。

デッドロックを回避する方法がいくつかあります。送受信するプロセッサのうち一方がmpi_send をコールし、他方がmpi_recvをコールする場合、システムバッファが十分でない場合でも通信は成功します(図6)。

Figure of Data Flow in Send-Recv Communication

図6 Send-Recv通信におけるデータの流れ

つまり、偶数/奇数番プロセッサがmpi_send/mpi_recvを先にコールするように変更すればよいことになります。または、これを自動的に行うmpi_sendrecvというサブルーチンも用意されています。

デッドロックを避ける他の方法として、非ブロッキング通信を行うサブルーチンを用いる方法があります。mpi_isendとmpi_irecvは非ブロッキング送受信を行うサブルーチンで、ただちに返ります。ブロッキング通信を行うサブルーチンを用いる理由は、それらがデータの正しさを保証してくれるためであり、逆に、非ブロッキング通信ではデータへの不正なアクセスを許してしまいます。これはmpi_irecvの挙動を見れば一目瞭然です。次のようにmpi_irecv(recvbuf,...) をコールし、直後で受け取るはずのデータrecvbufを参照するとしましょう。

call mpi_irecv(recvbuf,....)
tmp=recvbuf*recvbuf

mpi_irecvがすでに recvbufを受け取ったという保証はありませんので、recvbufには正しく値が代入されていないかもしれません。同様なことがmpi_isendにも起こり得ます。mpi_isendの直後に送信すべきデータにアクセスすると、思いがけずそのデータを上書きしてしまうかもしれません。mpi_irecv/mpi_isendがデータを受信/送信するまで待つためのサブルーチンがmpi_waitです。非ブロッキング通信サブルーチンを用いて、送受信する部分を書き直すと以下の様になります。

call mpi_isend(sendbuf,...,ireq1,...)
call mpi_irecv(recvbuf,...,ireq2,...)
call mpi_wait(ireq1,...)
call mpi_wait(ireq2,...)

リダクション演算

全質量の計算( $f$ の総和)は上記ソースコードではプロセッサrootprocで実行されていますが、これは $f$ が出力のためにすでに当該プロセッサ上に集められているためです。出力の必要がなく、 $f$ の全データを保持しているプロセッサがない場合、総和を計算するためにわざわざ全データを1つのプロセッサ上に集める必要はありません。各プロセッサが、保持している $f$ の局所的な和を計算し、それらを最後に足しあわせれば総和を計算することができます。このような種類の演算を、リダクション演算といい、MPIライブラリではリダクション演算を行うサブルーチンmpi_reduceを提供しています。このサブルーチンを使えば、仮に $f$ のサイズが非常に大きい場合、全データを集める通信に比べ、通信量を大幅に減らすことができます。総和を計算する部分(185行)は次のようになります。

total0=sum(f(npmin(pid):npmax(pid)))*dx
call mpi_reduce(total0,total,1,mpirealtype,mpi_sum,rootproc,mpi_comm_world,ierr)

派生型

MPIプログラミングにおいて最も重要なことは、いかにデータ通信量を減らすかということです。MPIサブルーチンコールには時間がかかりますので、コールの回数を減らすことは効果的です。例にあげたソースコードでは41-43行目にあるmpi_bcastをひとつにまとめることができそうです。(実際には、この変更はまったく効果はありませんが、派生型の例を示すために取り上げています。)

MPIライブラリでは、送受信されるデータをグループ化する枠組が提供されています。派生型と呼ばれるデータのグループを表すデータ型をユーザが定義し、それを一般のMPIのデータ型と同様に用いることができます。派生型を定義するためには、グループ化する各要素のデータ型、位置および大きさ(配列のサイズ)が必要になります。

図7 データのメモリレイアウト

図7にデータのグループの例を示します。例では、データグループは、2つのreal型(a, b)および長さが2のcharacter型(c)で構成されています。各要素の位置は、1番目の要素(a)からの変位で表されますので、a,b,cの位置は、それぞれ0, 16, 24ということになります。以下のサブルーチンで例に示しているソースコードで用いられる派生型(newtype)を定義します。

     1	subroutine inputtype(int1,int2,char1,newtype)

     2	  use mpi

     3	  implicit none

     4	  integer, intent(in) :: int1,int2
     5	  character (len=100), intent(in) :: char1
     6	  integer, intent(out) :: newtype
     7	  integer, parameter :: count=3
     8	  integer :: i
     9	  integer :: addr_start
    10	  integer :: blen(count)
    11	  integer :: disp(count)
    12	  integer :: tlst(count)

    13	  integer :: ierr

    14	  blen=(/ 1,1,100 /)
    15	  
    16	  tlst=(/ mpi_integer,mpi_integer,mpi_character /)
    17	  
    18	  call mpi_address(int1, disp(1),ierr)
    19	  call mpi_address(int2, disp(2),ierr)
    20	  call mpi_address(char1,disp(3),ierr)
    21	  addr_start=disp(1)
    22	  do i=1,count
    23	     disp(i)=disp(i)-addr_start
    24	  end do

    25	  call mpi_type_struct(count,blen,disp,tlst,newtype,ierr)
    26	  call mpi_type_commit(newtype,ierr)
    27	  
    28	  return

    29	end subroutine inputtype

データのグループは3つの要素(int1=NX0, int2=NT, char1=DATAFILE)で構成され、それぞれの位置、データ型、大きさは配列"blen", "tlst", "disp"に代入されています。mpi_addressはデータの位置を与えます。与えられた情報をもとに、mpi_type_structが新しい派生型(newtype)を構築し、その新しい派生型はmpi_type_commitによりMPIに通知されたあとに使用可能になります。作成された新しい派生型を用いて、1回mpi_bcastを呼ぶことでデータのブロック(NX0, NT, DATAFILE)を送信するには以下のように変更すればよいことになります。

call mpi_bcast(NX0,1,newtype,rootproc,mpi_comm_world,ierr)

この操作でNX0だけでなく、(NX0,NT,DATAFILE)が送信されます。

ソースコードなどのセット

ダウンロード (diffuse.tgz)

内容物

diffuse/README.jp: 説明書(日本語)
diffuse/README: 説明書
diffuse/Makefile: プログラムをコンパイル、実行し、結果を可視化するためのMakefile
diffuse/diffuse_serial.f90: シリアルFortranソースコード
diffuse/diffuse_mpi.F90: 並列化されたFortranソースコード
diffuse/diffuse.F90: Fortranプリプロセッサを用いて、シリアルとパラレルを選べるようにしたハイブリッド版Fortranソースコード
diffuse/diffuse2.F90: ちょっとだけMPIプログラミングの発展的内容を含むようにdiffuse.F90を変更したハイブリッド版ソースコード
diffuse/diff_pgplot.F90: PGPLOTによって結果をプロットするプログラムのFortranソースコード
diffuse/diff_gnuplot.f90: 結果をGNUPLOTによってプロットするためにunformattedデータをformattedデータに変更するプログラムのFortranソースコード
diffuse/ps2anim: 得られたポストスクリプトファイルから、PSUtilsとImageMagickを用いて動画ファイルに変換するPerlスクリプト

遊び方

"make"と入力すればdiffuse.psが作成され結果を見ることができます。,

# make

動画に変換するまでを実行するには、"make anim"と入力してください。使用するプロセッサ数を指定するには

# make MPI=16

とします。その他、makeの動作を変更するマクロがいくつか用意されていますので、詳細はREADMEを参照してください。

Written by "Ryusuke NUMATA" <rnumata at umd.edu>
http://www.glue.umd.edu/~rnumata/
Last Last modified: $Date: 2008-01-23 16:15:00 -0500 (Wed, 23 Jan 2008)$