From VPP Fortran to HPF

「VPP FortranからHPFへ」

荻野竜樹　名古屋大学太陽地球環境研究所
Tatsuki Ogino Solar-Terrestrial Environment Laboratory, Nagoya University

１．はじめに

1995年12月にUCLAを訪問した時、共同研究者であったDr. J. Raeder達が「日本は間違った 方向を選んだ。日本はベクトル並列機を選んだが、米国の多勢は超並列計算機へ向かっており、 自分たちはHPF(High Performance Fortran)を用いた同じプログラムで米国の主要な超並列計算 機を利用できて、しかも満足できる性能を得ている。HPFへの書き換えはかなり時間を要する ので、これで近い将来日本の研究者に大きく差をつけることができる」と言って私を驚かせま した。当時はcpu数の競争をしたCM-5などの超並列計算機が一世を風靡した後で、超並列機 のcpu数も100〜2000個で実質的な効率がより重要であると再認識されてきた頃でした。 CRAY-T3D、PARAGON、IBM-SPなどの新世代の超並列計算機が現れて、cpu数に比例した性 能がでるコンパイラーの開発と並列機の性能の指標としてスケーラビィリティーの重要性が指摘 されていました。

私もHPFに大変興味を持ち、HPFが動く超並列機の導入を計画するとともに、数社と協力して 私達の太陽風と地球磁気圏相互作用を解く３次元グローバル電磁流体力学（MHD）コードを Fortran 90とHPFに書き換えてテストする好機会を得ました。空間３次元を全て１次元配列と して扱いプログラムサイズを最小化した１次元化MHDコードは、CRAY Y-MPやVP100など のベクトル機で効率よく動きました。しかし、超並列機で並列化することは全て不可能という ことになり、空間３次元は３次元配列に戻した３次元MHDコードを用いて、改めてHPFへの 書き換えと計算効率を調べました。結果は、いずれも簡単なテストコードでは性能が出ました が、２次元と３次元MHDコードでは期待したほどの計算速度は得られませんでした。プログ ラムの書き換えを何度も試みなかったことはありましたが、特にcpu数が増えた場合の性能は 悲惨なもの（数％以下と記憶しています）でした。しかし、米国の研究者がHPFでかなりの成 果を上げていることは計算結果から見ても十分に信用できるものでした。

名古屋大学大型計算機センター（以後、大計センターと書く）のスーパーコンピュータでHPF が使えるようになったら教えてほしいと頼んでいましたが、昨年度あるメーカーの人がHPFの 開発が日本でスーパーコンピュータメーカー３社と科学技術庁の協力のもとに進められている と熱気を持った講演をされたのを聞く機会がありました。大変興味を持ったのと同時に、反面 本当にどこまでできるのであろうかと疑ってもいました。富士通のVPP Fortranで、既に高速計 算を実現しているのだからせめてその半分以上の性能を出せなければ結局は利用しないことに なります。しかし、VPP Fortranで性能がでていても、NECや日立のスーパーコンピュータで は高性能がでないままでは国内外の研究者と共同研究などをする上で大変困るというもう一つ の深刻な問題が現実にあります。

そうした状況下、この６月に富士通の講師による「HPF講習会」が大計センターで開かれまし た。この講習会を聞いて、VPP FortranからHPFへの書き換えが悪くても何とかなるという印 象を持ちました。その漠然とした理由は、「通信制御ができそうだ」と「配列の分割方向を高速 に切り換えれば、並列化（ベクトル化も含めて）の困難な部分は原理的に解決できる」と思っ たからでした。すぐに、富士通提供のサンプルプログラムのテスト、私達が持っている３次元 波動方程式などのテストプログラム及び３次元MHDコードのVPP FortranからHPF/JA（HPF の日本仕様拡張版）への書き換えを始め、約３週間で大体の書き換えを終えました。本来の目 標であった３次元グローバルMHDコードの計算結果確認には、その後更に２週間を要しまし たが、最終的にはHPF/JAの３次元MHDコードを用いて同じ実用的実行形式でVPP Fortranと 同等の計算速度を得ることができるようになりました。

２．いままでのコンパイラー利用の経験

これまでいくつかのFortranコンパイラーを利用してきましたが、周りの人がほとんど誰も使っ ていない早期に新しいコンパイラーを使い始めたことが何回もあったので、プログラムを作成 する時にいろいろなことを調べて考えると共に、様々なテストをしてきました。私がコンパイ ラーとプログラムに関して後で示すような考え方にどうしてなったのかを理解して頂くために、 印象が強く残っているコンパイラーに関して、その利用の中で考えてきたことと学んできたこ とをコンパイラーの利用の経験として以下にまとめてみます。

２−１　Fortran 77 コンパイラーの利用
Fortran 66によるプログラムなどで非常に不便を感じていたものに、サブルーチンの引数の引き 渡しがありました。サブルーチンの引数の大きさに制限があるのが普通であり、またサブルー チンに引数を増してサブルーチンを拡張するような場合、その都度新しいサブルーチンを作る ことになっていました。この不便さを何とかしたいと思っていた時、日本大学の川上一郎先生 が「外国ではサブルーチンの配列を dimension p(1) などと書くのがはやっている。これは間違 いではなく、 Fortran 77 では許されている便利な記述法である」と紹介されました。これは、 便利で面白いと思って、それからすべてのプログラムを１次元配列表現で書き改めました。そ の要点は、次のようにまとめることができます。
　（１）すべての配列（多次元でも）に１次元配列表現を用いる。
　（２）配列の大きさはすべてメインプログラムのパラメータ文で与える。
　（３）サブルーチンでは配列の大きさは与えず、すべて p(1) などと１次元１要素宣言で
与える。また、必要な変数の形（整数と実数）の１次元配列を用意する。

その利点として、
　（１）変数の受け渡しと入出力が簡単になる（書式無し入出力文の利用）。
　（２）変数（配列）の大きさを変える時、メインプログラムのパラメータ文の変更だけで
　 すむ。サブルーチンは全く考慮する必要がない。
　（３）サブルーチンで新たな引き数を増加して書き換える場合も、１次元引数の後にいく
　　　　らでも追加できるので、サブルーチンの形を保ったままで自由に書き換えることが
　 できる。
　（４）計算速度は１次元配列化によって確実に速くなる。

このように Fortran 77 の利用と１次元配列の利用はプログラムの高効率化と利便さに大き なメリットがあったと言えます。

２−２　CRAY CFT コンパイラーの利用
CRAY-1は1MWと主メモリが小さかったので、そのメモリ制限のもとでベクトル化の効率を 上げるために、すべての主要なdoループをベクトル化すべきであることが課題でした。それと Fortran 77で用いてきた１次元配列の利用に固執しました。３次元配列利用の方がベクトル化 が容易であることはわかっていましたが、１次元配列の方がやはり計算速度が少しでも速い ので、１次元配列を保ったままで、作業変数の最小化とすべての主要なdoループのベクトル 化を実現しました。地球磁気圏の３次元グローバルMHDコードにおいて、プログラム全体の サイズは必要不可欠な変数の1.3倍以下になりました。こうして、1MWの主メモリを持つ CRAY-1で、境界を含めて62*32*32 (必要不可欠な変数の容量は0.51MW)及び50*50*26(必要 不可欠な変数の容量は0.52MW)の３次元格子点を用いたMHDシミュレーションを行うことが できました。

作業変数を最小化する計算方法としては、z方向のある面でMHD方程式を解く場合作業領域は その上下面を加えた３面しか持たず、そのz方向の面を移動しながら計算を続ける方法を採用 しました。１次元配列を用いて、ベクトル化を行ったもう１つのポイントは、パラメータ文 による定数指定（データ領域のセパレーションを明示）と、変数配列の順序をうまく並べる こと（３章の(B1)の例を参照）で１次元配列のオーバーラップの計算が生じないことを、予 めコンパイラーに知らせることができることを利用しました。これらの方法を実現する上で プログラムの構造をきちんと決めることが重要であると学びました。

２−３　アデトランコンパイラーの利用
VPP Fortran にいく前に、京都大学の野木達夫先生が考案され、松下電器産業株式会社が試作し て、私達も受託研究で十分に利用する機会を得たADENARTのアデトランコンパイラーについ て簡単に紹介します。ADENARTは、ベクトル化機能は全く無しの超並列計算機の一種でした。その 並列計算の本質は、３次元などの場合各CPUに分割する空間方向を高速に切り換えてデータの 並び替え（通信）を行い、各CPUで高速計算することを駆使する方法でした。ADENARTを 有効利用するためには、空間の次元の方向とデータ分割の方向を一致させる必要があったため に１次元配列を完全に諦めて、３次元空間を３次元配列に戻しました。500MBとメモリはか なり大きかったのですが、配列分割方向の異なる実質的同じ変数を多重に定義する必要性があ ったために、何か工夫しないとプログラムサイズがあっという間に大きくなる傾向がありまし た。作業変数はできるだけ少なく取って、同時に並列化の効率を上げるにはプログラムの構造 を決めることは大変重要であると思いました。結果として、プログラムサイズは必要不可欠な MHD変数の４倍以内に収まり、計算速度も256プロセッサ利用で１個のプロセッサの利用の 220倍に達しました。即ち、256プロセッサのADENARTで85.9％の並列化効率が得られたこ とになります。

