メッセージ通信ライブラリを用いたプログラムの並列化例と計算速度および計算精度の評価

田村 克浩, 稲富 雄一, 長嶋 雲兵

Return

1 はじめに

Linpackベンチマークによる計算速度のコンピュータ世界ランキングTOP500 Super Computer Site[1]のランキングホルダはほとんどが256個以上のCPUをもつ並列計算機であり、上位の計算機の中には１万近くのCPUをもつ計算機もみられる。10年前のスーパーコンピュータは、高速な単一のCPUまたは数台程度で構成されていたが、ハードウエア性能の限界と開発コスト、およびソフトウエア上の並列プログラミング環境の向上から数百以上の複数のCPUを持つ並列計算機に取って代わられた。
スーパーコンピュータのような特殊なコンピュータだけでなく大学等の個々の研究室においても、コストパフォーマンスに優れたコモディティハードウエアを組み合わせたコンピュータ、いわゆるDOS/Vパソコンを、複数台ネットワークで接続し分散メモリマルチプロセッサ型クラスタを形成し、安価なスーパーコンピュータとして利用することもしばしば行われている。この場合、コストが非常に重要になるためLinuxまたはFreeBSD等 をOSとして採用し、MPICH、Parallel Virtual Machine (PVM)などのフリーのメッセージ通信ライブラリを利用した並列プログラミングが行われている。
並列計算は、並列台数に応じた高速化が図れることが最大の利点であるが、上記分散メモリマルチプロセッサ型クラスタではこれに加え、1）一台ではメモリ容量不足で扱えない問題を扱うことができる、2）大規模なシステム構築が容易である、3）共有メモリ型並列計算機（Shared　Memory　Parallel（SMP）machine）では大きな問題になるタイミングに依存するバグを内包しにくい等の利点を有している一方で、1）プログラミング時に通信をすべてスケジュールしなければならない、2）通信遅延のオーバヘッドが非常に大きいという欠点を有している。そのため、プログラミング時のコーディングの良し悪しによって、実行速度が大きく変化する。[2]
本報告では、メッセージ通信ライブラリMPIをもちいて、簡単な並列プログラムの実行速度の高速化および計算精度の挙動をAlta Technology AltaCluster、Hitachi　SR8000、IBM RS/6000 SP、SGI Origin2000の計4種類の並列計算機で測定したので報告する。

2 プログラムとMPIを用いた並列化

2進数における循環小数である10進数の0.1を109回分加えるFortranのプログラムを2種類作成した。一つは、0.1を109回逐次加え、総和を求めるプログラム（Program 1 in Figure 1）、もう一つは104回の部分和の和を105回加えることで、最終的に0.1 の109回の総和を計算するプログラム（Program 2 in Figure 2）を作成した。この２つのプログラムをMPIを用いてSPMD（Single Program Multiple Data）モデルによる並列化を行った（それぞれFigure 3（program 1）およびFigure 4（program 2））。Figure 3とFigure 4には並列化のために加えた命令を太字で示した。計算は計算誤差が現れやすいようすべて単精度で行った。

      PROGRAM main
      IMPLICIT NONE
      INTEGER nmax
      PARAMETER(nmax=1000000000)
      INTEGER i
      REAL sum
      sum=0.0e0
      DO i=1,nmax
        sum=sum+0.1e0
      ENDDO
      WRITE(6,200) sum
  200 FORMAT('The answer is ',e15.8,'.')
      END
Figure 1. Sequential Summation

      PROGRAM main
      IMPLICIT NONE
      INTEGER nmax,nunitsize
      PARAMETER(nmax=1000000000,nunitsize=10000)
      INTEGER i,j,nunit,namari
      REAL sum,psum,asum
      sum=0.0e0
      nunit=INT(nmax/nunitsize)
      namari=MOD(nmax,nunitsize)
      DO i=1,nunit
        psum=0.0e0
        DO j=1,nunitsize
          psum=psum+0.1e0
        ENDDO
        sum=sum+psum
      ENDDO
      IF(namari/=0) THEN
        asum=0.0e0
        DO i=1,namari
          asum=asum+0.1e0
        ENDDO
        sum=sum+asum
      ENDIF
      WRITE(6,200) sum
  200 FORMAT('The answer is ',e15.8,'.')
      END
