ＦＡＣＯＭ ＶＰＰ３００、５００での並列化 [Top][VPP2]
- 目次
- 導入
- ＶＰＰ上での実際の並列化作業
- ＶＰＰ上での具体的な作業（ｋ点に関する並列化）
- ＶＰＰシリーズの特徴
- ｋ点並列化
- グローバル配列の扱い方
- 重複ローカル配列、分割ローカル配列
- 具体的なグローバル配列
- 分割ローカル配列の運用
- 手続き分割
- ｋ点に関する並列化のまとめ
- バンドに関する並列化
- グラムシュミットの直交化ループの並列化
- バンドの並べ替え
- バンドに関わる部分以外の並列化
- バンドに関する並列化のまとめ
ここでは、筆者のＶＰＰ上での並列化作業の内容を比較的詳しく述べていき
たいと思います。
並列マシンは以外と昔からあり、筆者も１０年以上前にアライアントと言う
並列可能なマシンにちょっとだけ触った（？）ことがあります。ただ、本格的
に、並列化をしなければならないと思いはじめたのはここ１、２年のことです。
ここ１、２年（２、３年か？）で並列計算機に関しての状況も大きく変わりつ
つあります。ＶＰＰやＳＸ４のような本格的なベクトル並列型スーパーコンピュー
ターの出現（ベクトル並列型という意味ではそれ以前にもマシンは存在はして
いました。例えば、TITAN3000などがそうですが、あまり成功したとは言えま
せん。少なくともTITAN3000を製造したメーカーは潰れています。）、並列計
算機の世界の寵児ともいえるシンキングマシン社の倒産などが起こりました。
並列計算機（ソフトウェア環境も含めて）はいまだエンドユーザーにとって
扱い易いものとは言えません。いわゆる自動並列化コンパイラーは存在はしま
すが、いまだ玩具の領域を出ていません（つまり一般的には使いものにならな
い。当然、非常に特化した〔自動並列向きな〕状況では性能が出ることもあり
ます。）。
現在、並列計算機を使って性能を引き出そうとするなら、自分で自身の所有
するプログラムを並列化するしかありません（外注で頼むという手があるには
ありますが、それにはお金がかかり、期待するほどの性能が出せる確たる保証
もありません）。しかしその並列化も並大抵のこではできません。まず並列化
の仕方にいろいろあります。メーカー独自の仕様からＰＶＭ、ＭＰＩなどいく
つかの方法があります。おそらくメーカー独自仕様の並列化方法は今後廃れる
可能性が高いです。ＰＶＭ、ＭＰＩは今後の主流になるだろうとは思いますが、
これはライブラリ形式で並列化を実現するもので、並列化のための命令はサブ
ルーチン（FORTRANの使用を想定）を呼び出して行なわなければなりません。
この場合、（ＰＶＭ、ＭＰＩをインプリメントしていない）スカラーマシンで
の計算は不可能になります。
一方、ライブラリ形式ではなく、注釈行（並列化指示行）
の形式でプログラムの並列化を実現するものもあります。ここで扱うＶＰ
Ｐシリーズも注釈行形式での並列化が可能です（ＰＶＭ、ＭＰＩの実行もでき
ます）。この形式の利点は、他のシステムやマシンのコンパイラーでは注釈行
として解釈されるので、並列化したプログラムを別のスカラーマシンや（単一
ＣＰＵ）ベクトルマシン上で動かすことが可能です（並列化のために、プログ
ラム本体に大きな変更を加えている場合、速度が元より遅くなったり、最悪動
作しなくなる可能性もあります。）。
この他にもＨＰＦ（High
Performance FORTRAN）もあります。これは、言語仕様に並列化に関する要素
も取り込んだもので、並列化を意識してプログラムを書けば、ＨＰＦに対応し
ている並列マシン上でコンパイルすれば、そのままで並列動作が可能（なはず）
です。しかし、いまだＨＰＦに対応したシステムはほとんどなく（今年いくつ
かの並列マシン、システムが対応予定）、日進月歩のこのコンピューターの世
界でＨＰＦが必ず将来の主流になるという保証はどこにもありません。
いよいよ本題に入っていきます。筆者とベクトル並列型スーパーコンピュー
ターの代表とも言えるＶＰＰシリーズとの最初の出会いは、自分が通産省産業
融合領域研究所（融合研）の寺倉先生（博士課程のときの私の親分）のところ
と併任していた時（１９９４年頃）です。この時、融合研にはＶＰＰ５００が
導入されました。筆者もこの時ＶＰＰに触って（？）みました。この時は自分
のプログラムをＶＰＰ用に並列化することはあまり真剣には考えてはいません
でした。もっぱら単独のＣＰＵでの実行を行なっていました。それでも並列化
に関する講習会に参加したり、簡易マニュアル等を手に入れたりしてはいまし
た。
その後、１９９５年には古巣の物性研にもＶＰＰ５００が入ることとなり、
物性研は共同利用研で国研の職員も物性研のＶＰＰ５００を利用可能だったの
で、利用申請を行ない。１９９５年の１月（か２月）から年度末にかけての試
用期間からＶＰＰを利用していました。
この段階以降、このＶＰＰシリーズはいろいろな機関（高エネルギー研究所、
原研、理研などなど）に導入されました。またＶＰＰ以外にも、同じベクトル
並列型スパコンであるＳＸ４やスカラー並列のＩＢＭのＳＰ２やクレイのＴ９
０、Ｔ３Ｄ、Ｔ３Ｅ等、日立のＳＲ２２０１など多数の並列マシンが続々と登
場し、いろいろな機関に導入されました。
これをみて筆者は、並列化をいま始めないと時代に乗り遅れると思い、自身
の持つ第一原理分子動力学プログラムの並列化作業に本格的に乗り出しました。
既に融合研でＶＰＰには触っていたことと、物性研で開かれたＶＰＰ利用者講
習会（大変役に立ちました）で配られた資料をもとに構想半年、実労１カ月で
第一原理分子動力学プログラムの並列化ができました。
バンド計算（第一原理電子状態計算）プログラムは比較的並列化し易い部類
に入ります。特に、ブリュアンゾーンのサンプリングｋ点に関するループは、
非常に独立性が高く、大体の場合一番外側のループなので、並列化し易く、並
列化による効果も非常に高いです。ｋ点に関しての並列化での最大の問題は、
扱う系のサイズが大きくなると必要なｋ点の数が相対的に減ってしまうことで
す。これは、システムのサイズが大きくなるから、プログラムを並列化するこ
とによって高速化し、大規模な問題に対応することと相反してしまいます（ｋ
点数が系の大規模化で少なくなれば、並列化ができなくなる）。
但し、実際にはｋ点に関して並列化することは思ったほど非効率ではありま
せん。本当の意味で（少なくとも系が金属的なバンドの場合）ｋ点がΓ点（つ
まり１点だけ）で済むのは、系のユニットセル内の原子数が１０００以上の場
合です。実際は１００個程度の系で、Γ点だけの計算が多いですが、その系が
金属的な場合、本当のことを言えばΓ点だけというのは問題があると思います。
現実には計算資源と時間の問題を考えれば止むおえないのですが、ｋ点での並
列化ができれば、この問題はある程度解決します。
今後、擬ポテンシャル＋平面波の計算でも、遷移金属とその化合物や磁性の
問題を扱う場合もどんどん出てくるだろうと思います。このようになるとます
ます複数ｋ点を使用する必要が出てくると思われます。このようなことからも
ｋ点に関しての並列化にはいまだ十分な意味と価値があります。加えて、ｋ点
に関する並列化は後述するバンドや、逆格子点に関する並列化よりずっと少な
い作業量で、目的が達成できるので、並列化をまず最初に試してみるには丁度
良い課題ではないかと思います。
さてここからが本当の本題です。まずＶＰＰでの並列化の仕方について話し
ていきたいと思います。
気をとりなおして続けたいと思います。
尚、ここで私が話すことは、全て自身の（乏しい）経験から成り立っているの
で、舌足らずなところや、間違った記述があるかもしれません。その場合、遠
慮無く御指摘頂けるとありがたいです。（メイルアドレス： kobayashi.kazuaki-@-nims.go.jp ← "-@-"
は変なメイル対策）
以下に、ＶＰＰシリーズの（個人的な独断と偏見による）特
徴を挙げます。
- ベクトル並列型スーパーコンピューター（単一ＣＰＵでも速い）
- 分散メモリー型 （ＳＸ４〔日本電気〕は共有メモリー型）
- 並列化指示行（ディレクティブ）での並列化が可能
- ＳＥ、サポートが比較的良い
- 少なくとも国内ではこの手のマシンとしては良く売れて（普及して）いる
共同利用研である物性研にＶＰＰ５００があるので、バンド屋さんがこの機
械を触る機会は多いかと思います。ＶＰＰは分散メモリー型のベクトル並列型
スーパーコンピューターで、１ＣＰＵだけでも１．６ＧＦＬＯＰＳの速度があ
ります。従って、単独のＣＰＵのみでも扱う系によっては十分計算が可能です。
ただ、分散メモリー（各ＣＰＵ毎にメモリーが分散していて、他のＣＰＵに跨っ
た〔その跨われたＣＰＵは使用しないで〕形でのメモリーの利用はできません。）
なので、単独のＣＰＵのみでの利用では、そのＣＰＵに載っているメモリー以
上の大規模系の計算は不可能です。そのような場合には諦めるか、より速くて
メモリーを沢山積んだスパコン等を利用するか、あるいは並列化を目指すかの
いずれかです。そしてここでは並列化について話しているわけです。
ＶＰＰの並列化での一番の特徴（と私は思っています）は先のリストにも挙
げた、並列化指示行（ディレクティブ）による並列化が可能なことです。ＶＰ
Ｐでは!XOCLで始まる行が、並列化指示行になっています。
例 !XOCL PARALLEL REGION
!XOCL以降にあるのが、並列化に関しての具体的な命令です。先に挙げたＶ
ＰＰに関しての有用なリンク先に、次のような表（←既にアクセス不能）があ
ります。これに、並列化指示行の命令（ディレクティブ）の一覧があります。
この方法の大きな利点として、並列化指示行がＶＰＰ以外のマシンでは、コ
メント（注釈行）として扱われることです。従って、ｋ点に関する並列化のよ
うに、もとのバンド計算プログラムに並列化指示行を加えただけ（厳密にはそ
うではありません）だと、容易に他の非並列型マシン（単独ＣＰＵのスパコン、
ワークステーション、ＰＣ等）で、プログラムをそのままの形で実行すること
が可能です。これはバージョンの管理や、並列化に伴うデバッグを考える上で
非常に役立ちます。
バンド計算において、ｋ点のループに関しての並列化で、この一覧（←既に
アクセス不能）にある命令の内で、使用する可能性が高いものを以下に挙げま
す。これらの命令には必ず!XOCLが最初に付きます。
- PROCESSOR
- SUBPROCESSOR
- INDEX PARTITION
- PARALLEL REGION
- END PARALLEL
- GLOBAL
- LOCAL
- SPREAD DO
- END SPREAD
一覧には沢山命令がありますが、ｋ点での並列化で使う命令は、せいぜい上
に挙げた程度です。（プログラムの並列性能を上げるためにチューンする場合
には上記以外の命令を使う可能性があります。）ＶＰＰ（富士通）では並列化
は、手続きの分割と、データの分割の２つに大別することができるとしていま
す。手続きの分割は、つまりＤＯループを分割することです（当然ＤＯループ
以外の手続きの分割もある）。データの分割はデータつまりFORTRANの世界で
言うところの配列の分割です。分割によって分けた仕事（手続き）やデータ
（配列）が各ＣＰＵ（ＶＰＰの世界ではＰＥと言っています）で計算、実行さ
れます。
手続きの分割は、大抵ＤＯループの分割で、特にｋ点並列の場合、ＰＥ間
（ＣＰＵのこと）でのデータの通信の必要もほとんどなく、比較的楽に並列化
作業が行なえます。
ｋ点（ブリュアンゾーンのサンプリング点）のループは、非常に独立性が高
く、非常に並列化が難しくなる２重以上のループもありません。また、あるＰ
Ｅ上の計算が、他のＰＥ上のｋ点毎のデータを必要とすることはほとんどあり
ません。ただ、電荷密度や全エネルギーなどは、バンドとｋ点の足し上げで得
られる量なので、このような物理量を求める時には他のＰＥのデータを参照す
る必要がありますが、それは非常に簡単にできます。（むしろ、そのような参
照のための操作が必要かどうかを見極める方が難しいです。）
手続き分割のための命令は、先のリストの最後の２つSPREAD DOとEND
SPREADで行ないます。この２つの間にＤＯループ（手続き）が入ります（詳し
いことは後述します）。ｋ点並列ではこれだけで十分だと思います。
一方、データの分割はちょっと難しいと言えます。一般的に並列マシンにお
いて、手続きの分割よりデータの分割の方がずっと難しい（と私は思っていま
す）です。ＶＰＰでは更にＶＰＰ固有の難しさがあると思います。
ＶＰＰでは並列用にグローバル配列と分割ローカル配列があります。並列化
を考慮していない普通（従来通り）の配列は重複ローカル並列と言います（と
思っています）。ここで問題なのが、グローバル配列と分割ローカル配列の意
味と関係です。グローバル配列は、ＰＥ間に跨った広域的な配列です。跨って
はいますが、この配列は事前に設定した並列計算のためのＰＥ数分に分割され、
それぞれ分割されたものが各ＰＥのメモリー上に存在します。並列計算を行なっ
た時に、手続きの分割を行なっていなくてもグローバル配列に値を代入したり、
演算を施したりすると、並列計算している全ＰＥ上のグローバル配列に対して、
代入、演算が行なわれます。
つまりＡ（１０００）がグローバル配列だったとして、８ＰＥでの並列計算
が行なわれている場合（ＩＰ＝８、これは分割数）
!XOCL SPREAD DO /IP
DO 1000 I=1,1000
A(I)=0.0D0
1000 CONTINUE
!XOCL END SPREAD
でも
DO 1000 I=1,1000
A(I)=0.0D0
1000 CONTINUE
でも同じ計算が行なわれます。前者の例ではSPREAD DO文を使って、ＤＯルー
プの手続き分割が行なわれ、８つのＰＥで分割（この場合等分割を想定）して
実行が行なわれています。後者のルーチン例では、並列実行は行なわれず、８
個あるＰＥのどれか１つしか動きません。Ａ（１０００）は１から１２５、１
２６から２５０、２５１から３７５、３７６から５００、５０１から６２５、
６２６から７５０、７５１から８７５、８７６から１０００までと各ＰＥのメ
モリー上に、そのデータ（値）が割り振られています。並列実行している場合
は、各１２５個ずつの配列値が、各ＰＥ上で計算されますが、並列実行しない
時はどれか１つのＰＥ（ＰＥの選択はシステムが勝手に決めます）で１から１
０００までの計算が行なわれます。ＶＰＰは分散メモリーなので、直接他のＰ
Ｅ上のメモリーにアクセスできないので、これでは正しく計算が行なわれない
ようにも思われますが、グローバル配列は、ＰＥに跨って（自動的に）メモリー
にアクセスできる（厳密に言うと、ＰＥ間で通信しているのだと思います）よ
うにしてある、特殊なものになっています。
従って、どちらの場合でも、グローバル配列Ａに対するゼロ初期化は行なわ
れます。
ただし、グローバル配列はＰＥに跨った操作が可能ですが、その分操作に時
間が非常にかかります（ＰＥ間で通信しているから）。並列動作させてもとて
つもなく遅くなります。従ってグローバル配列を直接並列演算のためには使用
しません。高速な並列計算をするためには、分割ローカル配列を利用します。
分割ローカル配列は、各ＰＥ毎に分割された、局所的な配列です（つまり、自
分の担当する以外の部分はなにも起こりません）。これは個々のＰＥ上で分割
された部分の計算しかしませんから、先ほど述べたグローバル配列での速度の
遅さの問題はありませんが、このままでは他のＰＥ上のデータを扱うことも、
参照することもできません。例えば、和をとるループで、各ＰＥ毎に和の計算
を分担させると、各ＰＥ毎に分割されて計算された和が出てきます。これを全
てＰＥに関して足し上げて最終的な答が出てきます。このＰＥに関して、各Ｐ
Ｅ毎の部分的な和の結果を足すことは、分割ローカル配列だけでの操作では不
可能です。（ここらへんや、以下の知識はかなり生
半可なので、記述や、概念その他何か間違いがあれば御指摘頂くとありがたい
です。）
この問題のためにグローバル配列があります。グローバル配列は使用するＰ
Ｅ全てに跨っていて、かつあるＰＥから、他のＰＥをみることができます。た
だしグローバル配列は直接扱うと大変時間を浪費するので、そのまま使うわけ
にはいきません。そこで、グローバル配列と分割ローカル配列をイクウィバレ
ンス文（EQUIVALENCE文）でメモリーを共有するようにします。こうしておく
と各ＰＥ上で分割ローカル配列の内容が変われば、自動的にグローバル配列の
値も変わります（メモリーが共有されているから）。
ここでグローバル配列が定義の例として、最初のリスト項目にある”やさし
い（？）バンド計算プログラムの作り方”で説明した、インクルードファイル”
PACVPP”を示したいと思います。このファ
イルの最後の方で、グローバル配列が定義されています。（最初に、！XOCLが
付いているのですぐにわかると思います。）
以下に、その部分を抜粋して示します。
PARAMETER(IPARA=8)
!XOCL PROCESSOR PQ(IPARA)
DIMENSION ZFBB2(KEG,KNV3,KATM,6)
!XOCL INDEX PARTITION IG=(PQ,INDEX=1:KNV3,PART=BAND)
!XOCL INDEX PARTITION IH=(PQ,INDEX=1:KO,PART=BAND)
!XOCL GLOBAL SNL(:,/IG,:,:),SNL2(:,/IG,:,:),SNL3(:,/IG,:)
!XOCL GLOBAL ZAJ(:,:,/IG),ZFBB(:,/IG,:,:)
!XOCL GLOBAL EKO(:,/IG),RAK(:,/IG),OCCUP(:,/IG)
!XOCL GLOBAL ZFBB2(:,/IG,:,:)
C
!XOCL GLOBAL NGPT(:,/IH),NAPT(:,/IH),OP(:,:,/IH),TAU(:,/IH)
まず最初の行のパラメーター文で示してある、IPARAで使用するＰＥ数を指
定しています。次のPROCESSOR文で、実際にプログラム全体が使用するＰＥ数
を設定します。そして、INDEX PARTITION文で、ＰＥ数による分割の仕方を設
定します。その仕方の指標となるのがIGとIHです。IGはｋ点に関しての分割、
IHは系の対称性に関してのオペレーター数の分割に関わっています。
ここで、PART=BANDは等分割を意味します（そうでない分割の仕方もある）。
そして最後にグローバル配列を定義しています。これらの配列は既にPACVPPで
（コモン配列変数として）宣言されています。この宣言の時、コモン文内でグ
ローバル配列以外の配列を同時に定義できません。
正しい例
COMMON/PSSNL/ SNL(KNG1,KNV3,KTYP,10)
& ,SNL2(KNG1,KNV3,KTYP,9)
& ,SNL3(KNG1,KNV3,KTYP,4)
悪い例（余計な配列AC、BCがある）
COMMON/PSSNL/ SNL(KNG1,KNV3,KTYP,10),AC(KTYP,2),BC(KTYP,2)
& ,SNL2(KNG1,KNV3,KTYP,9)
& ,SNL3(KNG1,KNV3,KTYP,4)
この!XOCL GLOBALでは配列の定義をしているのではなく、具体的な分割を行
なっています。配列SNLを例にとると、SNLはSNL(KNG1,KNV3,KTYP,10)と定義さ
れています。このKNV3がｋ点数です。従って、このKNV3の部分が分割されます。
!XOCL GLOBAL SNL(:,/IG,:,:)
SNLは既に定義されているので、分割しない部分は：で表現されます。/IGは
ｋ点数KNV3に関して、８等分割することを示しています。因みに一度決めた分
割位置（分割軸）を途中で変えることは（他に別の分割位置で分割しておいた
ワーク用配列を使わない限り）できません。
ここで問題となるのは、ｋ点で並列化するとして、では一体どのデータ（配
列）を分割すれば良いかです。前にも述べましたが、データの種類として並列
化に関わるグローバル配列と分割ローカル配列、そして並列化と関係ない普通
の配列（重複ローカル配列）があります。ここで、重複ローカル配列は従来か
らある配列で、並列化に関して何も考慮していないものを指しています。この
配列はＶＰＰ上で並列計算を行なう時に問題があります。それは、この重複ロー
カル配列は各ＰＥ上でそれぞれ独立に定義されてしまいます。つまり８ＰＥ並
列計算を行なおうとして、Ａ（１０００）という重複ローカル配列を定義する
と、各ＰＥ毎にＡ（１０００）が定義されます。因みに、グローバル配列〔Ｇ
（１０００）とする〕、分割ローカル配列〔Ｂ（１０００）とする〕では８Ｐ
Ｅ並列（等分割として）なので、各ＰＥには１２５個分のデータしかメモリー
上に保持されません。これに対して、重複ローカルでは各ＰＥにまるまる１０
００個分のデータがメモリー上に保持されます。
ＶＰＰは分散メモリー型で、各ＰＥの保有するメモリーにが限りがあります
（ＶＰＰ３００、５００、７００で１ＰＥが持てる最大メモリー量は異なりま
すが、多くても２ＧＢまでです。実際には予算その他の理由から２５６ＭＢ、
５１２ＭＢ程度です。）。従って、重複ローカル配列が沢山あると、個々の配
列の大きさはそれほどでもなくても、全体としてはＰＥ上の限られたメモリー
を圧迫してしまいます。かといって、ｋ点に関して並列化する場合、ｋ点と何
ら関係のない配列を、安易に分割（グローバル配列や分割ローカル配列に）す
るわけにはいきません（場合によっては、分割できる場合もある）。
何がグローバル配列になり、何を分割ローカル配列にする（実は、並列化さ
れていない元プログラムには、分割ローカルに対応する配列〔データ〕はあり
ません。分割ローカル配列は大抵、グローバル配列と１対１に対応するように
新たに作ります。〔前述のようにEQUIVALENCE文でメモリーを共有する〕）か、
そしてどれが重複ローカル配列として残るかを考えなければなりません。
ｋ点並列の場合、分割の対象となるのはｋ点に関わる
ループと配列です。バンド計算においてｋ点の計算が関わっている部分は以下
のようになります。
- 非局所擬ポテンシャルの実空間から逆空間への変換部分
- 波動関数更新部分（対角化、ＳＤ、ＣＧ法など）
- フェルミ面決定部分
- 原子間に働く力の計算部分
- 電荷密度の計算部分
- 全エネルギー計算部分
- ストレス計算部分
他にも、こまごまとしたｋ点に関わるループ（ルーチン）がありますが、主
なものは以上に挙げた７つの部分です。並列化に関しては、この７つの部分に
対して行なえば、取り敢えずは十分かと思います。
ｋ点に関してのループは非常に大きいので、その中に、ＦＦＴや波動関数の
直交化操作の部分はすっぽりと入ってしまいます。またＦＦＴや直交化はｋ点
に対して完全に独立しているので、下の図にあるようにこれらの個々のルーチ
ンに対して並列化を考える必要がありません。
---- ｋ点のループ ------ <--- ここが並列化される
|
| -- バンドのループ ----
| |
| ------> ＦＦＴルーチン
| |
| ------> 直交化ルーチン
| |
| ----------------------
------------------------
次に、具体的なグローバル配列を示します。
!XOCL GLOBAL SNL(:,/IG,:,:),SNL2(:,/IG,:,:),SNL3(:,/IG,:)
!XOCL GLOBAL ZAJ(:,:,/IG),ZFBB(:,/IG,:,:)
!XOCL GLOBAL EKO(:,/IG),RAK(:,/IG),OCCUP(:,/IG)
!XOCL GLOBAL ZFBB2(:,/IG,:,:)
!XOCL GLOBAL NGPT(:,/IH),NAPT(:,/IH),OP(:,:,/IH),TAU(:,/IH)
（１）SNL、SNL2、SNL3は非局所擬ポテンシャルのグローバル配列です。スト
レスの計算を考えない場合（これが普通）はSNLだけです。
（２）ZAJは波動関数（固有ベクトル）のグローバル配列、ZFBB及びZFBB2は作
業用（力やストレスの計算で利用）のグローバル配列です。
（３）EKO（固有値）、OCCUP（バンド占有率）はフェルミ面決定の時に使用す
るグローバル配列です。
（４）RAKは１／（ｋ＋Ｇ）を格納したグローバル配列です。
（５）NGPT、NAPT、OP、TAUは、実はｋ点並列用のものではなく、系の対称性
に関しての並列化のためのグローバル配列です。
ｋ点並列で使われる主要な配列は、以上に挙げたもので大体尽くされている
と思います。当然、他にもｋ点に関わる配列データは存在しますが、PACVPPを
見てもらえばわかると思いますが、どれも小さい配列で数も多くないです。従っ
てそのようなこまごましたものは並列化の対象外としました。
（５）に関して、系の対称性の部分の計算では、ＢＣＣ、ＦＣＣのように対
称性の高い場合、非常に計算時間を要することがあります。そのため、系の対
称性に関しての並列化も行なっています。これは前述の重複ローカル配列に関
する対策にもなっていますが、この部分の配列は、あまり大きくなく効果はほ
とんど期待できません。
ここで全エネルギーの計算をしているルーチン（ｋ点に関して並列計算して
いる）を示します。
このサブルーチンでは、波動関数更新をしているルーチン（別のサブルーチ
ン）で求めた、固有値の値と更新される前と後の電荷密度とを使って、系の全
エネルギーを求めています（詳しくはバンド計算に関しての論文、参考書、自分のＤ論等参照）。まずは
ここで定義されている分割ローカル配列を見てみましょう。
!XOCL SUBPROCESSOR PS(IPARA)=PQ(1:IPARA)
!XOCL INDEX PARTITION IP=(PS,INDEX=1:KNV3,PART=BAND)
DIMENSION EKK(KEG,KNV3),OCCUU(KEG,KNV3)
!XOCL LOCAL EKK(:,/IP),OCCUU(:,/IP)
EQUIVALENCE (EKO,EKK),(OCCUP,OCCUU)
ここで必要なグローバル配列はEKOとOCCUPのふたつです。それに対応する分
割ローカル配列、EKKとOCCUUがここで定義（宣言）されています。グローバル
配列はコモンとしてインクルードファイルPACVPP内で定義されていますが、分
割ローカル配列は、あくまでローカルな配列としてここで定義されています。
EKK、OCCUUは、まずDIMENSION文でローカルな配列として定義された後、!XOCL
LOCAL文で分割ローカル配列としての分割が行なわれます。
並列化で使うＰＥ数をサブルーチン毎に定義しておく必要があり、そのため
の命令が上記のSUBPROCESSOR文です。それに対応するINDEX PARTITION文が下
の行にあります。使い方はPACVPP内でのPROCESSOR文、INDEX PARTITION文と同
じです。SUBPROCESSOR文ではPROCESSOR文で定義しておいたＰＥ数を越えない
範囲で自由なＰＥ数の設定ができますが、普通は同じ数で使用するのが良いか
と思います（ｋ点並列でいちいち変える必要は全くないです）。
話しは戻って、最後にグローバル配列と分割ローカル配列は、EQUIVALENCE
文を使ってメモリー空間が共有されます。
EQUIVALENCE (EKO,EKK),(OCCUP,OCCUU)
EKK、OCCUUは、このENERGY（全エネルギーを計算するサブルーチン）でロー
カルに定義されているだけです。他のサブルーチンでも固有値、占有数の分割
ローカル配列として、同じ名前で使われていますが、別の名前で定義しても一
向にかまいません（そうする必要は全くありませんが）。メモリー的にもグロー
バル配列とEQUIVALENCE文でメモリーを共有しているので、分割ローカル配列
定義による（メモリーの）損はありません。
では次に、サブルーチンENERGY上で行なわれている手続き分割について説明
して行きたいと思います。以下に、手続きで分割している部分を示します。
FFF = FLOAT(KV3)
TTT = 0.D0
!XOCL SPREAD DO /IP
DO 2 I=1,KV3
DO 100 IBAN=NBD1,NBD2
TTT = TTT + OCCUU(IBAN,I)*
& EKK(IBAN,I)
100 CONTINUE
2 CONTINUE
!XOCL END SPREAD SUM(TTT)
ETOTAL = ETOTAL + 2.D0*TTT/FFF
手続きの分割にはSPREAD DO文を使います。他にも、手続き分割の命令、方
法はありますが、ｋ点並列ではSPREAD DO文だけで十分です（と思っている）。
ここで重要なのは、必ず、SPREAD DO文は分割すべきDOループ（ここではｋ点
に関するDOループ）の直前に置くことです（注釈行などが入っていてもよいか
は試していないので不明）。
ここで一番重要なのは、SPREAD DO文に対応するEND SPREAD文です。普通、
各ＰＥ上で並列計算しているだけでは、END SPREADだけで良いのですが、この
全エネルギーの計算はｋ点に関しての和をとっています。和をとるので、各Ｐ
Ｅ毎にある部分和から全体としての和を求めなければなりません。ＶＰＰでは
これは非常に簡単に行なうことができます。上の例にもあるように、ＤＯルー
プで足されていく変数TTT（これは単なる変数で、並列化に関しての特別な定
義は一切必要としません。この場合、各ＰＥ毎にTTTのｋ点に関しての部分和
が存在します。）をEND SPREAD文のSUMオプションで指定してやればTTTの全体
としての和が自動的に求められます（確かグローバルサムという）。
これでデータ分割と手続き分割の仕方の大体のところを説明しました。尚、
非局所擬ポテンシャル部分のサブルーチンが擬ポテンシャルデータベースの説明のと
ころで例示されています。擬ポテンシャル
データベースの説明のところでは非局所擬ポテンシャルデータの入力と操作の
仕方の説明が行なわれましたが、このプログラムルーチン例のサブルーチン
KBMATでは、ｋ点に関しての部分の並列化が行なわれていますので、参考にな
るかと思います。並列化の内容としては先の全エネルギーのルーチンでの並列
化と大差ありません。
ｋ点に関しての並列化は、並列化すべきデータ（配列）の選択と、それに付
随するグローバル配列と分割ローカル配列の割り振り、手続き分割のためのＤ
Ｏループの分割、並列化されたＤＯループ内での和の計算のためのグローバル
サムの考慮、などを考えれば、比較的簡単で少ない作業量で、既存のバンド計
算プログラム（並列化を行なう人は対象となるプログラムを熟知していること
が大前提です。）を並列化することが可能です。
尚、筆者の第一原理電子状態計算コードは、ｋ点に関する並列化によって、
８ＰＥ並列の場合１ＰＥでの計算速度と比べて、最も良い場合で７倍近い高速
化が達成されますが、これはあくまで最高な場合で、大体５から６倍程度の高
速化が得られています。もっと徹底的な並列化に関する最適化を行なえば、一
層の高速化が計れると思われます（が、実際に時間も資源もない）。
では次に、バンド（エネルギーレベル）に関する並列化について話して行き
たいと思います。ｋ点の並列化が比較的容易なのに対して、バンドに関しての
並列化は少々難しくなります。何故ならば、ｋ点がループとして完全に独立し
ており、（普通のバンド計算では）あるｋ点のデータが、他のｋ点のデータを
参照することは（全エネルギーのように最終結果としての和をグローバルサム
する場合を除いて）全くありません。
一方、バンド（エネルギーレベル又は、状態数ともいう）に関するループは、
ｋ点と同じくらい外側の大きなループではありますが、独立して扱えない部分
が存在します。独立していないループは直交化に関する部分で、直交化を行なっ
ているバンド計算では大抵はグラムシュミットの方法が用いられています。
（例えば、従来通りの対角化手法を用いれば、固有ベクトルの直交化は対角化
に関しての数値演算ライブラリーのルーチンが自動的に行なってしまいますが、
大抵システムに備わっている演算ルーチンを利用するので、並列化が行なわれ
ていなかったり、不十分です。また、別の方法としては直交化を行なわないで、
波動関数を逐次更新する方法もありま
す。）
グラムシュミットの直交化ループはバンドに関しての２
重のＤＯループになっているため、そのままではうまく並列化が行なえません。
この部分を並列化させるにはいろいろな方法が考えられますが、先ほど示した、
富士通のＶＰＰに関するＷＷＷ情報（←既にアクセス不能）で、山崎さんが”
アトムテクノロジー研究体におけるＶＰＰ５００の応用事例”で、グラムシュ
ミットの直交化部分の並列化について詳しく述べておられます。非常に参考に
なると思います。
ＶＰＰの使い方の講習会の時に、富士通が配った資料の中に、グラムシュミッ
トの直交化部分のコードを並列化した参考例があります。しかし、これは勝手
に一般に公開するわけにはいかないので、興味のある方は、何とかして自力で
手に入れて（探して）下さい。ここで使われている並列化手法は、BROADCAST文という並列化指示行を使って、全てのＰＥに
必要な固有ベクトルの情報を送信します。そしてバンドに関する２重ループの
問題を回避して、並列化をしています（これだけでは何もわからない^^;）。
ただ、この場合、他のＰＥへのデータの転送が伴うため、並列化効率としては
損をしていることになります。
また先の山崎さんの説明にある、バンド並列から展開
係数（平面波の数、あるいはＧmaxまでの逆格子点Ｇに関して）に並列化の軸
を変える方法もあります。ただし、グローバル配列は一度決めた並列化の軸を
途中で変更はできません（グローバル配列G(/IP,:)を途中でG(:,/IP)とはでき
ない）。そうするためには、もうひとつ別の並列化の分割軸を持ったグローバ
ル配列を用意しなければならなくなります。そして、異なる分割軸を持つ配列
同士で通信が必要になります（メモリーと通信量を消費する、これは山崎さん
も指摘しています）
。グローバル配列の並列化に関する分割軸を途中で、
自由に変えられる（通信を必要としない）うまい手があると、余計な配列を用
意する必要もなく大変話しが簡単になるのですが、富士通では何かアイデアは
ないのでしょうか？
バンドに関する並列において、グラムシュミットの直交化の部分は、先にも
述べた富士通の配布資料をもとに並列化することができました。実際の並列化
作業において問題となったのはバンドの並べ替えの問題でした。筆者の持って
いる従来の第一原理電子状態計算プログラムはバンド（エネルギー準位、エネ
ルギーレベル、バンドレベルともいう）を小さい順に並べるようになっていま
す。対角化を行なう場合は、汎用の科学技術計算ライブラリルーチンを使う場
合は、大抵バンドを昇順、降順に並べることが指示できますが、波動関数（固
有ベクトル）と固有値（バンド）を逐次計算（カー・パリネロ法が代表例）で
求める場合、最初バンドを小さい順に並べておいても、計算途中でその順番が
保持される保証はそのままでは全くありません。
一方、グラムシュミットの直交化ではバンドが小さい順に並んでいることが
前提条件です。従って、直交化する前にバンドを並べ替えて、小さい順に並ぶ
ようにしておかなければなりません。この並べ替えを行なう場合、バンドだけ
ではなく、当然そのバンドに対応する固有ベクトル部分のバンドに関しての並
びも替えておかなくてはなりません。以下に、その部分を示します。
DO 262 IBAN=NBD1,NBD2-1
DO 262 JBAN=IBAN+1,NBD2
IF(EKK(JBAN,NNN).LT.EKK(IBAN,NNN)) THEN
EE =EKK(IBAN,NNN)
EKK(IBAN,NNN)=EKK(JBAN,NNN)
EKK(JBAN,NNN)=EE
DO 270 IV=1,IIBA
ZTV = ZZZ(IV,IBAN,NNN)
ZZZ(IV,IBAN,NNN) = ZZZ(IV,JBAN,NNN)
ZZZ(IV,JBAN,NNN) = ZTV
270 CONTINUE
END IF
262 CONTINUE
これは、計算量としては思った以上に大きなものです。バンドに関する２重
のループの中に、更に固有ベクトルのスワップ（平面波数分）をするループが
あり、グラムシュミットに匹敵するほどの量になりす（実際は、イタレーショ
ン毎に並びが替わるバンド部分は、そう多くはないのでバンド数の２乗かける
平面波数分という計算量よりはずっと少ないですが）。並べ替えもバブルソー
トなので効率は非常に悪いです。そして、この部分をバンドに関して並列化す
るのも、非常に難しいです（メモリーと計算量を犠牲にしないで並列化は不可
能）。
これに関して、結局わからず前述の山崎さんに相談してみると、この問題は
すぐに解決することを教えて貰いました。バンドに関して、その並びが問題と
なるのは、先に挙げたグラムシュミットの部分だけで、それ以外の部分で小さ
い順に並んでいることが必須である部分は、（少なくとも筆者の使っているプ
ログラムでは）存在しないことがわかりました。更に、上のルーチンのように
固有ベクトルそのものを演算する必要は全くなく、固有ベクトルのバンドの指
標のみを並べ替えるマッピングをする配列を考え、それを並べ替えるだけで良
いことがわかりました。そのように修正したループを以下に示します。
DO 262 IBAN=NBD1,NBD2-1
DO 262 JBAN=IBAN+1,NBD2
IF(EKB(JBAN,NNN).LT.EKB(IBAN,NNN)) THEN
II =IMAP(IBAN,NNN)
IMAP(IBAN,NNN)=IMAP(JBAN,NNN)
IMAP(JBAN,NNN)=II
EE =EKB(IBAN,NNN)
EKB(IBAN,NNN)=EKB(JBAN,NNN)
EKB(JBAN,NNN)=EE
C DO 270 IV=1,IIBA
C ZTV = ZZZ(IV,IBAN,NNN)
C ZZZ(IV,IBAN,NNN) = ZZZ(IV,JBAN,NNN)
C ZZZ(IV,JBAN,NNN) = ZTV
C 270 CONTINUE
END IF
262 CONTINUE
IMAPなるバンドとｋ点に関しての配列を用意し、それがバンドに連動して並
べ替えるようにし、２７０のＤＯループは計算する必要がなくなりました。た
だし、この並べ替えのループそのものはバンドに関しての並列化はいまのとこ
ろ行なっていません。本当はもっと良いやり方があると思いますが、いまのと
ころバンドに関しての並列化はこの段階までです。
当然、最初から直交化を必要としないようなバンド計算では、上で述べた意
味でのバンドの並び替えの必要はありません（別の部分で必要になるかもしれ
ませんが、筆者はこの方法を使ったことがないので、これ以上は何とも言えま
せん）。もしバンドの順番に計算が全く依らなければ、バンドの並列化もｋ点
並列化と同じ位の作業で遂行することができます（できるはずです）。
バンドの並列化を行なう場合、ｋ点では並列化できた部分の中に、並列化が
できなくなってしまうものがあります。先にｋ点並列のところで挙げたリスト
は、
- 非局所擬ポテンシャルの実空間から逆空間への変換部分 <---- バンド並列不可
- 波動関数更新部分（対角化、ＳＤ、ＣＧ法など） <------ バンド並列可
- フェルミ面決定部分 <------ バンド並列可
- 原子間に働く力の計算部分 <------ バンド並列可
- 電荷密度の計算部分 <------ バンド並列可
- 全エネルギー計算部分 <------ バンド並列可
- ストレス計算部分 <------ バンド並列可
となります。
非局所擬ポテンシャルを逆空間にフーリエ変換する部分は、バンドに関して
のループが存在しません。また、それに関わる配列SNL1、SNL2、SNL3もバンド
に全く依らない配列変数です。この部分は、大雑把に平面波数×動径方向のｒ
の点数×ｋ点数の計算量が必要で、特に筆者の計算では、動径方向の計算点数
が多いため、思った以上に計算量が必要になります。
ここの計算量を減らす方法としては、動径方向の計算（数値積分）部分にお
いて、（Ａ）解析的に計算できるように、非局所擬ポテンシャルを何らかの解
析関数で表現し直す。（Ｂ）補間（スプライン補間など）を使って、積分すべ
き点数を減らす（点数を少なくすればルジャンドルの積分法が使えるかもしれ
ません）。（Ｃ）動径方向の積分ループに関して並列化してしまう。などが考
えられます。
筆者は、取り敢えず（Ｃ）の方法を考えました。これは単なる好みの問題で、
特に強い理由があるわけではありません。この部分は計算量は非常に大きくは
なりますが、電子状態のみの計算を考える場合は最初の１回のみ計算するだけ
だったので、動径方向の積分は真正直に計算していました。ただ系が大きくな
ると１回でも、大変な計算時間を消費するので、もう少し高速化してみようと
いうことで、どうせならいろいろな並列化を試してみようということで動径方
向の積分ループに関しての並列化を行なってみました。
この部分の並列化は他のバンド並列化部分と競合することはありません。ｋ
点並列では非局所擬ポテンシャル（逆格子空間Ｇ用）の配列変数SNL1、SNL2、
SNL3は、並列化対象とはならないので、グローバル配列として宣言する必要も、
ローカル配列も用意する必要はありません。むしろ問題だったのは、もともと
ｋ点並列化されていたプログラムを出発点として、バンド並列化を開始したの
でｋ点並列からバンド並列への移行、前述のｋ点での並列化は可能だったが、
バンド並列はできない部分の修正、その修正したルーチンを、今度は動径方向
で並列化するためのプログラムの変更にちょっと手間取りました。そして、動
径方向に関しての並列化は行なえるようになったのですが、いまのところあま
り効率良く計算（高速化）できていません。むしろこの部分を並列化しないで
バンド部分のみ並列化していたバージョンの方が速いくらいです。原因は現在
究明中ですが、あまり捗っていません。
バンドの並列化は、扱う系のサイズが大きくなるほど並列化すべきバンドの
数はどんどん増えていくので、ｋ点並列での問題（系を大きくすると相対的に
必要なｋ点数が少なくなる）は全くありません。従って、大きな系に対しても
並列計算のためのＰＥ（ＣＰＵ）数を増やしていくことで対応することができ
ます。バンドのループはｋ点のループに次ぐ大きなループで、直交化の部分以
外ではバンド部分の計算の独立性も高いので、徹底的なバンドに関しての並列
化を遂行すれば、非常に効率的な並列計算が実現できます。
ただし、バンドの並列化に関しては、かなり苦労した部分があり（先の、バ
ンドに関する並べ替え）、筆者もｋ点並列が構想半年、実働１カ月でほぼ完成
したのに対し、バンド並列は構想１カ月、実働半年かかりました。実働のほと
んどの部分がデバッグです。はずかしながらバンドの並べ替えの部分は、何度
も間違いを繰り返し、正しい結果を得るまでに大変な労力と時間を費やしまし
た。
（最後に残った問題）が、先のバンド並べ替え
のループでIMAPのみ並べ替え、エネルギー固有値EKBの並べ替えを忘れていた
のに気付かず、何度も再計算を繰り返しても正しい結果が出なかったことです
（筆者はこれに大変悩みました）。
結局、ある日に何となくソースコードを見ているうちに、ふとIMAPを変更す
れば、それに付随して固有値も並べ替えておかなければならないことに気付き
（その時まで、IMAPを並べ替えるだけで良いと、信じて疑わなかった）、確か
にEKBを並べ替えると正しい結果を得ることができました。
バンドに関しての並列化は、まだ最適化の途中で速度的にはあまり高速化で
きていません。８ＰＥ並列で、せいぜい４倍程度です。今後、一層の高速化の
ためのプログラムに対する最適化が必要と思われます。
尚、再三に渡り恐縮ですが、特に並列化に関しては筆
者は素人に毛の生えたレベルなので、記述内容や、考え方、捉え方の間違いや
誤解を指摘してもらえると大変有難いです。（メイルアドレス、
kobayashi.kazuaki-@-nims.go.jp、
"-@-"は変なメイル対策）
[先頭][総目次][最初に戻る][ＶＰＰでの並列化２][ＳＸ４並列化][OpenMPによる並列化][並列化情報]