２−４　VPP Fortran コンパイラーの利用
先ず、予想通り私達はVPシリーズまで利用してきた１次元配列利用が不可能になり、すべ てのプログラムにおいて３次元配列に変更するための全面的書き換えが必要になりました。こ れはかなりショックで、１次元配列のままで何とかならないかとかなり考えてみましたが、結 局どうしようもないという結論になりました。言い訳しますと、これは CRAY Y-MP C90を用 いてのベクトル並列化（自動ベクトル化機能を利用）では全くプログラムを書き換えずに8cpu で十分な性能（8 processorsで4.88Gflops）が出ていました。こうして空間３次元を３次元配列 に戻してベクトル並列化を始めました。１次元配列をすべて３次元配列に戻すのは、非常に手 間のかかる作業でした。その後のベクトル並列化の作業は順調とはいきませんでしたが、大き な困難もありませんでした。領域分割による並列化は CRAY Y-MP の時に既に経験していた ことで、問題が生じた点と言えば、配列の分割領域の扱い（袖の宣言とdo変数の範囲の指定の 両立性）及び境界条件での一括転送の利用とベクトル化の保持でした。VPP Fortran ではグロー バル変数とローカル変数をequivalence文によってうまく結びつけることによってメモリの有効 利用を図って分散メモリにおける並列計算を実現しているわけですが、正にそのequivalence文 の利用が Fortran 90などの新しい並列計算用コンパイラーで使用制限しようとしているもので あることも、非常に気になる点でした。しかし、VPP Fortran の利用でMHDコードのベクトル 並列化はすべての主要な内側のdoループでベクトル化し、最も外側のdoループで並列化する ことに成功し、ベクトル並列化の効率は絶対速度で比較して70％以上の高性能を得ることがで きるようになりました。

３．並列計算について

１度に１命令を高速に実行するノイマン型のコンピュータを更に高速化する方法として、パイ プライン方式のベクトル型とcpuが沢山ある並列型のスーパーコンピュータが同時に考えられ ていたと聞いていました。ベクトル型で平均してノイマン型の50倍程度、並列型で１万倍程度 の高速計算が可能であろうと予想されていました。最初にベクトル型が選ばれたのは当時cpu が高価であったためであり、その後cpuの価格が下がったために並列機、そして超並列機が世 の中に出てきました。一方これは、スーパースカラーとも呼ばれるRISCプロセッサによるス カラー計算の高速化にも大きく依存しています。

プログラミングでも、最初のベクトル機、CRAY-1のFortranプログラムを書く時には、ベクト ル型であるか並列型のスーパーコンピュータであるかは区別しないで、将来どちらにも通用す るような基本的プログラムを作るにはどうしたらよいかを考えていましたし、CRAY-1のプロ グラミングの考え方のガイドブックにもその様に説明してあったのを憶えています。そのエッ センスはプログラムの中でdependency（依存性）があるかないかはどうして決まるのか。そし てそのdependencyがあるかないかを予めコンパイラーにどうやって認識させるかでした。 dependencyがないということは、例えばdoループの中で変数の順序を変えても計算結果が変わ らない保証があれば、成り立ちます。そのようなdoループを完全独立なdoループと呼ぶこと にします。次にdoループの中で、変数の変化に伴う逐次計算ではなく、変数の順序を守れば計 算結果が変わらない保証があれば、これも条件付きでdependencyがないことになり、これを順 序付き（条件付き）独立doループと呼ぶことにします。従って、doループを小さいdoループ に分解していくと、多くのdoループが完全独立か、順序付き独立doループの要素に分けられ ます。もちろん、関数の級数計算などのdoループの変数に伴う逐次計算が必要なdoループは最 後までdependencyが無くなりません。

順序付き独立doループと順序付き依存doループの違いを具体例で示します。

　（A1）順序付き依存doループ
　　　　do  10  i = 2, nx
　      f(i)=f(i-1)
　（A2）順序付き独立doループ
　　　　do  10  i = nx, 2, -1
 　     f(i)=f(i-1)

（A1）と（A2）は計算結果が異なり、（A1）は依存性を持つdoループの例となります。（A2） は変数の順序を守れば、まとめて計算しても結果が変わらないことが分かると思います。この 具体例を見られると多くの人が「なんだそんなことか」と思われるでしょうが、順序付き独立 かそうでないかの本質はこの簡単なdoループに含まれています。問題のポイントはコンパイラ ーが独立かどうかを認識できるかどうかです。dependencyがあるかないかを決めるもう一つの 要素はif文などの条件文です。しかし、条件文は通常何かの方法（マスク、リストベクトル、 gather/scatter等）で解決することができます。

完全独立なdoループは、変数をどんな順序で計算しても結果が変わらないのですから、ベクト ル化はもちろん並列化もできるのではないかと予想していました。doループを分解していき、 新たな計算補助配列を追加すると多くのdoループを完全独立なdoループに分解できます。 （A2）などもその例になります。ベクトル化はもちろん問題ありませんでしたが、並列化する 場合には共有メモリであるか分散メモリであるかに依存し、共有メモリでは並列化の問題はあ りませんが、分散メモリの場合は、計算に必要なすべての変数が同一cpuメモリ上にあるかど うかに依存します。従って、同一cpuメモリ上にない変数は、計算に先立ち、他のcpuメモリ から自分のcpuメモリに持ってこなければなりません。

順序付き独立doループの場合はどうでしょうか。これはいかにして順序付き独立doループで あるかを予めコンパイラーに知らせることができるかがポイントになります。３つのdoループ の例を示します。

（B1）do  10  i=1, ny
      f(i+nx)=f(i+nx)+p(i)
      f(i)=f(i)+q(i)
（B2）do  10  i=1, ny
      f(i)=f(i)+q(i)
      f(i+nx)=f(i+nx)+p(i)
（B3）do  10  i=1, ny
      f(i,2)=f(i,2)+p(i)
      f(i,1)=f(i,1)+q(i)

nx≧ny＞0の場合は（B1）〜（B3）は本質的に同じプログラムで、実質的には完全独立なdo ループとなりますが、コンパイラーがそれを自動的に認識してくれるでしょうか。そのために はnx≧ny＞0の条件をコンパイル時に知る必要があります。これはparameter文でnxとnyを定 数指示することによって可能となりますが、nxかnyのいずれかが変数であれば、コンパイラ ーはnx≧ny＞0の条件を全く認識できません。また、（B1）と（B2）にも微妙な違いが出てき て、（B1）ではコンパイラーが順序付き独立であることを認識できるが（B2）では認識できな いことがしばしばで、CRAYやVP等ではそうでした。これは考え方としては、（A1）と（A2） の違いと共通しています。即ち、（B2）では変数iが増加し、nx＞0の場合、増加方向のf(i+nx) が f(i)の後で計算されるために、順序付き独立であることを認識するのがより困難になると考 えられます。（B2）から（B1）のように式の順序を並べ換えるだけでコンパイラーが順序付き 独立であることを知ることが容易になります。これが、CRAY CFTコンパイラーで１次元配 列計算を大きくまとめてベクトル化するために用いた方法です。

次に（B1）と（B3）の比較ですが、これは配列を１次元配列にするか、２次元配列のままに しておくのかの違いです。CRAY Y-MPでは（B1）の１次元配列も全て自動並列化ができま したが、HPFやVPP Fortranによる並列化では（B1）の１次元配列は取り扱うことができず、 すべて（B3）のように元の多次元配列に戻して並列化を行いました。何故（B3）が並列化でき て、（B1）の１次元配列が並列化できないのでしょうか。これは配列のセパレーション又は ブロック化の方法と認識の仕方が関係していると思われます。分割データの場合、分割境界で のデータのオーバーラップ（VPP Fortranでの袖、HPF/JAでのshadow/reflect）をいかに実施し て計算するかもありますが、１次元配列はVPP FortranやHPFでは全く並列化が実施できな いという大きな問題があります。これが、３次元MHDプログラムの並列化でCRAY Y-MP、 VP-2600、SX-3、S820等のベクトル機で常用してきた１次元配列利用のプログラムを全て棚上げ して、空間多次元配列をそのまま使用したプログラムに書き換えなければならなかった最大の 理由です。

４．VPP FortranからHPFへ（全体的な印象）

数年前からHPFが使えるようになったらすぐに使いたいという希望を大計センターに伝えて あったのですが、この６月に富士通の講師によるHPF講習会があるから是非来ませんかと誘わ れました。欧米ではHPFがあまり成功したとは言えず、ユーザー離れが起こっているのを知っ ていたので、期待半分で出かけました。もちろん、他方では富士通はVPP Fortranの実績を基礎 にして、HPF/JAを開発したのですから、何とかなるのではという期待も持っていました。それ で、VPP Fortranの半分以上の性能が出るようならHPF/JAに本格的に書き換えようと考えまし た。前にも述べました通り、私達はいくつかのテストプログラムと地球磁気圏の３次元MHD コードでフルベクトル化とフル並列化を行ってきた経験を持っています。それでHPF講習会を 聞いて、VPP FortranからHPF/JAへの書き換えは悪くても何とかなると感じました。その要点 は次のようにまとめることができます。