Figure 2. Summation with Partial Sum

      PROGRAM main
      IMPLICIT NONE
      INCLUDE 'mpif.h'
      INTEGER nmax
      PARAMETER(nmax=1000000000)
      INTEGER i,j,k,nunit,istart,iend,
     &        ierr,myrank,isize,
     &        istatus(MPI_STATUS_SIZE)
      REAL sum,buf
      CALL MPI_INIT(ierr)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,isize,ierr)
      nunit=INT(nmax/isize)
      istart=nunit*myrank+1
      iend=nunit*(myrank+1)
      sum=0.0e0
      IF(myrank==isize-1) THEN
        DO i=istart,nmax
          sum=sum+0.1e0
        ENDDO
      ELSE
        DO i=istart,iend
          sum=sum+0.1e0
        ENDDO
      ENDIF
      IF(myrank==0) THEN
        DO i=1,isize-1
          CALL MPI_RECV(buf,1,MPI_REAL,i,MPI_ANY_TAG,
     &                  MPI_COMM_WORLD,istatus,ierr)
          sum=sum+buf
        ENDDO
      ELSE
        buf=sum
        CALL MPI_SEND(sum,1,MPI_REAL,0,1,
     &                MPI_COMM_WORLD,ierr)
      ENDIF
      IF(myrank==0) THEN
        WRITE(6,200) sum
  200   FORMAT('The answer is',e15.8,'.')
      ENDIF
      CALL MPI_FINALIZE(ierr)
      END
Figure 3. Parallelized Program of Sequential Summation.(Program1)

      PROGRAM main
      IMPLICIT NONE
     INCLUDE 'mpif.h' 
      INTEGER nmax,nunitsize
      PARAMETER(nmax=1000000000,nunitsize=10000)
      INTEGER i,j,ierr,myrank,isize,
     &        nunit,nunitamari,namari,istart,
     &        iend,istatus(MPI_STATUS_SIZE)
      REAL sum,psum,asum,buf
      CALL MPI_INIT(ierr)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,isize,ierr)
      nunit=INT(INT(nmax/nunitsize)/isize)
      nunitamari=MOD(INT(nmax/nunitsize),isize)
      namari=MOD(nmax,nunitsize)
      istart=nunit*myrank+1
      iend=nunit*(myrank+1)
      sum=0.0e0
      IF(myrank==isize-1) THEN
        DO i=istart,iend
          psum=0.0e0
          DO j=1,nunitsize
            psum=psum+0.1e0
          ENDDO
          sum=sum+psum
        ENDDO
        IF(namari/=0) THEN
          asum=0.0e0
          DO i=1,namari
            asum=asum+0.1e0
          ENDDO
          sum=sum+asum
        ENDIF
      ELSE
        IF(nunitamari>myrank) THEN
          iend=iend+1
        ENDIF
        DO i=istart,iend
          psum=0.0e0
          DO j=1,nunitsize
            psum=psum+0.1e0
          ENDDO
          sum=sum+psum
        ENDDO
      ENDIF
      IF(myrank==0) THEN
        DO i=1,isize-1
          CALL MPI_RECV(buf,1,MPI_REAL,i,MPI_ANY_TAG,
     &                  MPI_COMM_WORLD,istatus,ierr)
          sum=sum+buf
        ENDDO
      ELSE
        CALL MPI_SEND(sum,1,MPI_REAL,0,1,
     &                MPI_COMM_WORLD,ierr)
      ENDIF
      IF(myrank==0) THEN
        WRITE(6,200) sum
  200   FORMAT('The answer is',e15.8,'.')
      ENDIF
      CALL MPI_FINALIZE(ierr)
      END
Figure 4. Parallelized Program of Summation with Partial Sum. (Program2)

