図 2 グリッド乞の次の日製十 1 の量を決めるために要な情報 図 1 をグリッドに分割 2 〇 - Q 、 ちょっと注意しないといけないのは、両端のグリッドで 差分法による解法 す。たとえば、左端のグリッドを引算するときには、さら にその左側のグリッドの情報が必要になります。この情報 以下では、差分法を用いて拡散方程式を数値的に解きま を与えるには、両端の条件を決めなければなりません。 す。 差分法では、空間をグリッドに分割し、そのグリッド点 問題設定 での変化量を数イ缶 1 算します ( 図 1 ) 。 差分のとり方にはいろいろな方法がありますが、 はシンカレに、グリッドは等幅で、空間は 2 次精度としま す。また、時間は 1 次精度とし、陽解法で解くことにしま す。△ェ、△″、△ 2 をそれぞれムリ、 2 方向のグリッドのサ イズ、△を時間ステップの幅とすると、差うには次のよう になります。 こで、問題を定しておきましよう。まず、簡単にす るため、次元を落として 1 次元の拡散問題とします。また、 初期条件として、サイン波、 佖は , t = 0 ) = s ⅲ ( 2 冗工 ) を考えます。また、境界条件は周期境界条件とします。 れによって、左端のグリッドのさらに左側、右端のグリッ ドのさらに右側にもデータが与えられ、両端のグリッドの 1 算が可能になります。 この初期値問題の場合、以下のように変数分離法で簡単 に厳密解カ球められます。 = exp(—4T2t) s ⅲ ( 2 冗ェ ) 奴ェ , の に付いている添字は、下付きの乞 , んが空間のグリ 図 3 は、この進化を示したものです。初期に与えられた ッドの位置を、上付きの / カ塒間を表します。すなわち、 サイン波は時間の経過とともに減衰していき、最後 ( 厳密 朝 , んは空間のグリッド位置を、 / 十 1 は / の次の時間ステ には無限時間後 ) には振幅ゼロになります。 ップを表しています。 プログラム 次の時間 / 十 1 の値は、現在の時間 / の全グリッドの情 というわけで、実際にプログラムを作成してみましよう。 報を使って言算します。その引算を次々とおこなって、進 末尾のリスト 1 は、私が C 言言吾で書いたサンプル・プログラ 化を追っていくことになります。上記の差分式から明らか ムです。与えられた差分式はそれほと灘しくないので、そ ですが、グリッド 0 物 , ん ) の次の時間 / 十 1 の量を決める の他の言語で同種のプログラムを書くのは比較的容易だと ために必要な情報は、時間 / でのそのグリッドと、それに 思います。 隣接するグリッドです ( 3 次元の場合は、隣接するグリッ ドは全部で 6 つになります ) 。 1 次元であれば、図 2 のよ このサンプル・プログラムについて、簡単に解説してお うなイメージになります。 きましよう。 ~ 十 1 をツ , ん 20 を .7 , ん十・佖を一 1 , 去ん ( △の 2 2 佖 . ノん・十 0 朝一 1 , ん ( △の 2 2 佖乞ムん - ト a を朝 , ん一 1 ( △の 2 0 を十 1 , 去ん 佖ツ朝 , ん △ 0 を , j 十 1 , ん 十 十 0 を朝 , ん十 1 ( 1 ) 72 UNIX MAGAZINE 2006 . 3

表 1 グリッド数を変えた場合の計のヒ ( 104 ステッフ。の所要時間 ) ICPU 2CPU 4CPU グリッド数 4 , 096 16 , 384 65 , 536 262 , 144 1 , 048 , 576 4 は 94 , 304 time(s) 0.56 2.22 8.88 35.99 143.6 574.8 %cpu 96.4 99.5 99.4 99.0 99.7 99.5 time 5.48 7.76 9.89 25.4 81.25 305.6 ロ 0.6 %cpu 5.1 17.9 28.1 74.3 90.8 95.2 rank0 rankl rank2 rank3 ICPU 0.8 ロ tirne 6.91 7.04 10.20 15.95 48.21 167.3 図 7 single & MPI の上肆交 、ロ . ロ 0.4 0.2 ー 0.2 ー 0.4 0 ロ 戸 0.2 0.4 1 に対応しているのは偶数個の場合だけです。 76 になると ( なんと ) 48 分、 4CPU では 56 分もかかってい てしまいます。 ICPU では 73 秒で終った言 t 算が、 2CPU このケースでは、 MPI 化によって計算時間が遅くなっ ったく同じです。 らも明らかですが、答は 1 台の言 t 算機で引した場合とま 果を、、ロ " で、 MPI の結果を線で表しています。この図か は完全に一致します。そのため、さきほどの ICPU の結 並列化によって言算川印茅カ畯わることはないので、結果 列化して得られた結果を以前の結果と比較します ( 図 7 ) 。 ます、グリッド数を以前と同じ 1024 にし、 4CPU で並 のように実行します。 % mpirun -np 4 onedim-mpi プログラム onedim-mpi を、 % mpiCC ー 0 onedim—mpi onedim—mpi . C す ) で保存し、次のようにしてコンパイルします。 これを適当なファイル名 ( 今回は onedirn-mpi. c としま %cpu 16.4 18.8 58.9 78.2 87.4 ます。これは、計算よりも通信に多くの時間カ饋やされる からです。グリッド数が 1024 、すなわち 4 並列の場合、 ICPU あたりのグリッド数が 256 ではさすがに少なすぎ るようです。予想していたこととはいえ、ちょっと悲しい 結果ですね。 グリッド数を変えた場合 今度は、グリッド数を多くすれは並列化の効果が得られ るのかをみてみましよう。 このテストプログラムで使っている時間刻みの決め方 は、グリッド数の 2 乗の逆数に比例するようになっていま す。したがって、グリッド数を多くすると、時間刻みも小 さくなります。ある時間まで計算させるようにすると、グ リッド数が多いところではステップ数も増えて計算に時間 がかかってしまいます。そこで、今回は 10 , 000 ステップ の言 t 算に要した時間を測定してみました ( 物理的には、あ まり前未のない言 t 算になってしまいますが・ time コマンドを付けて mpirun を実行した場合の引算 時間を表 1 に示します。 この表のグリッド数を CPU 数で割った値が、各 CPU が受け持っグリッド数ということになります。たとえば、 グリッド数が 4 , 096 だとすると、 2CPU であれば ICPU あたりのグリッド数は 2 , 048 、 4CPU であれば 1 024 と なります。 この結果から分かるのは、今回の場合、グリッド数カ 十万以上あれば並列化の効果カ硼待できそうだということ です。グリッド数が 100 万以上になると、引時間はほば CPU 数に反比例するようになります。 グリッドが少ないと、ほとんどは通信を待っているら しく、 CPU の占有率がかなり低下します。グリッド数が 4 , 096 で 2CPU の場合をみると、 CPU はたった 5.1 % しか使われていません。これでは、いくら並列化しても、 UNIX MAGAZINE 2006.3

図 5 グリッドの分割 GRID CPUO 図 6 グリッド情幸医のイメー CPUO [ わ CPU3] CPUI CPU2 CPU3 CPU3 ン CPUI CPU2 MPl_Sendrecv のペア MPI Sendrecv のペア というラベルが付いている曲線 ) 。一見して明らかなよう に、時間刻みを長くしたものは厳密解から大きくずれてし まいます。 グリッド数についても同じようなことがあり、すくなく とも、調べたい現象に対して十分に空間的に分解するグリ ッド数が必要になります。 適切な時間刻みやグリッド数の決定は、いわは職人芸的 MPI による並列化 なものカ球められるところといえるでしよう。 まずは、言 fr 算の分割です。今回の例では、図 5 のように 並列化を考える 合や、グリッド数と効果の関連性を紹介しましよう。 たグリッド数が少ないからです。グリッド数を多くした場 い、もしくはまったくありません。これは、例として使っ 場合は、たとえ並列化しても、その効果があまり得られな ち、計算量はそれほど多くありませんでした。このような さきほと挙げた例は、計算時間がきわめて短い、すなわ Sendrecv を使った実装を紹介します。 ます。以下では、もっとも簡単な方法として、関数 MPI- 続いて、 MPI を用いてプログラムの並列化をおこない 74 [toCP 0 ] グリッドを CPU ぶんに分け、各 CPU はそれそれ分割さ れた 1 つの領域の計算を受け持つようにします。これで、 言 1 算量が CPU の数ぶんに分割されます。 さきほど述べたように、グリッドの 1 点を計算するに は、その両隣のグリッドのデータカ坏可欠です。分割され た境界のグリッド点を言 t 算するには、その隣の領域にある グリッドの情報が必要になりますから、隣の領域を担当す る CPU からそれを送ってもらわなければなりません。 れは、いわゆる、、袖 ( そで ) 輯医 " と呼ばれる方法で、図 6 は そのイメージを表したものです。この図は 4CPU 構成の 場合で、・は各 CPU カ授け持っデータです。また、デー タの送り先を矢印で、送った先のデータ ( コピー ) を〇で表 しています。 このように、今回の例の場合は境界のデータのみを通信 すればよいため、通信量はかなり少なく抑えることができ ます。とはいえ、グリッド数が少ないと通信に大きな影響 を受けるので、せつかく並列化しても通信コストによって 計算が遅くなることがあります。 並列化プログラム それでは、実際に MPI の関数を使って並列化してみま しよう。 UN 工 X MAGAZINE 2006.3

連載 / 天文学と UNIX 図 4 プログラムの結果と厳懈の上肆交 0.4 図 3 1 次元の拡散 1- 8 ( 0 4- つ」 0 0•J 4 ・ ( 0 8 ィー 0 analytic numerical rough -1 0•J 一 4 8 0.2 0 0 ー 0.2 ー 0.4 1 0.8 0.6 0.4 0.2 0 基本的に、 main 関数のなかでグリッド点のデータが入 言算の正当性が検証できます。あくまでも見た目にすぎま る配列 a を定義し、初期条件を作成する関数や数イ t 算を せんが、計算結果と厳密解はそこそこ一致しているので、計 おこなう関数を呼び出すようにしています。用意する要素 算はそう大きくは間違っていないようです。厳密な結果を 数は、実際に必要なグリッド数 (NGRIDREAL と de- 求めるのであれば、各グリッドで厳密解と数値解の差をと fine 文で定義 ) に、その両隣のぶんを足した数 (NGRID- るなどして、定量的に検証する必要があります。 REAL 十 2 と NGRID で定義 ) です。 このようなシミュレーションは、物理やプログラムに関 グリッドのデータは、配列インデックス I—NGRID- するちょっとした知識があれば、簡単におこなえます。と REAL のなかに入ります。 0 番目と NGRIDREAL 十 1 はいっても、あらかじめきちんと考えないと、観察したい 番目には、周期境界条件から、 NGRIDREAL 番目と 0 現象が再現されなかったり、計算に時間がかかりすぎて実 番目のデータを入れます。こういった、、のりしろ " のよう 際上は引算不能ということにもなるので、注意が必要です。 な領域を用意し、すべてのグリッドの計算を 1 つの引式 このあたりをよく考慮しなければならない点が、シミュレ でおこなうことにします ( そうしない場合は、両端とそれ ーションの難しいところではありますが・ 以外とで計算式を分ける必要があります ) 。 上記の問題にしても、本当はグリッド数や時間の刻み幅 このプログラムを適当なファイル名 ( 今回の例では、 について十分に検詞する必要があります。 onedim. c とします ) で保存し、 たとえば、時間の刻みを考えてみましよう。観察したい % gcc onedim. C ー 0 onedim ¯lm 現象の時間が決まっているのなら、刻みを長くしたほうが として実行ファイル onedim を作ります。 計算に必要なステッフ。数が少なくなります。その結果、計 コンパイルができたら、プログラムを実行してみましょ 算の所要時間も短くなります。ところが、時間の刻みをあ う。図 4 は、時間 t ー 0.02 まで計算した結果を表した まりに長くすると、現象が正しく再現されなくなる可能性 ものです。この引算では、グリッド数を 1024 としていま があります。プログラムの中では、この時間刻の大きさは、 す。実際にかかった計算時間は、 amanogawa0 で 73 秒 変数 KOEFF-DT でコントロールします。 KOEFF-DT カ賠になると、時間刻も倍、すなわち、必要なステップは半 でした。 分になります。先の結果は、 KOEFF-DT= 0.00390625 図 4 の実線は、式 1 から計算した厳密解です。これは、 としていますが、図 4 には、 KOEFF-DT= 0.125 という 式 1 に t = 0.02 を代入して得たものです。今回の場合の 時間刻みをかなり粗くした結果も示してあります (rough ように厳密解があるときは、厳密解と比較することで数値 1 0.8 0.6 0.4 0.2 0 73 UN 工 X MAGAZINE 2006 . 3

/ / 時間を進める t + = dt ; + + step , if( ! (step%100) ) { if (myrank== の { printf ("step=%d, return(step) ; int main (int argc , char **argv) { MPI_Comm—rank (MPI_COMM—WORLD , MPI—Comm-size (MPI-COMM-WORLD , MPI_Init (&argc , &argv) ; int nlocalgrid, tmpb; static double a CMAXARRAY] ; / / 以下で隣の CPU を決める while(t く tend){ / / 境界のテータを MPI を使ってやりとりする 十 = MPI_Sendrecv(&(a[nIoca1grid] ) , 1 , 1 , MPI_COMM_WORLD , &stat) ; MPI-Sendrecv(&(a[1] ) , 1 , MPI_DOUBLE , MPI_DOUBLE , MPI_DOUBLE , MPI_DOUBLE , right , left , left , right , 0 , 1 , &(aCn10ca1grid + 1] ) , 1 , MPI_COMM_WORLD , &stat) ; a[i] + = dtdx2inv*( a[i + 1] for(i=l ; 土く =nlocalgrid; i + + ) { / / 積分 step, &myr ank ) ; &numprocs) ; numprocs ; numprocs ; / / おまじない / / おまじない / / おまじない left =myrank ー 1 ; right=myrank + 1 ; if (left ー 1 ) left if (right== numprocs ) right / / 各 CPU でのグリッド数を計算 nlocalgrid = NGRIDREALTOTAL tmpb NGRIDREALTOTAL / numprocs ; % numprocs ; fprintf (stderr , "myrank=%d , myrank, nlocalgrid, init—sin—mpi(nlocalgrid, a) ; integrate—mpi(nlocalgrid, a) ; output—mpi (nlocalgrid, a) ; / / 積分をおこなう / / 初期条件を作る tmpb, left, right) ; N/nproc=%d, N%%nproc=%d , left=%d, / / 結果を出力 ( 各 CPU が独立に出力 ) / / おまじない MPI—FinaIize ( ) ; return(O) ; 80 right=%d\n" UNIX MAGAZINE 2006.3

SC 翡 ネットワーク・ テクノロジー 好評発売中 ! 前回紹介した関数 MPI-Send や MPI-Recv を使って プログラムを書くこともできます。これらはペアで使う必 要があります。たとえば、 CPUO と CPUI とのあいだ で通信する場合、 CPUO で MPI-Send を呼び出すとする と、 CPUI ではさきに MPI-Recv を呼ばなければなりま せん。さもないと通信が成立せず、デッドロック状態に陥 ってしまうからです。 CPU (MPI では rank) によって 振舞いを変えるには、 C 言語であれは次のように書く必要 があります。 if ( myrank & 0X1 ) { / * myrank が奇数なら実行 * / MPI_Send() ; MP 工 _Recv() ; else{ / * myr k が偶数なら実行 * / MPI_Recv() ; MPI_Send() ; ネットワーク - アクノロソー インターーネットを支える技術・様ロ陽一 8 ・・播ロ陽ー著 ・ A5 判、 448 ペーシ ・ ISBN 4-7561-4092-0 ・ 3 , 675 円 ( 税込み ) ネットワーク技術者必 ATM 、、′ T などーこついて、その 技術の“心”が分かる本 ネットワーク技術者必携 ! この例では、最初に奇数の rank が送信し、偶数の rank が受信するとします。そして、その後に逆のことをおこな 技術には、かならず核心となる部分があります。これ を理解しているかが、その技術を使いこなせるかどう います。 かの分かれ目といえるでしよう。 このようにプログラムを書いてもいいのですが、通信と 本書では、現在のインターネットを支えるさまざまな 送信の面倒を同時にみてくれる関数 MPI-Sendrecv を使 技術について、開発の背景から、原理、応用まで、 うと、より簡潔に謎できます ( この関数の詳細は、前回に 数多くの図版を用いて詳細に解説します。 挙げた参考文献を参照してください ) 。 この関数を用いて MPI 化したプログラムが、末尾のリ 目次から スト 2 です。 1 章 Ethernet 2 章リピータ、プリッジ、スイッチ このプログラムについても、簡単に解説しておきましょ 3 章ス / ヾニングッリー う。まず、 ma ⅲ関数のなかで MPI の、、おまじない " をし 4 章 VLAN ます。自分の rank を取得したあと、両隣の領域の rank 5 章経路制御の基礎 6 章経路検索アルゴリズム ( 変数名 left 、 right) を定します。基本的には myrank 7 章 ATM に対して土 1 したものですが、周期境界条件を考慮してい 8 章 UNI シグナリング ます。 9 章 LANE と MPOA 初期条件も、自分が受け持っ領域のデータのみを作成し 1 0 章アドレス変換ーー NAT と旧 masquerade 1 1 章 ADSL ます。関数 integrate-mpi のなかでは、グリッド点の計 参考文献 / 索引 / 英和用語対照表 算を始める前に MPI-Sendrecv を 2 回呼び出し、計算 に必要なデータを用意しています。それぞれの関数 MPI- 株式会社アスキー Sendrecv を呼び出した際に奐されるデータについては、 図 6 を参照してください。 このプログラムでは、実行時に一 np オプションの引数で CPU の数を指定できるようにしています。ただし、実際 〒 ] 02 ー 8584 東京都千代田区九段北ト 1 3-5 日本地所第一ビル 電話 (03) 6888 ー 5500 ( 営業局 ) 75 UNIX MAGAZINE 2006.3