（１） 領域分割（３次元のz方向）でデータ分割を行っているが、それはそのまま使える。ま
　　　た、分割データにおける境界でのデータのオーバーラップはVPP Fortranでは袖の利用
　　　で行っているが、これはHPF/JAのshadowとreflectの利用に置き換えることができる。
（２）並列化の手続き分割がindependent文によってVPP Fortranと同等に実行できる。
（３）データの一括転送をasynchronous文によって実行できる。
（４）データの一括転送を行うasynchronous文は、領域分割方向を切り換えたデータ分割配列
　　　間にも適用できる。
（５）従来、HPFで問題であった不要な通信をindependent、new(　　)節とon home、
      local(　　)節によって指示することができる。
（６）最大値や最小値などがreduction節によって並列計算できる。

（１）、（２）、（３）からVPP Fortranで書かれたプログラムはHPF/JAへの書き換えがほとんど そのままできそうに感じられました。また、HPFで従来問題であった、不要な通信を（５）で 完全に止めることができれば、並列化の効率をかなり上げることができます。これも並列化の 手続き分割によって計算結果が変わらない保証があれば、通信不要の指示を出すことができま す。一般的に、変数であればそれがいつ書き換えられるかわからないので、HPFコンパイラー が予め通信不要であることを知るのはほとんど不可能で、そのために従来のHPFではユーザー が明確に意識できない不要な通信が頻繁に発生して、結果として並列化の効率が上がらないと 言われていました。即ち、コンパイラーは通常最悪の場合を考えてプログラムのコンパイルを 行いますから、少しでも怪しい部分があればコンパイル時には解決されずに実行時に調べられ ることになり、全ての場合を網羅して調べるために思わぬcpu間の通信が発生してしまいます。

並列化が悪くても何とかなると思ったのは、（４）の領域分割の方向を高速に切り換えることが できれば、プログラムサイズは少々増加するものの、アデトランコンパイラーの利用の時と 同様にその方法をうまく使うだけで並列化を完全に実施できると確信しているからです。従っ て、境界条件で通常並列化が困難となる部分も領域分割の方向を切り換えれば、比較的簡単に 実現できます。前述（２−３）のアデトランコンパイラーの利用では正にその機能だけを利用 して３次元MHDコードのフル並列化を実現し、256cpuで85％以上の並列化効率を得ていまし た。また、reduction節の利用で最大値、最小値、それから総和などの計算を並列化することが できます。

以上のようなHPF/JAに関する全体的な印象を持ったので、早速テストプログラムと３次元 MHDコードのVPP FortranからHPFへの書き換えを始めました。

５．VPP FortranからHPFへ（具体例）

プログラムをVPP FortranからHPF/JAへ書き換える場合、具体的にどうするのかをみるために、 ３次元の波動方程式を周期的境界条件を課してModified Leap-Frog法で解いたプログラムの具 体例を示します。これは、VPP Fortranでフルベクトル化フル並列化されていたプログラムを HPF/JAに書き換えた例で、数式の部分は省略してありますが、書き換えの部分はすべて残して あります。HPF/JAとVPP Fortranの完全なプログラムはWeb（後述のURL参照）で公開して いますので、必要な方はそちらを見て下さい。３次元波動方程式のプログラムでは、境界を除 いた３次元の格子点を（nx,ny,nz）＝（100,100,100）とし、初期条件と周期境界条件を与えて、 密度の変動分と速度の３成分の合計物理量４成分の時間発展を解いています。従って配列の大 きさとしてはf(nx2,ny2,nz2,nb)となり、z方向に領域分割して並列化します。ここで nx2=nx+2,ny2=ny+2,nz2=nz+2及びnb=4は物理量４成分になります。プログラム全体の大きな 流れとしては、1-35行がプログラムのパラメータ設定、66行が初期条件の設定、71-290行の do 500が最も外側の時間発展のdoループ、295行のclockが経過時間を測るサブルーチン、308 行のainit1が初期条件を与えるサブルーチンとなっています。

ここから、メインプログラムのHPF/JAによる並列化を前から順次より詳細にみていきます。
7-8行：プロセッサ数の指定
22-24行：配列においてz方向にデータ分割を指定
25-27行：データ分割境界のデータを両方向に１個づつ拡張して持つことを一定、後のreflect と対で用いる。
28行：データの一括転送を行うときに用いる引数の指定
29行：VPP Fortranで使用したglobal変数を用いない
30行：VPP Fortranで使用したequivalence文は使用できない
79-94行：k方向での並列化、independent文による手続き分割(k方向)の指定
79行：new( )節に指示された変数は並列処理後の値は保証されない。
81行：on home(f(:,:,k,:))には１つの分割配列を書けばよい。
81行：local( )節には、parameter文で指定された定数を除いてそのdoループの中で使わ 　　　れるすべての変数を指定する。
96-100行：データの一括転送
データ分割方向へのデータ転送がある場合、asynchronous文を用いたデータ一括転 送を用いないと処理が極端に遅くなる。 108行：分割データ境界の両側のデータのPE間での相互転送
109-121行：k方向での並列化
　　　　 　数式は省略
134-188行：k方向での並列化
　　　　 　数式は省略（first stepの計算）
193-202行：k方向での並列化
　　　　 　これはasynchronous文を用いて書き換えることができる。
204-205行：分割データ境界の両側のデータのPE間での相互転送
206-272行：k方向での並列化
　　　 　　数式は省略（second stepの計算）
282-285行：経過時間を計算して、結果を書き出す。
301行：HPFで経過時間を測る関数


初期条件を与えるサブルーチンの説明は、メインプログラムの説明にすべて含まれているので、 省略します。

VPP FortranからHPF/JAへ書き換える場合、どのようなことに注意したらよいのかをこの３次 元波動方程式のプログラムを例にとって更に述べてみます。先ず内側のdoループでベクトル化 を保ち、外側のdoループで並列化を実現するのはVPP Fortranでも HPF/JAでも同じです。ま た並列化の外側のdoループは、最終プログラムとしてはまとめられるならできるだけ広い範囲 をまとめた方が並列化の効率は上がります。最終プログラムと書いたのは理由があります。 206-273行の並列化のdoループを見て下さい。数式の内容は示していませんが、そこには多く の数式が書かれています（Webで見ることのできる完全プログラム参照）。208行のon home、 local(　　)節で多くの（このdo文に含まれる、必要な全ての）変数が指示されています。この 変数の指示はコンパイラーが良くなるに従って減少するとメーカー側からは聞いていますが、 パラメータ文で定数指示がされていない、全ての必要な変数は網羅すべきと考えておいた方が 良いと思います。そうしなければ、PE(Processing Element)間で必ず不要な通信が発生し、並列 化の効率が悪くなるのと同時に211行の最も内側のdoループがベクトル化されません。こう書 くと、local(　　)節に書くべき１個の変数を見落としたことが、どのくらい効率が悪くなるか 想像できると思います。それでは見落とし変数指示をどうして零にすることができるでしょう か。このための１つの良い方法は、doループを一度小さいdoループに分解してみて、それぞ れの小さいdoループでlocal(　　)節に書くべき変数を調べることです。小さいdoループでは 変数の個数も少ないですから比較的簡単に変数を網羅することができますし、最も内側のdo ループがベクトル化されたかどうかを調べることでlocal(　　)節の変数が網羅できたかどうか を知ることができます。こうして分解した小さいdoループがフルベクトル並列化されたことを 確認して、並列化のdoループを大きくまとめることで、効率の高いベクトル並列化のプログラ ムを完成させていくことができます。

VPP Fortranで利用していた、分割データにおける境界のデータのオーバーラップと袖の機能は、 HPF/JAではshadowとreflect文の利用で対応することができます。また、HPFでは、VPP Fortran であったようなglobal変数の概念はありませんから、equivalenceでglobal変数とlocal変数を結 びつけることはできません。従ってlocal変数を用いて計算して、global変数を用いてreadと write を行うことはできませんが、適当な配列を用意して、local変数との間でasynchronousな どを用いて一括転送を行い、その用意した配列をreadとwriteに使用することによって高速に 入出力を行うことができます。もちろん、その用意した配列の分だけプログラム容量は各PE で増えますので、VPP5000のメモリにプログラムが納まる範囲内までの配列しか取ることがで きません。また、分割データを直接writeすることはできますが、readする場合にはdo形反復 に一部制限があります。しかし、いずれの場合も大きな配列の分割データをdo形反復で直接出 力に用いることは非常に長時間を要しますので、前に述べたように別の配列を用意して、その 配列をreadとwriteに、それもできれば書式なし入出力文を利用することをお奨めします。分 割した方向におけるデータの転送には、必ずasynchronous文などによる一括転送を用いるべき です。そうしないと、極端に長い処理時間がかかることになります。 もう１つasynchronous文の使い方にデータ分割の方向を変えて処理する方法があります。