用いた並列計算機は、産業技術総合研究所先端計算機センター(TACC)に設置されている並列計算機Alta Technology AltaCluster 、Hitachi　SR8000 、IBM RS/6000 SPおよびSGI Origin2000を用いた。計算機の詳細をTable 1に示した。なお、並列計算は1、2、4、8CPUを用いて行い、プログラムの実行時間は、MPIの初期化、ノード数およびランクの取得直後から終了処理の直前までの時間を5回計測し、最短時間を採用した。また、コンパイル時の最適化等のオプションは使用しなかった。

Table 1. Details of Computer Architecture
AltaClusterHitach SR8000IBM RS/6000 SPSGI Origin2000
CPUPentium III (katmai)SR8000POWER3R12000
CPU数851225664
クロック周波数(MHz)450250200400
メモリサイズ256MB/CPU8GB/8CPUs960MB/2CPUs2GB/CPU
１次CachSize(データ＋命令）16KB+16KB128KB+64KB64KB+32KB32KB+32KB
ネットワーク種ether (10Base-T)多次元クロスバーSPスイッチNUMA SMP
Fortranコンパイラ名PGI Fortran最適化F90コンパイラXL Fortrran for AIXMIPSpro F90
Fortranコンパイラバージョン3.101-037.17.3

3 結果と考察

Table 2にプログラムの実行時間および速度向上率を示した。Figure 5はTable 2の結果をグラフ化したものである。

Table 2. Execution Time and Speed-up Ratio* in Number of Processors
Number of Processors
1248
Execution timeSpeed-up ratioExecution timeSpeed-up ratioExecution timeSpeed-up ratioExecution timeSpeed-up ratio
sec---sec---sec---sec---
Alta Cluster
Program 114.4118.801.644.453.242.216.52
Program 217.781 8.882.004.444.002.218.05
SR-8000
Program 182.57141.112.0120.564.0210.308.02
Program 289.23144.512.0022.314.0011.158.00
RS/6000 SP
Program 134.72117.342.008.674.004.348.00
Program 234.73117.352.008.684.004.348.00
Origin2000
Program 125.63 112.762.016.374.023.198.03
Program 225.63112.832.006.414.003.217.98
* Speed-up ratio = execution time with 1 processor / execution time with N processors.


Figure 5. Execution Time and Speed-up Ratio in Number of Processors

実行速度は期待されるようにCPUに応じてほぼ理想的な高速化が図られた。これは、プログラム1、2共に計算量に対して通信量が非常に小さく、比較的低速なTCP/IPを利用するAltaClusterにおいてもMPIの通信オーバヘッドが問題にならないことおよび負荷分散がきちんとできているためである。
Table 3に２つのプログラムによる和の値および誤差率を示した。

Table 3. Calculated Value and Relative Error * in Number of Processors
Number of Processors
1248
Calculated valueRel. error*Calculated valueRel. error*Calculated valueRel. error*Calculated valueRel. error*
---%---%---%---%
All computers
Program 12.09715E+06-97.94.1943E+06-95.88.38861E+06-91.61.67772E+07-83.2
Program 29.99998E+07-0.000259.99995E+07-0.000509.99990E+07-0.000989.99980E+07-0.00198
*Rel. error = ( calculated value - exact value) / exact value00.0 ： (exact value=1.000E+8).

計算精度では、プログラム1および2においてそれぞれ4つの並列計算機で全く同じ結果が得られた。これは、現在製造されているほぼすべてのCPUは2進浮動小数点数の標準規格IEEE754に従っており、今回用いた計算機では浮動小数点数の取り扱いに違いがないために、結果が同一になった。また、コンパイラによる最適化を行わなかったことも、同一の結果を与えることに寄与している。
プログラム1は、非常に単純なプログラムであるにもかかわらず、良く知られているように全く正しい解を与えない。これは、“情報落ち”という誤差の累積によるものであることが広く知られている。Table 3をみると予想どおり並列処理を行うことによって情報落ちが少なくなるため、CPU数を増加させ処理を分散させることにより精度が向上しているが、８CPU程度では十分な計算誤差の回復には至っておらず、でたらめな答えとなっている。Figure 6 に示す回帰直線（y=-2.1x+100.0）から外装すると、50CPU程度を用いた並列処理を行うと正しい計算結果を与えるようになることが期待される。詳細な議論は行わないが、SGI Origin2000 60CPUまでのprogram 3の実行結果をFigure 7に示した。50CPU程度までは、直線的に精度が向上している。Figure 6の回帰直線からも示唆されるように50CPUよりも少ないところで正しい値(1.0E+08)となり、さらに用いるCPU数が60程度に大きくなるとまたエラーが増大するが、50CPU程度までにみられるような直線的増加ではない。


Figure 6. Relative Error (%) by Programs 1 and 2


Figure 7. Calculated Values of Program 1 on SGI Origin 2000

プログラム１とプログラム２は原理的に同一の解を与えるはずであるが、数値的には全く異なる解を与える。この例は、プログラムの書き方によって計算結果が異なってくることの良い例となっている。
プログラム2は、部分和を取ることにより情報落ちを回避するプログラムとなっているためほぼ正しい解が得られている。しかし、Table 3およびFigure 6から判るように、並列処理を行うことで、非常にわずかであるが見かけ上精度が悪くなっており、並列分散処理を行うと精度が向上するという一般的な傾向とはずれている。これは、以下のような考察により理解できる。まず、１つのCPUにおいて104回の部分和が数回計算されるが、この時丸め誤差を含んだ数を加えたことによる誤差が蓄積され、例えば0.1の104回の部分和は、正確な値（1000.0）にはならず、999.902…と誤差を含んだ値をとる。部分和の和をとっていくごく初期の段階では、部分和の和と各部分和の大きさがそれほど違わないため、丸めに起因する部分和の段階での誤差が次第に蓄積していく。ところが、部分和の総和が、各部分和の絶対値に比べて非常に大きくなっていくと、部分和段階での演算誤差が丸めにより切り上げられ、偶然真値（1000.0）として加えられる。そのため部分和の和をとるごく初期の段階を除けば、部分和の和に誤差が蓄積しない。逆に並列処理を行うと部分和の和のごく初期にたまった誤差が和になって表れ、CPU数が増すことにより誤差が比例的に増す。したがって、プログラム１および2の並列実行時にみられた現象は、MPIによる並列処理に起因する誤差でなく、部分和の取り方の違いによるよる現象であると考えられる。これは、1CPUを用いてMPIでの並列処理による和の取り方をシミュレートすることにより全く同様の結果を得たことからも確かめられた。

4 まとめ

2進循環小数の0.1を109回加えるプログラムを例にとり、MPIによる並列プログラムの実行速度の高速化および計算精度の挙動を4種類の並列計算機を用いて測定した。２つのプログラムにおいて、通信量は演算量に比べて非常に小さい。そのためいずれの並列計算機でも、ほぼCPU数に応じた高速化が図れた。演算精度はそれぞれのプログラムとも4つの並列計算機で全く同一であった。0.1を逐次加えていく場合（Program 1）では、並列処理により著しい精度向上が観測された。計算誤差の蓄積が非常に少ない104回の部分和を105回加える場合(Program 2)ではCPUが増加することにより、非常にわずかな精度劣化が見られた。これは部分和の取り方による誤差の挙動変化であり、この現象は偶然生じた現象であると考えられる。
並列処理をおこなうことで、これまでに扱うことが困難であった大規模な数値計算を短時間に行うことが可能となってきた。しかしながら演算回数が増えるにつれ、丸め誤差や情報落ちなど、数値誤差が入り込む部分も多くなるため、小規模な数値計算の際にはあまり問題にならなかったような場合においても、誤差が小さくなるような計算手順の工夫が必要である。

この研究は静岡県職員資質向上研修として行ったものである。計算は独立行政法人産業技術総合研究所先端情報計算センター(TACC)で行われた。

参考文献

[ 1] http://www.top500.org
[ 2] 湯淺,安村,中田編, はじめての並列プログラミング, bit別冊, 共立出版 (1998).

Return