　!hpf$  distribute fa(*,*,block,*) onto pe
　!hpf$  distribute fb(*,block,*,*) onto pe

２つのデータ分割の方向を変えた配列faと（z方向に分割）とfb（y方向に分割）を用意し、 faを用いて、z方向を除いて計算を行い、asynchronous文を用いてfaからfbに一括転送し、fb でz方向の計算を行って再びasynchronous文を用いてfbからfaに一括転送することで計算を完 了させる方法です。例えば、96-100行のz方向の境界条件の計算ではこのデータ分割の方向を 変えた配列を利用することができ、その応用は非常に広いのと同時に処理時間もあまり増えま せん。asynchronous文を入出力文やデータ分割の方向を切り換えるために用いた例はWebでみ ることができます。

VPP FortranからHPF/JAに書き換えたプログラムはlocal(　　)節を除けば、むしろより簡明に なり、見やすくなっています。そして、フルベクトル化フル並列化ができていますので処理速 度もVPP Fortranと比べて全く遜色がありません。またプログラムサイズもほとんど増えていま せん。特に大きな利点は、VPP Fortranでフルベクトル化フル並列化できているプログラムを HPF/JAに書き換える時、プログラムの構造をほとんど変える必要がないということです。この 利点のために地球磁気圏の３次元グローバルMHDコードのような大きなプログラムの書き換 えでは労力が大幅に軽減できます。

- VPP FortranからHPFへ（具体例）の補足説明 -

この章のこれまでの内容について富士通株式会社の方から補足説明を頂きましたのでそれらの 重要点をまとめて以下に記述します。

(1) local節の記述方法について local節には (...) を省略する記述方法があります。その場合81行のon homeは次のように記述 され、

    81  !hpf$  on home(f(:,:,k,:)),local begin

local( ... )の ( ... )を取り除くとbegin 〜end onまでの間の変数、配列の全てを記述した ことと同じになります。VPP Fortranのspread doで並列化したことのあるdoループなどの場合 に利用すると便利です。即ち、プログラマーがプログラムの内容を完全に掌握していて、ルー プ内での通信は絶対に必要ないと分かっている場合は、local( ... )の ( ... )を省略できま す。しかし、一般的には一つづ確認しながらnew節、local節に記述する方が安全で確実です。

(2) local指示をどう作るかについて local指示をどう作るかですが、2通りのアプローチが提案されています。「パラメータ文で定 数指示がされていない、全ての必要な変数は網羅すべき」場合であれば方法１で十分と考えら れます。一部の変数はループ内で通信が必要な場合には方法２がお勧めです。

方法１：変数なしlocal指示を使用する。

上記(1)の説明の通りです。ただしこの記述方法が使用できるのは、プログラマーがプログラ ムの内容を理解し、並列ループ内での通信は必要ないと断定できる場合に限ります。もし誤っ て使用した場合には、正しい動作は保証されません。これは VPP Fortran で言えば、プロセ ッサ上にないデータをローカル変数名で参照した場合に相当します。

方法２：まずnew指示を完全にし、コンパイラの自動local化に任せ、漏れたものだけをlocal 指示する。

HPF L00061（８月から使用できるようになったはずです）では、on homeで指定した変数や、 その変数と分割の一致する変数については、かなりの場合自動的にlocal化されるようになり ました。（ただし、shadow部分のアクセスについては自動ではlocal化されませんので、local 指示が必要です。）また、independent指示文のnewで指定された変数についても、local指示 なしでも自動的にlocal化されます。L00061では、new指示の漏れについて以下のメッセージ で検出するようになりました。

joo7102i-i
英文メッセージ
A scalar variable appears on the left side of an assignment statement within DO loop, but does not appear in the NEW clause nor the REDUCTION clause. (name: 変数名 ) 和文メッセージ
DOループ内にNEW節にもREDUCTION節にも現れていないスカラ変数に対する代入文があ ります. (name: 変数名 )

方法２の実際の手順は、以下のようになります。

1) new指示を分かる範囲で書き、local指示は書かないで、コンパイラを通してみる。
　 → new指示の漏れがあると上記メッセージが出る。
2) メッセージを参考にnew指示を修正し、メッセージが出なくなるようにする。
3) リスタ出力(-Wh,-Ltオプション)を参考に、自動的にlocal化されていない変数アクセス について、通信なしでよいと判断できるならlocal指示を加える。この手間が大きい場 合には、前述のdoループを徐々に融合していく方法を採用することが考えられる。

(3) asynchronous指示文の使用について
193-202行：k方向での並列化について、asynchronous指示文を用いて書き換えることもできる と書いていますが、プロセッサ間通信の必要のない代入文なので、independent+onを使う方が 高速です。independent+onではベクトルのロード／ストアになり、asynchronous使う場合には 実行時ライブラリの呼出しになります。

(4) readとwriteについて
read文がdo形反復を含む場合、分割配列を入力並びに書くことはできません。分割配列がdo型 反復の外に書かれている場合でも不可です。 大きな配列の分割データを入出力する場合は、前述のようにdo形反復を利用するので別の配列 を用意してその配列をreadとwriteに利用するのが効率的ですが、分割データがあまり大きくな い場合は次のように変数名だけで指定して出力する方法や書式なし出力文

	write(*,100) a,b,c
	write(10) a,b,c

を用いた方が記述も簡単でかつ出力時間もそれほどかかリません（asynchronousコピーと重複 配列の出力の時間を足したものと同程度の時間がかかるだけです）。 また、新たな配列を用意しないで、次のように大きな配列の分割データの一部を書式なし入出 力文を用いてreadとwriteをすることができます。

       do 174 m=1,nb
       do 174 k=1,nz2
       read(10) f(1:nx2,1:ny2,k,m)
  174  continue

       do 176 m=1,nb
       do 176 k=1,nz2
       write(ntap) f(1:nx2,1:ny2,k,m)
  176  continue

この方法を用いれば、大きな配列の分割データを入出力する場合でもプログラムサイズは増えず、 かつ新たな配列を用意してasynchronousコピー使う場合と同程度の時間しかかかりません。こう して、global変数とlocal変数を結びつけるequivalence文を使用できないHPF/JAのプログラムも VPP Fortranのプログラムと同じサイズにすることができます。

６．HPFとVPP Fortranでの計算速度の比較

プログラムのVPP FortranからHPFへの書き換えができたら、コンパイルと実行を行うことに なります。その例をVPP FortranとHPFの場合について示しますが、詳しくは参考文献の「新 スーパーコンピュータVPP5000/56の利用について」の説明を見て下さい。

先ず、VPP Fortranプログラムのコンパイルと実行に用いているqsubとシェルスクリプトの例 をプログラムは全てユーザーアカウントの直下のtestデレクトリィにあるとして示します。VPP Fortranプログラムは、prog90.fでそのコンパイル結果のソースリストはprog90listに出力されま す。

  qsub -q x -eo -o pcomp902.out pcomp902.sh
  qsub -q z -eo -lP 2 -o pexec90.out pexec90.sh
  <>
    cd test
    frt -Wx,-Lt prog90.f -Pdt -o prog90 -Z prog90list
  <>
    #  @$ -lt 1:30:00
    #  @$-q z  -eo
    cd test
    timex prog90

HPFプログラムのコンパイルと実行に用いているqsubとシェルスクリプトの例を次に示しま す。HPFプログラムは、proghpf.fでそのコンパイル結果のソースリストはhpflistに出力されま す。

  qsub -q x -eo -o pconphpf2.out pcomphpf2.sh
  qsub -q z -eo -lPv 2 -o pexechpf.out pexechpf.sh
  <>
    cd test
    frt -Wh,-Lt proghpf.f -Pdt -o proghpf -Z hpflist
  <>
    #  @$ -lt 1:30:00
    #  @$-q z  -eo
    cd test
    timex proghpf
    

VPP FortranからHPF/JAへ書き換えたプログラムが効率よく動くかどうかは、コンパイル結果
のソースリストhpflistの最も内側のdoループが全てベクトル化できたかどうかで先ず知ること
ができます。富士通提供のサンプラを使うのも有効な方法です。サンプラを用いるとPE間の
通信がどれだけ発生しているか知ることもできます。しかし、ベクトル化と並列化によってプ
ログラムの計算速度がどれだけ速くなったかを知るには、cpu時間を実際に計測する必要があ
ります。前の章で示した３次元の波動方程式のプログラムでは、計算のメインdoループの計算
時間を計測するようになっています。その３次元の波動方程式のプログラムは、VPP Fortranと
HPF/JAのいずれの場合もフルベクトル化とフル並列化ができていて、完全なプログラムはWeb
で見ることができます。

前にも書きましたように私達の本来の目的は、太陽風と地球磁気圏の相互作用を調べる３次元
グローバルMHDシミュレーションコードをVPP FortranからHPF/JAに書き換えて、その計算
速度を比較し、HPF/JAのプログラムでVPP Fortranの場合の50％以上の計算速度を出すことで
した。太陽風地球磁気圏相互作用の３次元MHDシミュレーションでは、MHD方程式と
Maxwell's equationを差分化してModified Leap-Frog Methodで初期値境界値問題として解き、太
陽風と惑星間磁場の変化に伴う地球磁気圏の応答を調べています。このために、外側の境界を
できるだけ遠ざけて境界の影響を少なくすると共に磁気圏尾部の計算領域を広げて遠尾部構造
を調べ、かつ空間分解能をできるだけ上げようとしていますので、どうしても３次元格子点の
数をシミュレーションの限界まで増やすことになります。そのMHDシミュレーションの例を
第１図と第２図に示します。第１図は、惑星間磁場が北向きでかつ夕向きの場合の定常状態の
地球磁気圏の３次元磁力線構造を示しています。磁気圏境界の磁気リコネクションが北半球で
は夕側、南半球では朝側の高緯度尾部領域で起こるために朝夕非対称の磁力線構造になってい
ます。第２図は、太陽と地球の間で太陽風と惑星間磁場を観測しているACE衛星の１分毎の
データを計算のインプットとして用いた場合に、３次元グローバルMHDシミュレーションか
ら得られた地球磁気圏のスナップショットの例を示しています。これは、国際宇宙天気研究プ
ロジェクトの一環として、時々刻々の太陽地球環境の変化を知るための数値モデル開発基礎研
究の一翼を担うものです。このような地球磁気圏の３次元グローバルMHDシミュレーション
は、VPP5000/56でも計算をつないで合計50〜300時間の計算時間を必要とします。

VPP FortranとHPF/JAで書かれた３次元グローバルMHDコードをFujitsu VPP5000/56で実行し
た場合の計算速度の比較を表１に示します。いずれのプログラムもフルベクトル化とフル並列
化が実施されていて、現在までの様々なテスト計算の結果を集めて最も効率が出せるようにし
て作った実際のMHDシミュレーションを行うプログラムを用いています。表には、PE数、境
界を除いた３次元の格子点、１回の時間ステップを進めるのに要するcpu time (sec)、計算速度
(Gflops)、１PE当たりの計算速度(Gflops/PE)をVPP FortranとHPF/JAの場合を並べて載せてい
ます。ちなみに使用した名大大計センターのFujitsu VPP5000/56の特徴は以下のとおりです。

　名大大計センターのFujitsu VPP5000/56の特徴
　　　・PE台数：　　　　　　　　　　　　　　　56PE
　　　・PEの理論最大性能：　　　　　　　　　　9.6Gflops
　　　・1PEの最大ユーザー利用可能メモリ容量： 7.5GB


アムダールの法則から知られているように、ベクトル化や並列化の効率を上げるためには、全
ての内側のdoループをベクトル化して、かつ全ての外側のdoループを並列化することが極め
て重要になります。１つでもベクトル化や並列化できないdoループが残っているとその効率は
著しく低下します。実際に行った人なら、最後の1個のdoループをベクトル化／並列化した時、
飛躍的に（数倍になることもあります）計算効率が上がった経験を持っていられることと思い
ます。配列の格子点数はnx*ny*nzで示し、並列化のデータ分割はz方向に行っていますので、
両側の境界を加えるとz方向のデータ配列の大きさはnz2=nz+2になります。従って、並列化の
効率をより上げるためにはz方向のデータ配列の大きさnz2を使用するPE数の整数倍にする
のがより望ましいことになります。また、仕事量を大きくする、即ち格子点数nx*ny*nzを増加
してプログラムサイズを大きくすると、一般にベクトル化／並列化の効率がよくなります。ま
た、比較のためにベクトル化をしないscalarモードとベクトル化のみ働く1PEでの計算速度も
示しています。このような観点に基づいて、表１の結果をみると得られた計算速度の理解が容
易になると思います。

表からすぐに気がつくことは、scalarモードと1PEのベクトルモードとの比較からベクトル化
だけで約40倍の高速化が得られ、全体的にVPP FortranとHPF/JAの計算速度がほぼ同程度で
あり、かつ増やしたPE数に比例して計算速度も増加するスケーラビリティーは、大体の傾向
としては成り立っていることです。もう少し詳しくみると、従来からのMHDシミュレーショ
ンではnx*ny*nz=500*100*200の格子点をよく用いていましたが、nz2=202=2*101となって、PE
数の整数倍にはなり難く並列化の効率を上げるという点からはあまり適切な格子点数の選択で
はなかったことが分かります。これは、PE数を増やした時、計算速度は上がってはいるものの、
その上がり方は比例して上がるよりかなり低下していることから知ることができます。

そのスケーラビリティーは、格子点数nx*ny*nz=200*100*478 (nz2=480=2**5*3*5)及び
800*200*670 (nz2=672=2**5*3*7)ではかなりよくなっていることが分かります。また、多数の格
子点nx*ny*nz=1000*500*1118及び1000*1000*1118 (nz2=1120=2**5*5*7)を用いた計算では、nz2
がPE数の整数倍であり、かつプログラムサイズが大きくなった効果もあって、VPP Fortranと
HPF/JAの両方とも56PEの計算では400Gflops以上の効率の高い計算速度が得られています。
これは、Fujitsu VPP5000/56の56PEを用いた時、HPF/JAで書き換えた３次元グローバルMHD
コードが、ベクトル化／並列化の総合絶対効率として最大76.5％に達したことを示しています。

７．並列計算の将来

cpuのクロック時間が限界に近づいてきている今日、更なる高速計算を実現しようとする時、
並列計算は必要不可欠なものです。その並列計算法に関して常日頃問題に思っている私見をこ
こに述べてみます。その第１点は、ベクトル計算機の利用の時までに全てのプログラムの形を
作っていた１次元配列のままで並列計算ができないかということです。VPP FortranやHPF
では配列のデータ分割指示文のところで最初から拒否されてしまいます。元来フォートランプ
ログラムで

  f(i,j,k,m); i=1,nx, j=1,ny, k=1,nz, m=nb
  f(i1); i1=i+nx*(j-1)+nx*ny*(k-1)+nx*ny*nz*(m-1)


の２つの表現は同じであると教わりました。１次元配列のメリットは目的の番地のメモリを
指示するのに、四則演算の回数を減らせるので、より計算の高効率化ができることです。配列
の大きさを決めるnx,ny,nz,nbが定数として明確に定義され、１次元配列のセパレーションとし
て機能するならば、１次元配列のままでデータ分割ができてもよいように思われます。ちなみ
にHPFの次の分割データの指示文は

  f(*,*,block,*)
１次元配列を
  f(i1):i1=i+nx*(j-1)+nx*ny*(m-1)+nx*ny*nb*(k-1)
  

として、k=1,k1,k=k1+1,k2,…,k=kn-1+1,nzの範囲にn分割するのと同等にすることができます。た
だし、kの取り得るそれぞれの範囲には1+nz/nまたはnz/nになります。次に、分割データの境
界のデータを隣り合うPEで２重に持つHPFのshadow/reflectの機能ですが、これも、nx*ny*nb
の大きさの隣り合うデータを持ち合うだけですから、それを可能にする一括データ転送の機能
があれば解決できると考えられます。こうして、１次元配列データでもフルベクトル並列化又
はフル並列化が何らかの方法でできないかとここ数年間その実現の期待を持ち続けています。

第２点は、物理的３次元空間を考えた時、領域分割の本質は何か、そしてそれがデータ分割と
どのように関係を持つ必要があるかということです。CRAY-1は、主メモリが１MWと小さか
ったために、折角高速計算ができるにもかかわらず、多くの格子点を用いる大きなプログラム
を動かすことができませんでした。そのためにCRAY-XMP、CRAY-YMPが開発されたのです
が、多くの格子点を持つ大容量のプログラムをベクトル多重化（後では自動並列化）して動か
す方法として領域分割（３次元でz方向に分割）を用い、実際に8cpuで予想通りの高速計算結
果を得ることができました。私達の１次元配列コードを用いたのですが、ほとんど３次元
MHDコードを書き換える必要は生じませんでした。即ち、空間３次元をスライスしてある方
向に領域分割するのは、原理的に簡単で大きな困難は生じないだろうとの結論を得ました。

しかし、それと同時に３次元空間をスライスする領域分割法の限界のことも考えていました。
それは、格子点が非常に多くなった時、あるいはプロセッサ数が非常に多くなった超並列計算
では問題が生じるのは明らかだったからです。例えば、３次元格子点が(nx,ny,nz)=(100,100,100)
だった時、PEが101個以上だとするとz方向の領域分割だけではcpuが余ってしまいますし、
また領域分割の両側の分割データを隣り合うcpuが二重に持つとしますと、cpuのメモリの1/3
しか有効に使えないことになり、隣り合うcpu間の通信量も大幅に増加してしまいます。

これらの問題は、３次元空間をスライスする領域分割では解決できないので、もう１次元領域
分割の方向を増すか、３次元空間を直方体の集合に分割するような領域分割を用いて解決する
方法が考えられます。即ち、３次元空間を１つのcpuのメモリにおさまるような直方体のセル
に空間領域分割を行い、周囲の６面を境界条件として隣り合うcpuとの間でデータを二重に持
てば、立体空間の領域分割として並列計算ができると考えられるからです。もちろん、１回の
時間ステップを進める毎に、境界条件のデータは更新する必要があり、次の時間ステップの計
算を始める前にcpu間で通信して予め必要な境界のデータを１つのcpuに集める必要がありま
す。MHDの偏微分方程式や流体の偏微分方程式などを差分化して時間発展を解く場合、解く
べき領域の１面外側のデータ（隣のcpuにある）を予め６面集めれば、時間ステップを１回進
める間は全く他のcpuとの間で通信を発生させずに計算を実行することができますので、高速
並列計算が実現できます。その６面のデータは、各cpuにとっては境界条件の扱いになり、こ
れまでのプログラムで通常扱ってきたものとなんら変わるものではありません。隣り合うcpu
間で境界のデータを転送するために生じる通信量も、cpuが増えた場合並列化して通信できれ
ば実際の通信時間は増加しないことになります。即ち、３次元空間を直方体のセルに分割して
計算する方法は、MHDモデル等の差分解法による初期値境界値問題では、非常に望ましい計
算法になり得るということができます。

８．MHDシミュレーション発展

２次元と３次元のグローバルMHDシミュレーションから太陽風と地球磁気圏の相互作用を約
２０年間調べてきましたが、その研究の歴史を振り返ると同時に、実際のシミュレーション
実行のデータを基にして将来を展望してみたいと思います。先ず、各年度でどのくらいの２
次元と３次元の総格子点数を用いたMHDシミュレーションを実行することができたかを示し
たのが第３図です。同図で●▲印は他の研究グループのMHDシミュレーション実行例、○△
×は私達の実行例を示します。私達は1976年にトロイダルプラズマの３次元MHD不安定解析
のために３次元MHDシミュレーションを始め、使用できた格子点数も境界を加えて18*18*12
及び26*26*18と小さなものでした。しかし、計算時間はFujitsu M-100 及びM-200を用いて
50〜200時間が必要で、大型MHDシミュレーションの実行は計算速度によって制限されていま
した。

1982年からCRAY-1を用いた地球磁気圏の３次元グローバルMHDシミュレーションに着手し、
使用格子点数は、62*32*32、50*50*26と増大しましたが、CRAY-1では主メモリの大きさ
（1MW）が制限となっていました。この頃100の３乗の空間格子点を用いたMHDシミュレー
ションを実行できたらどんなにすばらしいことができるだろうと考えたものでした。その
後、国産のベクトル機が開発されて、Fujitsu VP-100, VP-200及びVP-2600の使用によって
100の３乗を越す３次元MHDシミュレーションを行うことができるようになりました。更に
1995年からベクトル並列機VPP500が、そして2000年にVPP5000が使用できるようになり、
1000の３乗の格子点数を用いた３次元MHDシミュレーションも現実のものとなりました。
利用できる総格子点数からみたMHDシミュレーションの年毎の発展拡大の傾向をみますと、
この24年間破線で示すように、平均して３年で４倍あるいは５年で１０倍の割合で増大し
てきています。

第４図にMHDシミュレーションを実行する場合、時間ステップを１回だけ進めるのに必要
とする計算時間を、第５図にMHDシミュレーションで用いる総格子点数に対する計算機メ
モリの必要総容量を示します。計算機の種類と計算方法にも依存しますが、それが決まる
と総格子点数に比例して計算時間は長くなり必要計算機メモリは増大します。地球磁気圏
のMHDシミュレーションでは、数千回から数万回の時間ステップを繰り返し計算しますか
ら、１回の時間ステップを進めるのに必要とするcpu timeは大体10秒以下でないと実用に
なりません。こうして、計算機の計算速度から使用できる３次元格子点数の大きさが決ま
ってきます。同時に、その３次元格子点数を用いた時の３次元MHDプログラムサイズが計
算機主メモリに収まる必要があります。

MHDモデルでは、密度、速度３成分、圧力、磁場３成分の８成分の時間発展を解きますか
ら、総格子点数をNgとした時独立な変数は８Ngで、容量は単精度で32Ng Bytesになりま
す。MHDモデルの必要不可欠の変数量は 2 step Lax-Wendrotf法（2LW）を用いた場合32Ng 
Bytesが、Modified Leap-Frog法（MLF）を用いた場合、２つの時刻の変数を同時に保持
するために64Ng Bytesが必要です。CRAY-1で開発して多種類のベクトル機で実行してい
た3次元MHDコード(2LW, MLF)では、作業変数を徹底的に減らした結果、プログラムサイ
ズは必要不可欠な変数量の1.3倍以内に収まっていました。並列機になって、ベクトル機
の場合と同じように作業変数を減らすことはほとんど不可能になり、ADENART、VPP500、
VPP5000ではプログラムサイズがかなり増大しています。しかし、VPP500やVPP5000では
その増加の割合も、MLF法使用の場合で必要不可欠の変数量の約2.5倍(160Ng=2.5*64Ng)
に収まっています。こうした総格子点数に対する計算速度と必要計算機メモリの関係を
調べておくと、シミュレーションで解明するテーマを設定した時、どれだけ大きな3次元
MHDシミュレーションを実行することができるか、そして将来の見通しはどうかをかなり
正確に予測することができます。

MHDシミュレーションで用いる総格子点数に対する種々の計算機による絶対処理能力の比
較をまとめたものを第６図に示します。ベクトル機からベクトル並列機への移行で計算
機の処理能力は飛躍的に高くなりました。最近の高性能のワークステーションやパーソナ
ルコンピュータは、一昔前のスーパーコンピュータの計算速度（0.3 Gflops程度）に達し
ています。しかし、最大級のスーパーコンピュータは常にそれらの約千倍の高速計算処理
能力を有しています。多数の3次元格子点を用いた大容量シミュレーションでは、メモリ
不足やRISC計算機でのcache不足などによって、更に決定的な差が生じます。フルベクト
ル化／フル並列化されたプログラムを用いて、VPP5000等のスーパーコンピュータで100時
間かかる3次元MHDシミュレーションは、スーパーコンピュータ以外の通常のワークステー
ションやパーソナルコンピュータで実行することは全く不可能です。地球磁気圏の3次元
グローバルMHDシミュレーションなどの大型コンピュータシミュレーションを実行するた
めには、最大級のスーパーコンピュータの利用は常に必須の条件です。

９．おわりに

富士通VP-2600、日立S820、NEC SX-3、CRAY Y-MPなどのベクトル計算機の時代までは、フ
ルベクトル化された３次元MHDコードを用いてFortranコンパイラーが載っている全ての計算
機を利用して、太陽風と地球磁気圏相互作用の３次元グローバルMHDシミュレーションを実
行することができました。このFortranプログラムの汎用性のために、私達は３次元MHDコー
ドを世界中どこででも動かすことができ、MHDコードの配布などを通して世界中の多くの研
究者と共同研究を行うことができました。しかし、ベクトル並列機と超並列機がコンピュータ
シミュレーションの世の中に現れてくるや否や、Fortranプログラムの並列化の効率を上げるた
めにコンピュータに依存した様々の異なった手法を採らなければならなくなりました。多くの
シミュレーション研究者は共通のプログラム言語を失い、特定のメーカーの機種でしか並列化
の効率を上げられない、方言のFortranプログラム言語を使わざるを得なくなりました。Fujitsu 
VPP Fortranもその１つです。即ち、VPP Fortranでフルベクトル化フル並列化されているプログ
ラムも、日立やNECのスーパーコンピュータではその威力を発揮できませんし、欧米で主流
となっている超並列計算機ではなおさらです。
私達シミュレーション研究者は、スーパーコンピュータの世界でどのスーパーコンピュータで
も効率よく動く共通プログラム言語が１日も早く再び使えるようになることを強く望んでいま
す。その共通プログラム言語の候補としてHigh Performance Fortran (HPF)とMessage Passing 
Interface (MPI)があると言われてきました。HPFは、米国の共同研究者が共通プログラム言語と
して優れていると言っている反面、多くの大型プログラムで性能が十分に出ないという批判の
声も出ていました。

こうした中、この６月から名大大計センターのFujitsu VPP5000/56でHPF/JA(JAHPFによるHPF
の日本拡張版）が使えるようになりました。私達はすぐに富士通提供のサンプルプログラムの
テスト、私達が準備していた３次元波動方程式などのテストプログラム、そして３次元MHD
コードのVPP FortranからHPF/JAへの書き換えを始めました。書き換えの具体的内容はこれま
でみてきたとおりです。書き換えは意外とスムーズに進み、全部のプログラムの書き換えとテ
ストを３週間で行うことができ、MHDコードの計算結果の最終確認には更に２週間が必要で
したが、全プログラムの書き換えを無事に完了することができました。VPP FortranからHPF/JA
へのプログラムの書き換えがスムーズに行えた最大の理由は、並列化するための形式に基本的
な共通性があったためです。こうして、56PEまでの多数のProcessing Elementsを用いた場合も
含めて、HPF/JAに書き換えたプログラムは、VPP Fortranと同等の性能を発揮できるようにな
りました。

今回の書き換えから得た１つの結論は、VPP Fortranでフルベクトル化フル並列化されている流
体コードやMHDコードは、比較的容易にHPF/JAに書き換えることができて、更にそのHPF/JA
のプログラムはVPP Fortranと同等の性能を得ることが期待できるということです。３次元
MHDコードによる最大性能はFujitsu VPP5000/56 32PEで230Gflops以上、56PEで400Gflops
以上を実現しています。私達はHPF/JAで書き換えた３次元MHDコードが、日立のSR8000
やNECのSX-5のようなスーパーコンピュータで近い将来効率よく動くようになるのを
とても強く期待しています。そして、コンピュータシミュレーションの世の中で、HPF/JAが世
界中の共通プログラム言語の１つになれば、大変すばらしいと思っています。さらに、これは
地球磁気圏の３次元グローバルMHDコードを国内外の共同研究者に提供して、シミュレーシ
ョン共同研究を新しいしかも高度な次元で再出発させるためにはどうしても必要な基本条件と
なります。日本のスーパーコンピュータのメーカー３社の開発に大いに期待しています。

本稿で述べてきた内容を補助するため、HPF/JAとVPP Fortranで書かれた３次元波動方程式の
完全なプログラムや境界条件を含む並列化の主要部分を示した３次元グローバルMHDコード
を次のURLで公開していますので、必要に応じてご覧頂けると幸いです。

　http://gedas.stelab.nagoya-u.ac.jp/simulation/hpfja/hpf01.html

謝辞　

本稿のコンピュータシミュレーションは名古屋大学大型計算機センターを利用してなされたも
のです。また、VPP FortranからHPF/JAへの書き換えでは多くの助言を頂いた名古屋大学大型
計算機センターの津田知子助手と富士通株式会社の方に感謝いたします。

参考文献

荻野竜樹、「太陽風と磁気圏相互作用の電磁流体力学的シミュレーション」，
プラズマ・核融合学会誌，CD-ROM 特別企画（解説論文），
Vol.75, No.5, CD-ROM 20-30, 1999.

荻野竜樹、「太陽風磁気圏相互作用の計算機シミュレーション」，
名古屋大学大型計算機センターニュース, Vol.28, No.4, 280-291, 1997.

津田知子、「新スーパーコンピュータVPP5000/56の利用について」、名古屋大学大型計算機セ
ンターニュース、Vol.31, No.1, 18-33, 2000.

High Performance Fortran Forum, 「High Performance Fortran 2.0 公式マニュアル」、
Springer, 1999.

富士通株式会社、「HPFプログラミング　〜並列プログラミング（HPF編）〜」、UXP/V 第1.0
版、2000年6月.

富士通株式会社、「VPP FORTRANプログラミング」、UXP/V, UXP/M 第1.2版、1997年4月.

プログラム１．　HPF/JAに書き換えた３次元波動方程式のプログラムの例。HPFの部分は全
て示していますが、数式の部分は省略されています。完全なプログラムはWebで公開していま
す。
     1  c
     2         program hwave3
     3  c
     4  cc     hwave3.f  main program of 3-dimensional wave equation
     5  c      modified leap-frog method  2000.06.17
     6         implicit real*8 (a-h,o-z)
     7         parameter(npe=2)
     8  !hpf$  processors pe(npe)
     9         parameter(nx=100,ny=100,nz=100)
    10         parameter(nb=4,iip0=8,iiq0=1,iir0=1,last=4)
    11         parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
    12         parameter(nz1=nz+1,nz2=nz+2)
    13         parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
    14         parameter(n4=n3*nb)
    15         parameter(thx=0.5,tam=0.5,vis0=0.10)
    16  c
    17         dimension f(nx2,ny2,nz2,nb),u(nx2,ny2,nz2,nb)
    18  c      dimension gf(nx2,ny2,nz2,nb),gu(nx2,ny2,nz2,nb)
    19         dimension v(nx2,ny2,nz2,nb),p(nb)
    20         dimension ppin(10)
    21  c      real*8 zt0,zt1,zt2,zt
    22  !hpf$  distribute f(*,*,block,*) onto pe
    23  !hpf$  distribute u(*,*,block,*) onto pe
    24  !hpf$  distribute v(*,*,block,*) onto pe
    25  !hpf$  shadow     f(0,0,1:1,0)
    26  !hpf$  shadow     u(0,0,1:1,0)
    27  !hpf$  shadow     v(0,0,1:1,0)
    28  !hpf$  asyncid id1
    29  c      !hpf$ global gf,gu
    30  c      equivalence (gf,f),(gu,u)
    31         common /blk/f
    32  c
    33  c           ppin xl   yl   zl   dxl  dyl  dzl  dn  dv   9  10
    34         data ppin/62.0,62.0,62.0,10.0,10.0,10.0,1.0,0.0, 9, 10/
    35  c

    58  c      initial parameters
    59  c
    60         write(6,122) nx,ny,nz,last,nx2,ny2,nz2
    61         write(6,124) t1,hx,hy,hz,dx1,dy1,dz1
    62    122  format(1h ,10i8)
    63    124  format(1h ,8(1pe10.3))
    64  c
    65  c      initial conditions
    66         call ainit1(ppin)
    67  c
    68  c      start of calculation
    69         call clock(zt0)
    70  c
    71         do 500 ii=1,last
    72         call clock(zt1)
    73         do 300 iir=1,iir0
    74         do 200 iiq=1,iiq0
    75         do 100 iip=1,iip0
    76  c
    77  c      boundary conditions
    78  c
    79  !hpf$  independent,new(i,j,k)
    80         do 22 k=2,nz1
    81  !hpf$  on home(f(:,:,k,:)),local(f,i,j,m) begin
    82  c
    83         do 222 m=1,nb
    84         do 221 j=2,ny1
    85         f(1,j,k,m)=f(nx1,j,k,m)
    86         f(nx2,j,k,m)=f(2,j,k,m)
    87    221  continue
    88  c
    89         do 222 i=1,nx2
    90         f(i,1,k,m)=f(i,ny1,k,m)
    91         f(i,ny2,k,m)=f(i,2,k,m)
    92    222  continue
    93  !hpf$  end on
    94     22  continue
    95  c
    96  !hpf$  asynchronous(id1),nobuffer begin
    97         f(1:nx2,1:ny2,  1,1:nb) = f(1:nx2,1:ny2,nz1,1:nb)
    98         f(1:nx2,1:ny2,nz2,1:nb) = f(1:nx2,1:ny2,  2,1:nb)
    99  !hpf$  end asynchronous
   100  !hpf$  asyncwait(id1)
   101  c
   102  c      case of iip=1  2-step Lax-Wendroff method
   103         if(iip.eq.1) then
   104         t=0.5*t1
   105         dx=0.5*dx1
   106         dy=0.5*dy1
   107         dz=0.5*dz1
   108  !hpf$  reflect f
   109  !hpf$  independent,new(i,j)
   110         do 32 k=1,nz1
   111  !hpf$  on home(u(:,:,k,:)),local(f,u,i,j,m) begin
   112         do 321 m=1,nb
   113         do 321 j=1,ny1
   114         do 321 i=1,nx1
   115         u(i,j,k,m)=0.125*(f(i,j,k,m)+f(i+1,j,k,m)+

   119    321  continue
   120  !hpf$  end on
   121     32  continue
   122  c
   123         else
   124         t=t1
   125         dx=dx1
   126         dy=dy1
   127         dz=dz1
   128         endif
   129  c
   130  c
   131  c      first step
   132  c
   133  c      !hpf$  reflect f
   134  !hpf$  independent,new(i,j)
   135         do 40 k=1,nz1
   136  !hpf$  on home(u(:,:,k,:)),local(f,u,p,i,j,m,dx,dy,dz) begin
   137         do 401 j=1,ny1
   138         do 401 i=1,nx1
   139  c
   140         p(1)=0.125*(f(i,j,k,1)+f(i+1,j,k,1)+

   186    401  continue
   187  !hpf$  end on
   188     40  continue
   189  c
   190  c      preparation of second step
   191  c      second step
   192  c
   193  !hpf$  independent,new(i,j,m)
   194         do 62 k=1,nz2
   195  !hpf$  on home(v(:,:,k,:)),local(f,v,m,i,j) begin
   196         do 621 m=1,nb
   197         do 621 j=1,ny2
   198         do 621 i=1,nx2
   199         v(i,j,k,m)=f(i,j,k,m)
   200    621  continue
   201  !hpf$  end on
   202     62  continue
   203  c
   204  !hpf$  reflect u
   205  !hpf$  reflect v
   206  !hpf$  independent,new(i,j,k)
   207         do 60 k=2,nz1
   208  !hpf$  on home(f(:,:,k,:)),local(f,u,v,p,i,j,m,dx,dy,dz,
   209  !hpf$*    dx1,dy1,dz1,dx2,dy2,dz2) begin
   210         do 601 j=2,ny1
   211         do 601 i=2,nx1
   212  c
   213         p(1)=0.125*(u(i,j,k,1)+u(i-1,j,k,1)+

   271    601  continue
   272  !hpf$  end on
   273     60  continue
   274  c
   275  c      end of 1 time step advance
   276  c
   277    100  continue
   278    200  continue
   279    300  continue
   280  c
   281         call clock(zt2)
   282         zt1=zt1-zt0
   283         zt2=zt2-zt0
   284         zt=zt2-zt1
   285         write(6,402) ii,zt0,zt1,zt2,zt
   286    402  format(1h , i6,1pe12.3,3(0pf12.5))
   287  c
   288  c      write the output data
   289  c
   290    500  continue
   291      9  continue
   292  c
   293         stop
   294         end
   295         subroutine clock(ti)
   296         implicit real*8 (a-h,o-z)
   297  c      real*8 ti,ti1
   298         ti=1.0d0
   299         ti1=1.0d0
   300  c      call gettod(ti1)
   301         call fjhpf_gettod(ti1)
   302         ti=1.0d-6*ti1
   303  c      x=0.0
   304  c      y=secnds(x)
   305  c      ti=1.0d0*y
   306         return
   307         end
   308         subroutine ainit1(ppin)
   309         implicit real*8 (a-h,o-z)
   310         parameter(npe=2)
   311  !hpf$  processors pe(npe)
   312         parameter(nx=100,ny=100,nz=100)
   313         parameter(nb=4,iip0=8,iiq0=1,iir0=1,last=4)
   314         parameter(nx1=nx+1,nx2=nx+2,ny1=ny+1,ny2=ny+2)
   315         parameter(nz1=nz+1,nz2=nz+2)
   316         parameter(n1=nx2,n2=n1*ny2,n3=n2*nz2,noinp=30)
   317         parameter(n4=n3*nb)
   318  c
   319         dimension f(nx2,ny2,nz2,nb)
   320  c      dimension gf(nx2,ny2,nz2,nb)
   321         dimension ppin(10)
   322  !hpf$  distribute f(*,*,block,*) onto pe
   323  !hpf$  shadow     f(0,0,1:1,0)
   324  c      equivalence (gf,f)
   325         common /blk/f
   326  c

   339  !hpf$  independent,new(i,j,k)
   340         do 10 k=1,nz2
   341  !hpf$  on home(f(:,:,k,:)),local(f,i,j,m,x,y,z,dn1,dv1,ax1,
   342  !hpf$*        hx,hy,hz,dn) begin
   343         z=0.5*hz*(2*k-nz2-1)
   344         do 101 j=1,ny2
   345         y=0.5*hy*(2*j-ny2-1)
   346         do 101 i=1,nx2

   356    101  continue
   357  !hpf$  end on
   358     10  continue
   359         return
   360         end


表１．VPP FortranとHPF/JAで書かれた太陽風地球磁気圏相互作用の３次元MHD
コードによるFujitsu VPP5000/56での計算速度の比較

Table 1.  Comparison of computer processing capability between VPP Fortran
and HPF/JA in a 3-dimensional global MHD code of the solar wind-magnetosphere
interaction by using Fujitsu VPP5000/56.
 Number    Number of     VPP Fortran               HPF/JA
 of PE     grids        cpu time  Gflops  Gf/PE   cpu time  Gflops  Gf/PE

  1PE     200x100x478    119.607 (  0.17) 0.17 (scalar)
  1PE     200x100x478      2.967 (  6.88) 6.88       3.002 (  6.80) 6.80
  2PE     200x100x478      1.458 ( 14.01) 7.00       1.535 ( 13.30) 6.65
  4PE     200x100x478      0.721 ( 28.32) 7.08       0.761 ( 26.85) 6.71
  8PE     200x100x478      0.365 ( 55.89) 6.99       0.386 ( 52.92) 6.62
 16PE     200x100x478      0.205 ( 99.38) 6.21       0.219 ( 93.39) 5.84
 32PE     200x100x478      0.107 (191.23) 5.98       0.110 (186.13) 5.82
 48PE     200x100x478      0.069 (297.96) 6.21       0.074 (276.96) 5.77
 56PE     200x100x478      0.064 (319.53) 5.71       0.068 (299.27) 5.34

  1PE     500x100x200      2.691 (  7.94) 7.94       2.691 (  7.94) 7.94
  2PE     500x100x200      1.381 ( 15.47) 7.73       1.390 ( 15.37) 7.68
  4PE     500x100x200      0.715 ( 29.97) 7.47       0.712 ( 29.99) 7.50
  8PE     500x100x200      0.398 ( 53.65) 6.71       0.393 ( 54.38) 6.80
 16PE     500x100x200      0.210 (101.87) 6.37       0.202 (105.74) 6.61
 32PE     500x100x200      0.131 (163.55) 5.11       0.120 (175.50) 5.48
 48PE     500x100x200      0.100 (214.48) 4.46       0.091 (231.69) 4.82
 56PE     500x100x200      0.089 (239.48) 4.28       0.086 (244.85) 4.37

  4PE     800x200x670      7.618 ( 30.06) 7.52       8.001 ( 28.62) 7.16
  8PE     800x200x670      3.794 ( 60.36) 7.54       3.962 ( 57.81) 7.23
 12PE     800x200x670      2.806 ( 81.61) 6.80       3.005 ( 76.21) 6.35
 16PE     800x200x670      1.924 (119.00) 7.44       2.012 (113.85) 7.12
 24PE     800x200x670      1.308 (175.10) 7.30       1.360 (168.44) 7.02
 32PE     800x200x670      0.979 (233.85) 7.31       1.032 (221.88) 6.93
 48PE     800x200x670      0.682 (335.62) 6.99       0.721 (317.80) 6.62
 56PE     800x200x670      0.595 (384.61) 6.87       0.628 (364.87) 6.52

 16PE   1000x500x1118      9.668 (123.52) 7.72       9.619 (125.50) 7.84
 32PE   1000x500x1118      5.044 (236.73) 7.40       4.992 (241.83) 7.56
 48PE   1000x500x1118      3.550 (336.40) 7.01       3.479 (346.97) 7.23
 56PE   1000x500x1118      2.985 (400.04) 7.14       2.935 (411.36) 7.35
 32PE  1000x1000x1118      9.979 (239.33) 7.48       9.813 (243.37) 7.61
 48PE  1000x1000x1118      7.177 (332.79) 6.93       7.028 (339.85) 7.08
 56PE  1000x1000x1118      5.817 (410.55) 7.33       5.794 (412.23) 7.36
------------------------------------------------------------------------
 : Mflops is an estimated value in comparison with the computation by
   1 processor of CRAY Y-MP C90.



図説


第１図　惑星間磁場が北向きでかつ夕向きの場合の定常状態の地球磁気圏の３次元磁力線構造。
磁気圏境界の磁気リコネクションが北半球では夕側、南半球では朝側の高緯度尾部領域で起こ
るために朝夕非対称の磁力線構造になる。

Figure 1  Structure of 3-dimensional magnetic field lines in the steady earth's
magnetosphere for northward and duskward interplanetary magnetic field, where
asymmetry appears in the magnetosphere due to occurrence of magnetic reconnection
on the dusk side in the northern hemisphere and on the dawn side in the
southern hemisphere near the magnetopause.

第２図　太陽と地球の間で太陽風と惑星間磁場を観測しているACE衛星の１分毎のデータを
計算のインプットして用いた場合に、３次元グローバルMHDシミュレーションから得られた
地球磁気圏のスナップショットの例。

Figure 2 An example of snapshots of the earth's magnetosphere obtained by
3-dimensional global MHD simulation, when the observational data of the solar wind
and interplanetary magnetic field by ACE satellite was used as input of simulation.

第３図　総格子点数でみた２次元と３次元MHDシミュレーションの年毎の発展、総格子点数は
３年で４倍あるいは５年で１０倍に増加している。

Figure 3  Development of the 2 and 3-dimensional MHD simulations in view of total 
grid number used, the total grid number has increased 4 times every 3 years or 10 
times every 5 years.

第４図　MHDシミュレーションの総格子点数に対する１回の時間ステップを進めるために必要
とする計算時間。
　

Figure 4  Comparison of computation times corresponding to one time step advance 
versus total grid number in the MHD simulation codes.

第５図　MHDシミュレーションの総格子点数に対する必要な計算機メモリの総容量。

Figure 5  Comparison of total capacity of the necessary computer memory versus 
total grid number in the MHD simulation codes.

第６図　MHDシミュレーションの総格子点数に対するコンピュータの絶対処理能力の比較。

Figure 6  Comparison of the computer processing capability versus total grid 
number in the MHD simulation codes.










 Simulation Home Page.




  GEDAS Home Page
 


  STELab Home Page.